add chili, proxy for insolation

easysnowdata/topography.py

@@ -7,7 +7,7 @@
 import pystac_client
 import planetary_computer
 import os
-
+import ee
 import odc.stac
 odc.stac.configure_rio(cloud_defaults=True)  
 
@@ -68,52 +68,126 @@ def get_copernicus_dem(bbox_input: gpd.GeoDataFrame | tuple | shapely.geometry.b
 
     return cop_dem_da
 
-# find a way to say given bounding box, best and available DEMs?
-#https://github.com/OpenTopography/OT_3DEP_Workflows
-#https://github.com/OpenTopography/OT_BulkAccess_COGs/blob/main/OT_BulkAccessCOGs.ipynb
-
-def get_3dep_dem(bbox_input: gpd.GeoDataFrame | tuple | shapely.geometry.base.BaseGeometry | None = None,
-                 dem_type: str = 'DSM',
+def get_chili(bbox_input: gpd.GeoDataFrame | tuple | shapely.geometry.base.BaseGeometry | None = None,
 ) -> xr.DataArray:
     """
-    Fetches 3DEP DEM data from Microsoft Planetary Computer.
+    Fetches Continuous Heat-Insolation Load Index (CHILI) data for a given bounding box. 
 
-    This function retrieves the 3D Elevation Program (3DEP) Digital Elevation Model (DEM) data
-    for a specified bounding box and DEM type.
+    Description:
+    CHILI is a topographic index that quantifies the combined effect of solar radiation and surface temperature.
+    It is derived from the ALOS World 3D - 30m (AW3D30) dataset and is available globally between 70°N and 70°S.
+    The values range from 0 to 1, with classifications warm (0.767,1], cool [0,0.448), and neutral [0.448,0.767].
 
     Parameters
     ----------
-    bbox_input : geopandas.GeoDataFrame or tuple or Shapely Geometry
+    bbox_input : geopandas.GeoDataFrame or tuple or shapely.Geometry
         GeoDataFrame containing the bounding box, or a tuple of (xmin, ymin, xmax, ymax), or a Shapely geometry.
-    dem_type : str, optional
-        The DEM type, either 'DSM' (Digital Surface Model) or 'DTM' (Digital Terrain Model). Default is 'DSM'.
 
     Returns
     -------
-    pystac_client.ItemCollection
-        A STAC ItemCollection containing the search results for the 3DEP DEM data.
-
-    Raises
-    ------
-    ValueError
-        If the dem_type is not 'DSM' or 'DTM'.
+    xarray.DataArray
+        CHILI DataArray for the specified region.
+
+    Examples
+    --------
+    >>> import geopandas as gpd
+    >>> easysnowdata
+    >>> 
+    >>> # Define a bounding box for an area of interest
+    >>> bbox = (-122.5, 47.0, -121.5, 48.0)
+    >>> 
+    >>> # Fetch CHILI data
+    >>> chili_data = easysnowdata.remote_sensing.get_chili(bbox)
+    >>> 
+    >>> # Plot the data
+    >>> chili_data.plot(cmap='viridis')
 
     Notes
     -----
-    This function currently returns the search results rather than the actual DEM data.
-    Further processing would be needed to load and process the DEM data from the search results.
+    - CHILI data is only available for latitudes between 70°N and 70°S.
+    - The function uses Google Earth Engine to access the dataset.
 
-    The 3DEP program provides high-quality elevation data for the United States.
+    Data citation:
+    Theobald, D.M., Harrison-Atlas, D., Monahan, W.B., Albano, C.M. (2015).
+    Ecologically-Relevant Maps of Landforms and Physiographic Diversity for Climate Adaptation Planning.
+    PLoS ONE 10(12): e0143619. https://doi.org/10.1371/journal.pone.0143619
     """
-    if dem_type != 'DSM' and dem_type != 'DTM':
-        raise ValueError("3DEP DEM type is available as DSM and DTM. Please select either DSM or DTM.")
+
+    # Initialize Earth Engine with high-volume endpoint
+    ee.Initialize(opt_url='https://earthengine-highvolume.googleapis.com')
 
     # Convert the input to a GeoDataFrame if it's not already one
-    bbox_gdf =  convert_bbox_to_geodataframe(bbox_input)
+    bbox_gdf = convert_bbox_to_geodataframe(bbox_input)
+    # Access the CHILI dataset
+    image = ee.Image("CSP/ERGo/1_0/Global/ALOS_CHILI")
+    image_collection = ee.ImageCollection(image)
 
-    catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",modifier=planetary_computer.sign_inplace)
-    search = catalog.search(collections=[f"3dep-lidar-{dem_type.lower()}"],bbox=bbox_gdf.total_bounds)
-    #dep_dem_da = odc.stac.load(search.items(),bbox=bbox_gdf.total_bounds,chunks={})
-    #dep_dem_da = dep_dem_da.rio.write_nodata(-32767,encoded=True)
+    # Get projection information
+    crs = image.projection().getInfo()['crs']
+    transform = image.projection().getInfo()['transform']
+
+    # Load the data using xarray and Earth Engine
+    chili_da = xr.open_dataset(image_collection, engine='ee', geometry=tuple(bbox_gdf.total_bounds), projection=ee.Projection(crs=crs,transform=transform))['constant'].rio.set_spatial_dims(x_dim='lon', y_dim='lat')
+
+    # Clip the data to the bounding box
+    chili_da = chili_da.rio.clip_box(*bbox_gdf.total_bounds, crs=bbox_gdf.crs).squeeze()
+
+    # Check if data is available for the specified region
+    if chili_da.dropna(dim='lat',how='all').dropna(dim='lon',how='all').shape == (0,0):
+        print("No CHILI data available for this location. CHILI data is only available for latitudes between 70°N and 70°S.")
+
+    chili_da = (chili_da - chili_da.min()) / (chili_da.max() - chili_da.min())
+
+    chili_da.attrs['data_citation'] = "Theobald, D.M., Harrison-Atlas, D., Monahan, W.B., Albano, C.M. (2015). Ecologically-Relevant Maps of Landforms and Physiographic Diversity for Climate Adaptation Planning. PLoS ONE 10(12): e0143619. https://doi.org/10.1371/journal.pone.0143619"
+
+    return chili_da
+
+# find a way to say given bounding box, best and available DEMs?
+#https://github.com/OpenTopography/OT_3DEP_Workflows
+#https://github.com/OpenTopography/OT_BulkAccess_COGs/blob/main/OT_BulkAccessCOGs.ipynb
 
-    return search
+# def get_3dep_dem(bbox_input: gpd.GeoDataFrame | tuple | shapely.geometry.base.BaseGeometry | None = None,
+#                  dem_type: str = 'DSM',
+# ) -> xr.DataArray:
+#     """
+#     Fetches 3DEP DEM data from Microsoft Planetary Computer.
+
+#     This function retrieves the 3D Elevation Program (3DEP) Digital Elevation Model (DEM) data
+#     for a specified bounding box and DEM type.
+
+#     Parameters
+#     ----------
+#     bbox_input : geopandas.GeoDataFrame or tuple or Shapely Geometry
+#         GeoDataFrame containing the bounding box, or a tuple of (xmin, ymin, xmax, ymax), or a Shapely geometry.
+#     dem_type : str, optional
+#         The DEM type, either 'DSM' (Digital Surface Model) or 'DTM' (Digital Terrain Model). Default is 'DSM'.
+
+#     Returns
+#     -------
+#     pystac_client.ItemCollection
+#         A STAC ItemCollection containing the search results for the 3DEP DEM data.
+
+#     Raises
+#     ------
+#     ValueError
+#         If the dem_type is not 'DSM' or 'DTM'.
+
+#     Notes
+#     -----
+#     This function currently returns the search results rather than the actual DEM data.
+#     Further processing would be needed to load and process the DEM data from the search results.
+
+#     The 3DEP program provides high-quality elevation data for the United States.
+#     """
+#     if dem_type != 'DSM' and dem_type != 'DTM':
+#         raise ValueError("3DEP DEM type is available as DSM and DTM. Please select either DSM or DTM.")
+
+#     # Convert the input to a GeoDataFrame if it's not already one
+#     bbox_gdf =  convert_bbox_to_geodataframe(bbox_input)
+
+#     catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",modifier=planetary_computer.sign_inplace)
+#     search = catalog.search(collections=[f"3dep-lidar-{dem_type.lower()}"],bbox=bbox_gdf.total_bounds)
+#     #dep_dem_da = odc.stac.load(search.items(),bbox=bbox_gdf.total_bounds,chunks={})
+#     #dep_dem_da = dep_dem_da.rio.write_nodata(-32767,encoded=True)
+
+#     return search