Intergration With 3rd-Party Extensions#

3rd-party Xarray geospatial extensions may leverage XProj in different ways:

simply consume the API exposed via the “proj” Dataset and DataArray accessors.
register a custom Xarray accessor that implements the proj accessor interface (example below)
implement one or more methods of the proj index interface in a custom Xarray index (example below)

import pyproj
import xarray as xr
import xproj

xr.set_options(display_expand_indexes=True);

CRS-aware Xarray index#

Here below is a basic example of a CRS-aware index, i.e., a custom Xarray index that adds some CRS-dependent functionality (via the proj index interface) on top of Xarray’s default PandasIndex.

Note

The ProjIndexMixin class can be used to mark an Xarray index as formally implementing the proj index interface. However, XProj doesn’t require an Xarray index to explicitly inherit from this mixin class to be recognized as CRS-aware.

import warnings


class GeoIndex(xr.indexes.PandasIndex, xproj.ProjIndexMixin):

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self._crs = None

    def sel(self, *args, **kwargs):
        if self._crs is not None:
            warnings.warn(
                f"make sure that indexer labels have CRS {self._crs}!",
                UserWarning,
            )

        return super().sel(*args, **kwargs)

    @property
    def crs(self):
        return self._crs

    def _proj_set_crs(self, spatial_ref, crs):
        # note: `spatial_ref` is not used here
        print(f"set CRS of index {self!r} to crs={crs}!")

        self._crs = crs
        return self

    def _copy(self, deep=True, memo=None):
        # bug in PandasIndex? crs attribute not copied here
        obj = super()._copy(deep=deep, memo=memo)
        obj._crs = self._crs
        return obj

    def _repr_inline_(self, max_width=70):
        return f"{type(self).__name__} (crs={self._crs})"

    def __repr__(self):
        return f"{type(self).__name__} (crs={self._crs})"

Let’s see it in action by reusing the example dataset above, to which we replace the default indexes of the “lat” and “lon” coordinates with instances of the GeoIndex defined above.

ds_geo_wgs84 = (
    ds_wgs84
    .drop_indexes(["lat", "lon"])
    .set_xindex("lat", GeoIndex)
    .set_xindex("lon", GeoIndex)
)

Note that the CRS of the “lat” and “lon” indexes aren’t yet initialized despite the presence of the “spatial_ref” coordinate:

Mapping the CRS of “spatial_ref” to the “lat” and “lon” geo-indexed coordinates has to be done manually using xarray.Dataset.proj.map_crs():

ds_geo_wgs84 = ds_geo_wgs84.proj.map_crs(spatial_ref=["lat", "lon"])

set CRS of index GeoIndex (crs=None) to crs=epsg:4326!
set CRS of index GeoIndex (crs=None) to crs=epsg:4326!

The CRS of the “lat” and “lon” geo-indexed coordinates is updated via the proj index interface implemented in GeoIndex. Data selection is now CRS-aware! (just a warning is emitted below).

GeoIndex has a crs property (as required by ProjIndexMixin), which is possible to access also via the proj accessor like so:

ds_geo_wgs84.proj("lat").crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Intergration With 3rd-Party Extensions

Contents

Intergration With 3rd-Party Extensions#

CRS-aware Xarray accessor#

CRS-aware Xarray index#

Caveat#