Volume statistics

[fmaussion]

Hi! Thanks for the dataset - we are looking into making this the new reference for GlacierMIP4 since this is the only available dataset for RGI7.

I have a few questions regarding total volume statistics

Glacier per glacier

I'm looking at

iceboost/iceboost_deploy.py

Lines 918 to 920 in 118ec4f

	# Calculate volume from produced data points
	dataID = data.loc[data["polygon"] == glacierID]
	f = 0.001 * areaID / len(dataID)

And trying to make sense of how the total glacier volume is computed on the file attributes. In Farinotti et al., 2019, the computational domains are the single glacier entities and the total volume of a glacier would be:

thickness.sum() * dx**2  # with dx the spatial resolution

This is done that way because this is how the individual models worked out the volume (in OGGM we would compute the flowline volume first and then distribute it on the grid ensuring that mass is conserved despite of the area mismatch.

It's correct that the total glaciated area in the gridded product might be different from the actual (vector) area of the glacier. I suppose that this is what the division in your code above is trying to achieve? But if I try to reproduce this myself I get slightly different results.

Take RGI60-06.00002 for example (I just picked one at random).

The total (naive) volume would be:

ds = rioxr.open_rasterio('../ice_thickness/iceboost_v2/RGI62/rgi6/RGI60-06.00002.tif')
# Take the thickness, multiply by the grid point area (here 100^2), then sum it and convert to km3
vol_naive = (ds.isel(band=0).astype(np.float64) * ds.attrs['resX']**2).sum()* 1e-9  
vol_naive

Which yields 0.07769982 km3.

If I try to apply a linear correction factor for area mismatch (which, btw, is quite an assumption since the correction also applies to the thickest areas of the glacier while the errors are more likely to be spread at the borders, but this is not too important here):

area_rgi = ds.attrs['area'] # Area RGI (km2) - here 1.8970186320105116, in the RGI file 1.897, so close enough
area_grid = (~ds.isel(band=0).isnull()).sum()*ds.attrs['resX']**2*1e-6 # Area grid (km2) - here 2.37
vol_corr = vol_naive * area_rgi/area_grid
vol_corr

Which yields 0.06219325 km3.

However, the attributes of the file say:

ds.attrs['volume']

Which yields 0.06530186680162677 km3, so slightly different from the above still. What am I doing wrong? What is f = 0.001 * areaID / len(dataID) doing in the code above?

Finally, I wonder if the 0.001 is not forgetting that some glaciers have finer grid resolutions?

Regional / glacier complex volumes

I think it would be really really beneficial to have glacier complex products for RGI7 as well if you can. As of now, models will have to stitch the files together, which would be a pain to do an likely to create errors for some.

Still, this leads me to the other question: the individual glacier volumes computed above are a simple "by product", while the actual volume should always be the one computed on the original grid (the glacier complex one). Did you check that the sum of the parts is the total volume, or is this what the code above is doing?

Thanks!

[fmaussion]

To add some context to what difference it makes globally (removing the CL2 glaciers in Greenland):

Area:

• total RGI6 area: 705738 km2
• total gridded area in Iceboost v2: 767361 km2 (+9%)

Volume:

• total volume as per the file attributes: 149847 km3
• total volume (sum of grid points): 154143 km3 (+3%)
• total volume (linear area correction method): 147236 km3 (-1.5%)

Looking at the geotiff dataset it looks like the choice of the grid was conservative (every pixel touched is providing ice thickness), which is fine as a choice but needs to be documented somewhere perhaps.

OGGM thicknesses:

Iceboost thicknesses:

[nmaffe]

1. Point generation

The points are generated (the model runs on points, not on a grid) in generate_points() inside utils_metadata.py. At least 1000 points need to be generated at 100x100m in any glacier. If not, the spacing is decreased at 1-meter step until the condition is met. The generated points need to be inside the RGI polygons (and outside nunataks). This set of generated points, where the model is queried, is used to calculate the volumes (rather than the grid).

In magenta the generated points used to query the model (and to calculate the volume).

2. Areas as tif attributes

The glacier areas saved in the tif attributes are calculated using pyproj.geod.Geod inside the add_regional_features() method.
area, perimeter = geod.geometry_area_perimeter(geometry) and
area = abs(area) * 1e-6, in km2.
It is these areas that are used to compute glacier volumes (and saved as attributes, 1.897018632 in your example), rather than using the grid areas from dx**2, which suffers from border effects (2.37 in your example).

I expect these areas in the attributes to be very close to RGI (1.897 in your example), if the latter uses a method similar to geod.geometry_area_perimeter(geometry). The motivation of using this method instead importing RGI areas is to make the code accept any polygon.

3. Individual glacier volumes as tif attributes

As the model works on points rather than on a grid natively, glacier volumes are calculated by:
f = 0.001 * areaID / len(dataID)
volID = dataID["thickness"].sum() * f
The area of each pixel is the glacier area areaID (the one described above) divided by number of generated points inside the glacier: areaID / len(dataID). Essentially the procedure is a Riemann sum in 2d.
The factor 0.001 converts thickness meters (output of the model ) to kilometers (to have the volume in km3).

Finally, I wonder if the 0.001 is not forgetting that some glaciers have finer grid resolutions?

No. In L 917 the grid resolution comes inside len(dataID), which is the number of generated points, so the finer grid is accounted for.

4. Geotiff dataset and Regional / glacier complex volumes

Looking at the geotiff dataset it looks like the choice of the grid was conservative (every pixel touched is providing ice thickness)

Correct, no pixel is removed: line 911, all_touched=True.
The generated grid is produced not to have empty spaces inside the polygons (by nearest neighbor interpolation), but maybe it can be made fancier. I get that having the areas and the volumes as attributes rather than having those calculated from the tif directly is not optimal.

A glacier complex can be achieved by:
from rioxarray import merge
glacier_complex = merge.merge_arrays(list_of_individual_glacier_rasters, method='max')
I can produce complexes within some computational limits. I'm thinking for example RGI13 there are probably too many glaciers connected to produce one single complex, but I'll try.

Did you check that the sum of the parts is the total volume, or is this what the code above is doing?

Not sure if I understood correctly - the correct volumes (and areas) are those in the attributes. The mismatch is due to the area of the grid which is not the same as the generated points (see the figure above), and also because of the policy of keeping all touched pixels.

[fmaussion]

Thanks for the reply - fine by me, but this may need documented somewhere if possible!

The point regarding RGI7C is important. Users using gridded models will be quite keen to have access to thickness grids on the RGI7C product as provided here: https://www.glims.org/rgi_user_guide/products/glacier_complex_product.html.

Did you check that the sum of the parts is the total volume, or is this what the code above is doing?

This was referring to RGI7C. Say a glacier complex is made of RGI7G glaciers 1, 2, and 3, then you would expect the final RGI7C product to be the sum of the volume of the three sub-glaciers. However given your replies above and that final product users cannot recompute the "true" volumes themselves this is perhaps not as important as I thought.