made data functions more transparent

Author: lukegre

cryogrid_pytools/data/from_earth_engine.py

@@ -1,11 +1,11 @@
 from functools import lru_cache
 
-import xarray as xr
 import ee
-import wxee  # noqa
 import rioxarray  # noqa
-
+import wxee  # noqa
+import xarray as xr
 from loguru import logger
+
 from .utils import _decorator_dataarray_to_bbox
 
 
@@ -27,14 +27,16 @@ def gee_auth():
 
 
 @_decorator_dataarray_to_bbox
-def get_modis_albedo_500m(bbox_WSEN) -> xr.DataArray:
+def get_modis_albedo_500m(bbox_WSEN, res=500, **kwargs) -> xr.DataArray:
     """
     Retrieve MODIS albedo data at 500m resolution for a given bounding box.
 
     Parameters
     ----------
     bbox_WSEN : list or tuple
         Bounding box in the format [west, south, east, north].
+    res : int
+        Resolution in meters. Default is 500m.
 
     Returns
     -------
@@ -69,7 +71,7 @@ def get_modis_albedo_500m(bbox_WSEN) -> xr.DataArray:
 
     # Convert the Earth Engine image to an xarray dataset
     ds = (
-        image_10th_pct.wx.to_xarray(scale=500, region=rio, progress=False)
+        image_10th_pct.wx.to_xarray(scale=res, region=rio, progress=False)
         .pipe(lambda x: x * 0.001)  # Scale factor for albedo
         .isel(time=0, drop=True)["Albedo_BSA_shortwave_p10"]
         .rename("modis_albedo_BSA_shortwave")
@@ -98,14 +100,16 @@ def get_modis_albedo_500m(bbox_WSEN) -> xr.DataArray:
 
 
 @_decorator_dataarray_to_bbox
-def get_aster_ged_100m_v3(bbox_WSEN):
+def get_aster_ged_100m_v3(bbox_WSEN, res=100, **kwargs):
     """
     Retrieve ASTER GED data at 100m resolution for a given bounding box.
 
     Parameters
     ----------
     bbox_WSEN : list or tuple
         Bounding box in the format [west, south, east, north].
+    res : int
+        Resolution in meters. Default is 100m.
 
     Returns
     -------
@@ -132,7 +136,7 @@ def get_aster_ged_100m_v3(bbox_WSEN):
 
     # Convert the Earth Engine image to an xarray dataset
     ds = (
-        image_clipped.wx.to_xarray(scale=100, region=roi, progress=False)
+        image_clipped.wx.to_xarray(scale=res, region=roi, progress=False)
         .isel(time=0, drop=True)
         .assign_attrs(period="2000-2008")
     )
@@ -141,7 +145,7 @@ def get_aster_ged_100m_v3(bbox_WSEN):
 
 
 @_decorator_dataarray_to_bbox
-def get_aster_ged_emmis_elev(bbox_WSEN):
+def get_aster_ged_emmis_elev(bbox_WSEN, res=100, **kwargs):
     """
     Retrieve ASTER GED emissivity and elevation data for a given bounding box.
 
@@ -161,7 +165,7 @@ def get_aster_ged_emmis_elev(bbox_WSEN):
     bbox_WSEN = tuple(bbox_WSEN)
 
     # Retrieve ASTER GED data
-    ds = get_aster_ged_100m_v3(bbox_WSEN)
+    ds = get_aster_ged_100m_v3(bbox_WSEN, res=res)
 
     # Process emissivity bands
     emissivity_bands = [f"emissivity_band{b:02d}" for b in range(10, 15)]

cryogrid_pytools/data/from_planetary_computer.py

@@ -52,7 +52,7 @@ def search_stac_items_planetary_computer(collection, bbox, **kwargs) -> list:
 
 @utils._decorator_dataarray_to_bbox
 def get_dem_copernicus30(
-    bbox_WSEN: list, res_m: int = 30, epsg=32643, smoothing_iters=2, smoothing_size=3
+    bbox_WSEN: list, res: int = 30, epsg=32643, smoothing_iters=2, smoothing_size=3
 ) -> _xr.DataArray:
     """
     Download DEM data from the STAC catalog (default is COP DEM Global 30m).
