Merge pull request #11 from ESDS-Leipzig/edge_size_m

Units for edge_size

README.md

@@ -172,6 +172,41 @@ da = cubo.create(
 )
 ```
 
+### Using different units for `edge_size`
+
+By default, the units of `edge_size` are pixels. But you can modify this using the `units` argument:
+
+```python
+da = cubo.create(
+    lat=4.31,
+    lon=-76.2,
+    collection="sentinel-2-l2a",
+    bands=["B02","B03","B04"],
+    start_date="2021-06-01",
+    end_date="2021-06-10",
+    edge_size=1500,
+    units="m",
+    resolution=10,
+)
+```
+
+> [!TIP]
+> You can use "px" (pixels), "m" (meters), or any unit available in [`scipy.constants`](https://docs.scipy.org/doc/scipy/reference/constants.html#units).
+
+```python
+da = cubo.create(
+    lat=4.31,
+    lon=-76.2,
+    collection="sentinel-2-l2a",
+    bands=["B02","B03","B04"],
+    start_date="2021-06-01",
+    end_date="2021-06-10",
+    edge_size=1.5,
+    units="kilo",
+    resolution=10,
+)
+```
+
 ### Using another endpoint
 
 By default, `cubo` uses Planetary Computer. But you can use another STAC provider endpoint if you want:

cubo/__init__.py

@@ -1,6 +1,6 @@
 """cubo - On-Demand Earth System Data Cubes in Python"""
 
-__version__ = "2024.1.0"
+__version__ = "2024.1.1"
 __author__ = "David Montero Loaiza <[email protected]>"
 __all__ = []
 

cubo/cubo.py

@@ -7,6 +7,7 @@
 import rasterio.features
 import stackstac
 import xarray as xr
+from scipy import constants
 
 from .utils import _central_pixel_bbox, _compute_distance_to_center
 
@@ -19,6 +20,7 @@ def create(
     end_date: str,
     bands: Optional[Union[str, List[str]]] = None,
     edge_size: Union[float, int] = 128.0,
+    units: str = "px",
     resolution: Union[float, int] = 10.0,
     stac: str = "https://planetarycomputer.microsoft.com/api/stac/v1",
     gee: bool = False,
@@ -45,11 +47,23 @@ def create(
     bands : str | List[str], default = None
         Name of the band(s) from the collection to use.
     edge_size : float | int, default = 128
-        Size of the edge of the cube in pixels. All edges share the same size.
+        Size of the edge of the cube in the units specified by :code:`units`. All edges share the same size.
 
         .. warning::
            If :code:`edge_size` is not a multiple of 2, it will be rounded.
 
+    units : str, default = 'px'
+        Units of the provided edge size in :code:`edge_size`. Must be 'px' (pixels), 'm' (meters), or a unit
+        name in https://docs.scipy.org/doc/scipy/reference/constants.html#units.
+
+        .. versionadded:: 2024.1.1
+
+        .. warning::
+           Note that when :code:`units!='px'` the edge size will be transformed to meters (if :code:`units!='m'`).
+           Furthermore, the edge size will be converted to pixels (using :code:`edge_size/resolution`)
+           and rounded (see :code:`edge_size`). Therefore, the edge size when :code:`units!='px'` is just an approximation
+           if its value in meters is not divisible by the requested resolution.
+
     resolution : float | int, default = 10
         Pixel size in meters.
     stac : str, default = 'https://planetarycomputer.microsoft.com/api/stac/v1'
@@ -106,45 +120,49 @@ def create(
     ... )
     <xarray.DataArray (time: 27, band: 3, x: 128, y: 128)>
     """
+    # Harmonize units to pixels
+    if units != "px":
+        if units == "m":
+            edge_size = edge_size / resolution
+        else:
+            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
     )
 
     # Use Google Earth Engine
     if gee:
-        
+
         # Try to import ee, otherwise raise an ImportError
         try:
-            import xee
             import ee
+            import xee
         except ImportError:
             raise ImportError(
-                    '"earthengine-api" and "xee" could not be loaded. Please install them, or install "cubo" using "pip install cubo[ee]"'
-                )
+                '"earthengine-api" and "xee" could not be loaded. Please install them, or install "cubo" using "pip install cubo[ee]"'
+            )
 
         # Initialize Google Earth Engine with the high volume endpoint
-        ee.Initialize(opt_url='https://earthengine-highvolume.googleapis.com')
+        ee.Initialize(opt_url="https://earthengine-highvolume.googleapis.com")
 
         # Get BBox values in latlon
-        west = bbox_latlon['coordinates'][0][0][0]
-        south = bbox_latlon['coordinates'][0][0][1]
-        east = bbox_latlon['coordinates'][0][2][0]
-        north = bbox_latlon['coordinates'][0][2][1]
+        west = bbox_latlon["coordinates"][0][0][0]
+        south = bbox_latlon["coordinates"][0][0][1]
+        east = bbox_latlon["coordinates"][0][2][0]
+        north = bbox_latlon["coordinates"][0][2][1]
 
         # Create the BBox geometry in GEE
-        BBox = ee.Geometry.BBox(west,south,east,north)
+        BBox = ee.Geometry.BBox(west, south, east, north)
 
         # If the collection is string then access the Image Collection
-        if isinstance(collection,str):
+        if isinstance(collection, str):
             collection = ee.ImageCollection(collection)
 
         # Do the filtering: Bounds, time, and bands
         collection = (
-            collection
-                .filterBounds(BBox)
-                .filterDate(start_date,end_date)
-                .select(bands)
+            collection.filterBounds(BBox).filterDate(start_date, end_date).select(bands)
         )
 
         # Return the cube via xee
@@ -154,17 +172,21 @@ def create(
             geometry=BBox,
             scale=resolution,
             crs=f"EPSG:{epsg}",
-            chunks=dict()
+            chunks=dict(),
         )
 
         # Rename the coords to match stackstac names, also rearrange
-        cube = cube.rename(Y="y",X="x").to_array("band").transpose("time","band","y","x")
+        cube = (
+            cube.rename(Y="y", X="x")
+            .to_array("band")
+            .transpose("time", "band", "y", "x")
+        )
 
         # Delete all attributes
         cube.attrs = dict()
 
         # Get the name of the collection
-        collection = collection.get('system:id').getInfo()
+        collection = collection.get("system:id").getInfo()
 
         # Override the stac argument using the GEE STAC
         stac = "https://earthengine-stac.storage.googleapis.com/catalog/catalog.json"
@@ -216,13 +238,17 @@ def create(
             if attribute in cube.attrs:
                 del cube.attrs[attribute]
 
+    # Rounded edge size
+    rounded_edge_size = cube.x.shape[0]
+
     # New attributes
     cube.attrs = dict(
         collection=collection,
         stac=stac,
         epsg=epsg,
         resolution=resolution,
-        edge_size=edge_size,
+        edge_size=rounded_edge_size,
+        edge_size_m=rounded_edge_size * resolution,
         central_lat=lat,
         central_lon=lon,
         central_y=utm_coords[1],

docs/changelog.rst

@@ -1,6 +1,12 @@
 Changelog
 =========
 
+v2024.1.1
+---------
+
+- Added the :code:`units` argument to :code:`cubo.create()`.
+- Added support for :code:`scipy.constants` units.
+
 v2024.1.0
 ---------
 

docs/conf.py

@@ -27,7 +27,7 @@
 author = "David Montero Loaiza"
 
 # The full version, including alpha/beta/rc tags
-release = "2024.1.0"
+release = "2024.1.1"
 
 
 # -- General configuration ---------------------------------------------------

docs/index.rst

@@ -200,6 +200,42 @@ Main function: `create()`
       resolution=10,
    )
 
+Using different units for `edge_size`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+By default, the units of `edge_size` are pixels. But you can modify this using the `units` argument:
+
+.. code-block:: python
+
+   da = cubo.create(
+      lat=4.31,
+      lon=-76.2,
+      collection="sentinel-2-l2a",
+      bands=["B02","B03","B04"],
+      start_date="2021-06-01",
+      end_date="2021-06-10",
+      edge_size=1500,
+      units="m",
+      resolution=10,
+   )
+
+You can use "px" (pixels), "m" (meters), or any unit available in
+ `scipy.constants <https://docs.scipy.org/doc/scipy/reference/constants.html#units>`_.
+
+.. code-block:: python
+
+   da = cubo.create(
+      lat=4.31,
+      lon=-76.2,
+      collection="sentinel-2-l2a",
+      bands=["B02","B03","B04"],
+      start_date="2021-06-01",
+      end_date="2021-06-10",
+      edge_size=1.5,
+      units="kilo",
+      resolution=10,
+   )
+
 Using another endpoint
 ~~~~~~~~~~~~~~~~~~~~~~
 

docs/tutorials.rst

@@ -9,4 +9,5 @@ Tutorials
     tutorials/using_collections.ipynb
     tutorials/visualization_lexcube.ipynb
     tutorials/using_gee.ipynb
-    tutorials/stackstac.ipynb
+    tutorials/stackstac.ipynb
+    tutorials/edge_size_units.ipynb