Trouble with Applying Glacier Mask in BlockwiseCoreg for DEM Coregistration

[cwywilson]

Hi Everyone,

I'm trying to use blockwise_coreg = xdem.coreg.BlockwiseCoreg(xdem.coreg.NuthKaab(), subdivision=16) to coregister my photogrammetry-derived DEM (which contains a lot of noise) to a reference DEM with a glacier mask.

However, after several attempts, it seems that BlockwiseCoreg doesn't properly apply the glacier mask during the fit and apply steps, and it's still performing coregistration on the glacier-covered areas.

The glacier mask works correctly when using xdem.coreg.NuthKaab() as a global coregistration.

Also, whenever I use BlockwiseCoreg, it takes up a lot of memory. Could this be a memory limitation affecting the process size? Is there an alternative way to handle this, such as using Dask to manage the memory load?

Could you help me understand how to fix this issue?

[rhugonnet]

Thanks @cwywilson for opening this! Might be linked to #604 that @quinto-clara is looking at right now.
Could you share a short code example, and print the content of your blockwise_coreg._meta["step_meta"] and blockwise_coreg.stats() (which has the number of inliers per block), so that we can reproduce this?

The blockwise functionality is undergoing a bit of a "revamp" as this code was on the older side, and less tested than the rest. It should become more robust in the future 😉.
We're also working on implementing Dask support for this type of method (#525), but it's not there yet.
For reducing current memory usage, maybe defining the global parameters linked to the parallelization of the fit() through concurrent.futures.ThreadPoolExecutor() would solve the issue... Any advice on this @erikmannerfelt?

[cwywilson]

Hi @rhugonnet,

Apologies for the delay in getting back to you. The blockwise code takes a significant amount of time to run, and due to the large memory requirements, it often crashes at the end.

However, I noticed an issue: the blockwise_coreg.stats().inlier_count is higher than np.count_nonzero(inlier_mask), and I’m unsure why this discrepancy occurs. I’d appreciate any insights you might have. Here’s an overview of the code I’m using:

I have refdem, demtobealign, and glacier_outline defined as:
glacier_outlines = gu.Vector("OUTLINE.shp")
I also have a slope mask:

slope = xdem.terrain.slope(reference_dem)

# Create masks for slopes between 0-5 degrees and slopes over 30 degrees
slope_mask_low = (slope.data <= 5)  # Mask for slopes between 0-5 degrees
slope_mask_high = (slope.data >= 30)  # Mask for slopes over 30 degrees

# Combine the masks (0-5 degrees or >= 30 degrees)
slope_combined_mask = slope_mask_low | slope_mask_high

The inlier mask is defined as:

glacier_mask = glacier_outlines.create_mask(reference_dem)
combined_mask = slope_combined_mask | glacier_mask.data
inlier_mask = ~combined_mask

When I count the inliers using np.count_nonzero(inlier_mask), I get 84,965,169. However, after running blockwise_coreg, the reported inlier_count is 105,454,194.

If needed, I can upload both DEMs and the notebook for you to review.

Thank you for your help!

Wilson

[rhugonnet]

Thanks @cwywilson,

EDIT: After more investigation, I found the bug! Indeed the mask was ignored 😄
The inlier_mask variable is used twice for two different purposes in this function, overwriting the previous one:

xdem/xdem/coreg/base.py

Line 2964 in 5bde36e

inlier_mask = groups == i

It's a 2-line fix!

For memory usage: Until we finalize Dask support of this, the best would likely be for you to manually split your DEM into tiles (with some overlap), and save those as smaller files. It will then be extremely fast to co-register those with a simple Nuth and Kääb, and later merge them back together (for example with: https://geoutils.readthedocs.io/en/stable/gen_modules/geoutils.raster.merge_rasters.html#geoutils.raster.merge_rasters) 😉.

[rhugonnet]

Opened #624 and resolved in #625.
It'll take a bit of time (probably at least a month) until we can publish a new release, because we also have to fix #584 for BlockwiseCoreg (this was not an issue for you running the code on the published xDEM versions, it's only a recent development issue on the latest GitHub main due to other changes).

[rhugonnet]

Thanks again for reporting this @cwywilson, good catch!
As I mentioned above, the BlockwiseCoreg was an old functionality with not so many tests. We are improving it to be more modular, robust (this will come soon), and computationally efficient (this will take a while longer but is on the way too).

[cwywilson]

I’m glad we could identify the issues—thank you for your efforts. XDEM is indeed a powerful tool and has great potential, as demonstrated in this study: https://tc.copernicus.org/articles/18/3195/2024/, as well as in historical image reconstruction. I would love to contribute to the XDEM project in the future. However, as I’m in the final year of my PhD, I need to focus on completing my current work before taking on new tasks.