Ideas for xDEM structure

Initial global structure of package:

Editable/commentable document here: https://docs.google.com/presentation/d/1kqa4QYZHIikIJMtLien_eWq33GZ7CtsDz3Au67STumA/edit?usp=sharing

coreg.py

Each coregistration function should have the following attributes (names to be discussed):

• fit: calculates the transformation between a reference DEM and a source DEM (inputs: dem_ref, dem_tba, kwargs - Outputs: save result as transformation matrix + optional ramp)
• transform_matrix: the calculated transformation matrix (except for deg >1 deramping)
• apply_to_dem: apply the previously calculated transformation to a given DEM (input: dem_tba - output: coregistered DEM)
• apply_to_points: apply the transformation to a point/camera, i.e. X/Y/Z and possibly orientation (input: camera position/rotation=array of size 3 or 12 - output: transformed camera)

Ideally, the coregistration methods should have a method to chain/add them. Below is an example of a typical use case:

# Calculate a horizontal shift
coreg_nk = xdem.coreg.Coregistration("nuth_kaab")
coreg_nk.fit(dem_ref, dem_tba)
dem_coreg = coreg_nk.applyt_to_dem(dem_tba)

# Remove a tilt
coreg_ramp = xdem.coreg.Coregistration("deramping")
coreg_ramp.fit(dem_ref, dem_coreg)

# Sum up the two transformations
coreg_final = coreg_nk + coreg_ramp

# Apply transformation to all similar DEMs
for dem in list_dems:
    dem_coreg = coreg_final.apply_to_dem(dem)

# Apply transformation to a camera
camera = (X, Y, Z)
camera_coreg = coreg_final.apply_to_points(camera)

filters.py

Contains a set filter operations to 1D and 2D arrays. These belongs to mostly two categories: outlier filtering and interpolation.

2D filter

To be used directly on array/raster objects. Input: 2D np array or xarray Output: 2D np array or xarray

Filtering

• median filter: kwargs = kernel size, threshold
• NMAD filter: kwargs = threshold
• others?

Interpolation

• bilinear interpolation (use rasterio's?)
• higher degree interpolation (polynomial, spline...)
• krigging: requires spstats functionalities

1D filter

To be used on the output of the hypsometric method for example. Input: 1D np array or xarray Output: 1D np array or xarray

Filtering

• median filter
• NMAD filter
• other? Romain?

Interpolation

• linear
• polynomial
• LOWESS
• bayesian (based on e.g. regional hypsometric results)

dDEM / DEMCollection

dDEM could be a wrapper around some functionalities described above. Here some ideas of use case:

# Initialize object (load attributes etc)
dDEM = xdem.dDEM(dem1, dem2)

# Reproject (or coregister?) both DEM to a common grid and calculate difference in self.dh (?)
dDEM.reproject_common(ref_grid)
dDEM.subtract()

# Reject all pixels whose value differ by more than 3 NMAD
dDEM.filter("NMAD", factor=3)

# Remove all pixels that differ by more than 50 m from their neighbors
dDEM.filter("median", kernel_size = (5,5), threshold=50)

# Fill voids less than 10 pixels large
dDEM.interpolate("bilinear", max_size=10)

# Calculate median dh in 100 m bins for individual glaciers using the local hypsometric approach
# Bins with less than 50% coverage are excluded
# implace forces to return an object rather than replacing the dh values directly
hypso_output = dDEM.interpolate("local_hypsometric", outlines=rgi_outlines, bin_height=100, coverage_thresh =0.5, inplace=False)

# Interpolate missing bins
hypso_output.interpolate("polynomial")

# Calculate geodetic volume change for each glacier
glaciers_volume_change = hypso_output.get_dv()

# Use to fill dh voids, e.g. to visualize
dh_filled = hypso_output.as_raster(ref_grid)
dDEM.filed_data = dh_filled

# Fill larger gaps using the regional hypsometric approach
dDEM.interpolate("regional_hypsometric", outlines=rgi_outlines, bin_height=100)

# Calculate regional geodetic volume change
regional_volume_change = dDEM.get_dv()