Interpolation issues on the margins of retreating glaciers

[adehecq]

This is more a glaciological issue than a code-related issue.

The normalized regional hypsometric interpolation is great for filling large gaps. However, for dramatically retreating glaciers, interpolating the elevation change can create artifacts near the margins of the glacier where the elevation change is bounded by the ice thickness. Here is an example over Aletsch glacier from @erikmannerfelt's results with the Swiss TerrA images. The left figure shows the interpolated dh, while the right figure shows the results obtained with topographic maps. One can see that the elevation change is less negative near the margins of the tongue in the right image, which is not reproduced by the interpolation.

An idea could be to look at the relationship "old elevation" vs "new elevation" as opposed to "elevation change" vs "new elevation". Then the elevation change can be calculated as "new elevation" - "interpolated old elevation", which would never be larger than the ice thickness.

[erikmannerfelt]

Hmm, perhaps I'm missing something fundamental, but I don't get how that is different?

If what you mean is simply that the interpolated values should fit the measured ones, that flag exists in the NRHI. But I guess this is a different approach altogether, no?

[adehecq]

Well, at the moment, one calculates new_elevation - old_elevation and aggregate it over altitude bins. But within an altitude bin, elevation change might be relatively homogeneous near the center of the glacier, but smaller near the margins, as seen on the figure above. However, within that bin, the glacier elevation of the old period should be more or less uniform (if the glacier shape did not change too much over time.

A simple example is easier. Let's take a glacier cross-section:

• old elevation could be : 2000, 2001, 2002, 2003, 2002, 2001, 2000 (assuming a rounded glacier shape)
• new elevation could be : 1990, 1981, 1982, 1983, 1982, 1981, 1990 (assuming a thinning of 20 m, except on the margins, where we reach bedrock at 1990 m elevation). Hence dh = 10, 20, 20, 20, 20, 20, 10. With the current approach, one gets mean_dh = 17.5.

With the suggestion, one would take all values within e.g elevation bin 1980-1990 m (in the new elevation), and take the mean. One gets dh_mean = 20. Actually... that does not work, because one would also include the values of the margins from the elevation bin above, which also have smaller values...

I need to think more about that.

[adehecq]

Another idea could be to calculate the mean elevation change (per elevation bin) only for pixels in the most recent outline. Then we linearly interpolate dh values in-between the recent and older outlines.

[rhugonnet]

Fully agree, I think this is what we discussed last time. Linearly interpolating the two outlines, or the dh within the two outline seems like a reasonable solution :)

[erikmannerfelt]

This should work in most cases, but for cases of extreme change and/or splitting of glaciers into multiple bodies, there may be some complications. Maybe we can work it to fix those as well, just saying that there are some difficult ones!

[erikmannerfelt]

Further thought: There should be some threshold whereby linear interpolation should not be made. For example, I sometimes work with 80% area reductions, meaning 80% of the glacier by area would be linearly interpolated instead of hypsometrically. If there was a threshold so only glaciers without extreme change were affected, this would be a good approach!

Also, based on the above statement: A more pessimistic view is that your suggested approach only works on glaciers with limited extent change

[rhugonnet]

Agreed. Summarized, I think we have to keep in mind that this problem is mostly exacerbated for:

• Some types of imagery that has gaps in ablation area (for most those are not an issue to resolve even for moderate-resolution photogrammetry), such as terrestrial photogrammetry or incomplete historical imagery,
• For when there is a lot of outliers in the data and one wishes to filter the dh observations per elevation band using, let's say, 3NMAD or 5NMAD from the median or mean dh. This will artificially create gaps in the margins where the dh is ~0 and fill them with too large dh values. So the filtering can have a huge impact.

Some work with simulated gaps should help us contrain this. When are we doing this "best interpolation" paper again? :stuck_out_tongue:

[adehecq]

Further thought: There should be some threshold whereby linear interpolation should not be made. For example, I sometimes work with 80% area reductions, meaning 80% of the glacier by area would be linearly interpolated instead of hypsometrically. If there was a threshold so only glaciers without extreme change were affected, this would be a good approach!

Well if you have an 80% retreat and no observation of it, then you can only guess... Let's keep in mind with this approach that this is only meant to interpolate gaps, so the issue doesn't arise where there are observations. But it was my mistake. What I meant with "linearly interpolate" is rather a standard spatial interpolation (e.g. bilinear or IDW), so if you have lots of observations to constrain your interpolation, that should not be a problem.

[adehecq]

For when there is a lot of outliers in the data and one wishes to filter the dh observations per elevation band using, let's say, 3NMAD or 5NMAD from the median or mean dh. This will artificially create gaps in the margins where the dh is ~0 and fill them with too large dh values. So the filtering can have a huge impact.

Good point for filtering too. One should filter pixels inside/outside of the most recent outline differently.

[adehecq]

Some more thoughts on the problem from this Berthier et al. (2010) paper and an associated presentation slide below.

[erikmannerfelt]

Well if you have an 80% retreat and no observation of it, then you can only guess... Let's keep in mind with this approach that this is only meant to interpolate gaps, so the issue doesn't arise where there are observations. But it was my mistake. What I meant with "linearly interpolate" is rather a standard spatial interpolation (e.g. bilinear or IDW), so if you have lots of observations to constrain your interpolation, that should not be a problem.

I don't mean when we have no observations. I mean when we have sparse observations with extreme outline differences.

I'm intrigued by the idea, and I'm sure it may work in many cases! Just not in all cases (which no approach does anyway).

[erikmannerfelt]

Some more thoughts on the problem from this Berthier et al. (2010) paper and an associated presentation slide below.

I was playing with that idea for the TerrA project actually! Basically interpolating polynomials perpendicular to glacier flow. One just needs a centreline, and a fallback approach for the top of the glacier where things can get messy. But it's definitely interesting to pursue!

[adehecq]

I'm intrigued by the idea, and I'm sure it may work in many cases! Just not in all cases (which no approach does anyway).

Well, of course it won't be perfect in very difficult cases, but it will most likely still work better than the hypsometric approach, which assume that thinning is purely dependent on the elevation!

[erikmannerfelt]

I'm just brainstorming now, but one could perhaps get that effect by combining the general hypsometric signal with a proximity distance to the edge. Basically a two-variable regression instead of one.

[rhugonnet]

Hypsometric + Proximity should certainly work. However only if you account for changes in proximity over time. So one always need to interpolate between the old and new outlines!