Skip to content

Add dem Xarray accessor mirroring DEM class#656

Open
rhugonnet wants to merge 13 commits intoGlacioHack:mainfrom
rhugonnet:add_xr_accessor
Open

Add dem Xarray accessor mirroring DEM class#656
rhugonnet wants to merge 13 commits intoGlacioHack:mainfrom
rhugonnet:add_xr_accessor

Conversation

@rhugonnet
Copy link
Member

@rhugonnet rhugonnet commented Nov 21, 2024

This PR adds the dem Xarray accessor mirroring the DEM class, since the finalization of the rst accessor in GlacioHack/geoutils#446, and Dask/Multiproc support for terrain attributes and vertical CRS transformation.

More complex implementations (coregistration, uncertainty analysis) will be done in separate PRs.

The accessor allows to access all attributes and run all methods already implemented for DEMs from a xarray.DataArray object (e.g., ds.dem.slope(), ds.dem.coregister_3d()), including functionalities from rasters inherited from the rst accessor (e.g., ds.dem.reproject()).
The setup follows a double inheritance between the RasterAccessor and DEMBase, so that the dem accessor can also run reproject(), etc.

Resolves #392
Resolves #370
Resolves #933
Resolves #934
Resolves #935

Summary of changes

To understand the changes below, first see changes related to the rst accessor in GeoUtils: GlacioHack/geoutils#446

As for the rst accessor, most of the contents of the DEM class are moved into a DEMBase parent class, inherited by both the DEM class and dem accessor class.

In this case, however, inheritance is more delicate because DEM also needs to inherit functionalities from Raster (that are not in RasterBase) and the dem accessor needs to inherit functionalities from rst (that are not in RasterBase either).
Therefore, the DEM class has dual-inheritance of Raster and DEMBase, while the dem accessor has dual inheritance of rst and DEMBase. This makes things slightly more complex when explicitly calling super(), which is only needed in a few occurences of core functions or __init__, but otherwise easy for other functionalities!

Similar as for the GeoUtils PR, the from_array functions allow to return a same-class object. In this case, only DEM needs to override with its own from_array (as dem already returns a xarray.DataArray through rst).

TO-DO

Code

  • Ensure dual support of masked arrays and NaN arrays for all DEM methods,
  • Add Dask support for all terrain attributes,
  • Add Dask/Multiprocessing support for to_vcrs,
  • Add tests specific to DEMBase,
  • Add tests specific to DEMAccessor,
  • Add tests specific to chunked implementation in terrain and vcrs modules.

Other Dask support to add

Mostly left is to add Dask support to coregistration and uncertainty analysis (both heteroscedasticity and variography/autocorrelation). This is related to old #525, and will be facilitated by ongoing #759

@rhugonnet rhugonnet marked this pull request as draft November 21, 2024 00:58
@rhugonnet
Copy link
Member Author

rhugonnet commented Mar 12, 2026

@belletva @marinebcht This PR might remove the short info() we added not too long ago, I'm curious of your opinion.

As having a separate volatile vcrs attribute had quite a big impact on how to re-build an array (copy() and from_array() methods, which are duplicated in the case of an accessor), I used the opportunity of this PR to reconcile vcrs and crs.
Now, the vcrs can still be passed by a user (on opening, or set later), but it simply updates the 3rd dimension of the crs, as was discussed in #370.

This means that we don't need a specific copy(), from_array() (and potentially info()) method in xDEM; it all links back to GeoUtils naturally.

To give you an idea of how it looks for info(), it simply matches the first if loop that @marinebcht had previously written in case the CRS was 3D, which used directly pyproj.CRS.name, for example for a 3D CRS:

  • "['Horizontal: ETRS89 / UTM zone 33N; Vertical: EGM96 height']". (we actually need to update GeoUtils that print to_string() which is a bit too long; we should just show the name)

The "problem" is that, if no vertical CRS is defined, the CRS simply displays as a classic CRS (no "horizontal" and "vertical"). If we think that is too implicit, I could be convinced of keeping our info() override in DEMBase to always show "Horizontal: ... ; Vertical: ..." (and print "None" when no vertical reference exists, for instance).

What do you think?

@rhugonnet rhugonnet marked this pull request as ready for review March 12, 2026 23:07
@rhugonnet
Copy link
Member Author

@belletva @adehecq @adebardo @marinebcht This PR is pretty much done: All tests passing locally. You can also review, and we'll be able to merge once a new GeoUtils is released.

To avoid big conflicts, it'll be better to merge this first, then the coregistration/uncertainty PR #759, then change the internal structure of some Coreg functions to support Dask/Multiprocessing more easily.

@marinebcht
Copy link
Contributor

marinebcht commented Mar 17, 2026

Yes, you right, as the info() function in xDEM is only used to add a layer to handle .vcrs when they are (were!) separate from .crs objects, it is not really usefull now.

Plus, the idea to correct the error reported in GlacioHack/geoutils#805 (dataset 2) is to move the componed crs treatment into geoutils to handle CRSs when have directly a vertical axis.

My only suggestion s about your last point (crs with no vertical axis that I have mention in #805), either we print classical EPSG:XXXX for exemple or we print [ EPSG:XXXX, None ]. Perhaps, we could add a xxx.crs_to_string() function (and call directly it in infos()) in geoutils without the explicit ; NONE], and with it in xDEM, since we are dealing with actual elevation data. The only parameter I see is if the use of the VCRS become the norm in the futur.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

2 participants