@@ -61,8 +61,8 @@ def get_dem_copernicus30(
     ----------
     bbox_WSEN : list
         The bounding box of the area of interest in WSEN format.
-    res_m : int
-        The resolution of the DEM data in meters.
+    res : int
+        The resolution of the DEM data in EPSG units.
     epsg : int, optional
         The EPSG code of the projection of the DEM data. Default is
         EPSG:32643 (UTM 43N) for the Pamir region.
@@ -82,13 +82,8 @@ def get_dem_copernicus30(
 
     check_epsg(epsg)
 
-    assert res_m >= 30, (
-        "The resolution must be greater than 30m for the COP DEM Global 30m dataset."
-    )
-    res = res_m / 111_111 if epsg == 4326 else res_m
-
     _logger.info(
-        f"Fetching COP DEM Global data from Planetary Computer (cop-dem-glo-30 @ {res_m}m)"
+        f"Fetching COP DEM Global data from Planetary Computer (cop-dem-glo-30 @ {res:.2g})"
     )
     items = search_stac_items_planetary_computer("cop-dem-glo-30", bbox_WSEN)
     da_dem = _stackstac.stack(
@@ -112,16 +107,16 @@ def get_dem_copernicus30(
 
 
 @utils._decorator_dataarray_to_bbox
-def get_esa_land_cover(bbox_WSEN: tuple, res_m: int = 30, epsg=32643) -> _xr.DataArray:
+def get_esa_land_cover(bbox_WSEN: tuple, res: int = 30, epsg=32643) -> _xr.DataArray:
     """
     Get the ESA World Cover dataset on the target grid and resolution.
 
     Parameters
     ----------
     bbox_WSEN : tuple
         Bounding box in the format (West, South, East, North).
-    res_m : int, optional
-        Resolution in meters. Defaults to 30.
+    res : int, optional
+        Resolution in EPSG units. Defaults to 30 (meters for UTM).
     epsg : int, optional
         EPSG code for the coordinate reference system. Defaults to 32643.
 
@@ -168,10 +163,10 @@ def get_land_cover_classes(item):
     check_epsg(epsg)
 
     # get the units in the projection
-    res = get_res_in_proj_units(res_m, epsg, min_res=10)
+    res = get_res_in_proj_units(res, epsg, min_res=10)
 
     _logger.info(
-        f"Fetching ESA World Cover (v2.0) data from Planetary Computer (esa-worldcover @ {res_m}m)"
+        f"Fetching ESA World Cover (v2.0) data from Planetary Computer (esa-worldcover @ {res}m)"
     )
     items = search_stac_items_planetary_computer(
         collection="esa-worldcover",
@@ -210,7 +205,7 @@ def get_land_cover_classes(item):
 
 
 @utils._decorator_dataarray_to_bbox
-def get_esri_land_cover(bbox_WSEN: tuple, res_m: int = 30, epsg: int = 32643):
+def get_esri_land_cover(bbox_WSEN: tuple, res: int = 30, epsg: int = 32643):
     from .utils import _long_string_processor
 
     items = search_stac_items_planetary_computer(
@@ -222,7 +217,7 @@ def get_esri_land_cover(bbox_WSEN: tuple, res_m: int = 30, epsg: int = 32643):
     da = _stackstac.stack(
         items,
         epsg=epsg,
-        resolution=res_m,
+        resolution=res,
         bounds_latlon=bbox_WSEN,
     ).isel(time=0, band=0)
 
@@ -269,7 +264,7 @@ def get_sentinel2_data(
     bbox_WSEN: tuple,
     years=range(2018, 2025),
     assets=["SCL"],
-    res_m: int = 30,
+    res: int = 30,
     epsg=32643,
     max_cloud_cover=5,
 ) -> _xr.DataArray:
@@ -284,8 +279,8 @@ def get_sentinel2_data(
         Range of years to fetch data for. Defaults to range(2018, 2025).
     assets : list, optional
         List of assets to fetch. Defaults to ['SCL'].
-    res_m : int, optional
-        Resolution in meters. Defaults to 30.
+    res : int, optional
+        Resolution in EPSG units. Defaults to 30 (meters for UTM).
     epsg : int, optional
         EPSG code for the coordinate reference system. Defaults to 32643.
     max_cloud_cover : int, optional
@@ -299,11 +294,9 @@ def get_sentinel2_data(
     from ..utils import drop_coords_without_dim
     from .utils import (
         check_epsg,
-        get_res_in_proj_units,
     )
 
     check_epsg(epsg)
-    res = get_res_in_proj_units(res_m, epsg, min_res=10)
 
     da_list = []
     for year in years:
@@ -348,7 +341,7 @@ def get_sentinel2_data(
 
 @utils._decorator_dataarray_to_bbox
 def get_snow_melt_doy(
-    bbox_WSEN: tuple, years=range(2018, 2025), res_m: int = 30, epsg=32643
+    bbox_WSEN: tuple, years=range(2018, 2025), res: int = 30, epsg=32643
 ) -> _xr.DataArray:
     """
     Calculate the snow melt day of year (DOY) from Sentinel-2 SCL data for a given bounding box and years.
@@ -359,7 +352,7 @@ def get_snow_melt_doy(
         Bounding box coordinates in the format (West, South, East, North).
     years : range, optional
         Range of years to consider. Defaults to range(2018, 2025).
-    res_m : int, optional
+    res : int, optional
         Spatial resolution in meters. Defaults to 30.
     epsg : int, optional
         EPSG code for the coordinate reference system. Defaults to 32643.
@@ -371,7 +364,7 @@ def get_snow_melt_doy(
     """
 
     da = get_sentinel2_data(
-        bbox_WSEN, years=years, res_m=res_m, epsg=epsg, max_cloud_cover=10
+        bbox_WSEN, years=years, res=res, epsg=epsg, max_cloud_cover=10
     )
 
     _logger.info("Calculating snow melt day of year (DOY) from Sentinel-2 SCL data")

cryogrid_pytools/data/utils.py

@@ -26,8 +26,12 @@ def _decorator_dataarray_to_bbox(func):
         The wrapped function with the additional functionality of handling
         DataArray and reprojecting the output.
     """
+    from functools import wraps
+
+    import numpy as np
 
     @_cached
+    @wraps(func)
     def wrapper(*args, **kwargs):
         if len(args) >= 1:
             bbox_or_target = args[0]
@@ -38,6 +42,10 @@ def wrapper(*args, **kwargs):
         if isinstance(bbox_or_target, _xr.DataArray):
             da = bbox_or_target
             bbox = da.rv.get_bbox_latlon()
+            res = np.abs(da.rio.resolution()).mean()
+            epsg = da.rio.crs.to_epsg()
+            kwargs.update(res=res, epsg=epsg)
+
             out = func(bbox, *args, **kwargs)
             out = out.rio.reproject_match(da)
 

cryogrid_pytools/outputs.py

@@ -234,7 +234,7 @@ def read_OUT_regridded_files(
             "Depths are the same for all gridcells. Setting depth as the dimension."
         )
         ds = ds.swap_dims(level="depth")
-        ds = ds.reset_coords("elevation")
+        ds = ds.reset_coords("elevation").astype("float32")
 
     return ds