FEAT: Added Equi7Grid

cubo/cubo.py

@@ -25,6 +25,7 @@ def create(
     stac: str = "https://planetarycomputer.microsoft.com/api/stac/v1",
     gee: bool = False,
     stackstac_kw: Optional[dict] = None,
+    e7g: bool = False,
     **kwargs,
 ) -> xr.DataArray:
     """Creates a data cube from a STAC Catalogue as a :code:`xr.DataArray` object.
@@ -78,6 +79,11 @@ def create(
 
         .. versionadded:: 2024.1.0
 
+    e7g : bool
+        Whether to use Equi7Grid instead of UTM.
+
+        .. versionadded:: 2024.6.0
+
     kwargs :
         Additional keyword arguments passed to :code:`pystac_client.Client.search()`.
 
@@ -128,8 +134,8 @@ def create(
             edge_size = (edge_size * getattr(constants, units)) / resolution
 
     # Get the BBox and EPSG
-    bbox_utm, bbox_latlon, utm_coords, epsg = _central_pixel_bbox(
-        lat, lon, edge_size, resolution
+    bbox_proj, bbox_latlon, proj_coords, epsg = _central_pixel_bbox(
+        lat, lon, edge_size, resolution, e7g
     )
 
     # Use Google Earth Engine
@@ -194,7 +200,7 @@ def create(
     else:
 
         # Convert UTM Bbox to a Feature
-        bbox_utm = rasterio.features.bounds(bbox_utm)
+        bbox_proj = rasterio.features.bounds(bbox_proj)
 
         # Open the Catalogue
         CATALOG = pystac_client.Client.open(stac)
@@ -226,7 +232,7 @@ def create(
             items,
             assets=bands,
             resolution=resolution,
-            bounds=bbox_utm,
+            bounds=bbox_proj,
             epsg=epsg,
             **stackstac_kw,
         )
@@ -251,8 +257,8 @@ def create(
         edge_size_m=rounded_edge_size * resolution,
         central_lat=lat,
         central_lon=lon,
-        central_y=utm_coords[1],
-        central_x=utm_coords[0],
+        central_y=proj_coords[1],
+        central_x=proj_coords[0],
         time_coverage_start=start_date,
         time_coverage_end=end_date,
     )

cubo/utils.py

@@ -7,11 +7,73 @@
 from pyproj.database import query_utm_crs_info
 
 
+def _get_epsg(
+    lat: Union[float, int],
+    lon: Union[float, int],
+    edge_size: Union[float, int],
+    resolution: Union[float, int],
+    e7g: bool,
+) -> int:
+    """Gets the EPSG code of the UTM Zone or the Equi7Grid Zone.
+
+    Parameters
+    ----------
+    lat : float
+        Latitude.
+    lon : float
+        Longitude.
+    edge_size : float
+        Buffer distance in meters.
+    resolution : int | float
+        Spatial resolution to use.
+    e7g : bool
+        Whether to use Equi7Grid instead of UTM.
+
+    Returns
+    -------
+    int
+        EPSG code.
+    """
+    if e7g:
+
+        # Try to import ee, otherwise raise an ImportError
+        try:
+            import equi7grid_lite
+        except ImportError:
+            raise ImportError(
+                '"equi7grid-lite" could not be loaded. Please install it, or install "cubo" using "pip install cubo[e7g]"'
+            )
+
+        # EPSG list for Equi7Grid
+        e7g_epsgs = dict(
+            AN=27702, NA=27705, OC=27706, SA=27707, AF=27701, EU=27704, AS=27703
+        )
+
+        # Get the EPSG from latlon
+        grid_system = equi7grid_lite.Equi7Grid(min_grid_size = edge_size * resolution)
+        grid_polygon = grid_system.lonlat2grid(lon, lat)
+        epsg = e7g_epsgs[grid_polygon.zone.values[0]]
+
+    else:
+
+        # Get the UTM EPSG from latlon
+        utm_crs_list = query_utm_crs_info(
+            datum_name="WGS 84",
+            area_of_interest=AreaOfInterest(lon, lat, lon, lat),
+        )
+
+        # Save the CRS
+        epsg = utm_crs_list[0].code
+
+    return int(epsg)
+
+
 def _central_pixel_bbox(
     lat: Union[float, int],
     lon: Union[float, int],
     edge_size: Union[float, int],
     resolution: Union[float, int],
+    e7g: bool,
 ) -> tuple:
     """Creates a Bounding Box (BBox) given a pair of coordinates and a buffer distance.
 
@@ -25,23 +87,16 @@ def _central_pixel_bbox(
         Buffer distance in meters.
     resolution : int | float
         Spatial resolution to use.
-    latlng : bool, default = True
-        Whether to return the BBox as geographic coordinates.
+    e7g : bool
+        Whether to use Equi7Grid instead of UTM.
 
     Returns
     -------
     tuple
-        BBox in UTM coordinates, BBox in latlon, and EPSG.
+        BBox in Equi7Grid coordinates, BBox in latlon, and EPSG.
     """
-    # Get the UTM EPSG from latlon
-    utm_crs_list = query_utm_crs_info(
-        datum_name="WGS 84",
-        area_of_interest=AreaOfInterest(lon, lat, lon, lat),
-    )
-
-    # Save the CRS
-    epsg = utm_crs_list[0].code
-    utm_crs = CRS.from_epsg(epsg)
+    # Get EPSG
+    epsg = _get_epsg(lat, lon, edge_size, resolution, e7g)
 
     # Initialize a transformer to UTM
     transformer = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True)
@@ -52,19 +107,19 @@ def _central_pixel_bbox(
     )
 
     # Transform latlon to UTM
-    utm_coords = transformer.transform(lon, lat)
+    proj_coords = transformer.transform(lon, lat)
 
     # Round the coordinates
-    utm_coords_round = [round(coord / resolution) * resolution for coord in utm_coords]
+    proj_coords_round = [round(coord / resolution) * resolution for coord in proj_coords]
 
     # Buffer size
     buffer = round(edge_size * resolution / 2)
 
     # Create BBox coordinates according to the edge size
-    E = utm_coords_round[0] + buffer
-    W = utm_coords_round[0] - buffer
-    N = utm_coords_round[1] + buffer
-    S = utm_coords_round[1] - buffer
+    E = proj_coords_round[0] + buffer
+    W = proj_coords_round[0] - buffer
+    N = proj_coords_round[1] + buffer
+    S = proj_coords_round[1] - buffer
 
     # Create polygon from BBox coordinates
     polygon = [
@@ -79,7 +134,7 @@ def _central_pixel_bbox(
     polygon_latlon = [list(inverse_transformer.transform(x[0], x[1])) for x in polygon]
 
     # Create UTM BBox
-    bbox_utm = {
+    bbox_proj = {
         "type": "Polygon",
         "coordinates": [polygon],
     }
@@ -90,7 +145,7 @@ def _central_pixel_bbox(
         "coordinates": [polygon_latlon],
     }
 
-    return (bbox_utm, bbox_latlon, utm_coords, int(epsg))
+    return (bbox_proj, bbox_latlon, proj_coords, epsg)
 
 
 def _compute_distance_to_center(da: xr.DataArray) -> np.ndarray: