diff --git a/README.md b/README.md index 01361ea..74ea9f8 100644 --- a/README.md +++ b/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: diff --git a/cubo/__init__.py b/cubo/__init__.py index 4bb0d3a..d83a5a8 100644 --- a/cubo/__init__.py +++ b/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 " __all__ = [] diff --git a/cubo/cubo.py b/cubo/cubo.py index fd96319..dfc4668 100644 --- a/cubo/cubo.py +++ b/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,6 +120,13 @@ def create( ... ) """ + # 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 @@ -113,38 +134,35 @@ def create( # 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], diff --git a/docs/changelog.rst b/docs/changelog.rst index 4ebbe9d..a018b28 100644 --- a/docs/changelog.rst +++ b/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 --------- diff --git a/docs/conf.py b/docs/conf.py index b7a834b..d260d5e 100644 --- a/docs/conf.py +++ b/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 --------------------------------------------------- diff --git a/docs/index.rst b/docs/index.rst index 4a2e203..0402000 100644 --- a/docs/index.rst +++ b/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 `_. + +.. 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 ~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/tutorials.rst b/docs/tutorials.rst index 582804f..a9bd8f0 100644 --- a/docs/tutorials.rst +++ b/docs/tutorials.rst @@ -9,4 +9,5 @@ Tutorials tutorials/using_collections.ipynb tutorials/visualization_lexcube.ipynb tutorials/using_gee.ipynb - tutorials/stackstac.ipynb \ No newline at end of file + tutorials/stackstac.ipynb + tutorials/edge_size_units.ipynb \ No newline at end of file diff --git a/docs/tutorials/edge_size_units.ipynb b/docs/tutorials/edge_size_units.ipynb new file mode 100644 index 0000000..0bd2ee6 --- /dev/null +++ b/docs/tutorials/edge_size_units.ipynb @@ -0,0 +1,2579 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Change the units of `edge_size`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this tutorial we will learn how to change the units of `edge_size`." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import cubo" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By default, the units of `edge_size` are pixels:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/dmontero/anaconda3/envs/cubo/lib/python3.9/site-packages/stackstac/prepare.py:364: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.\n", + " times = pd.to_datetime(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'sentinel-2-l2a' (time: 2, band: 3, y: 64, x: 64)>\n",
+       "dask.array<fetch_raster_window, shape=(2, 3, 64, 64), dtype=float64, chunksize=(1, 1, 64, 64), chunktype=numpy.ndarray>\n",
+       "Coordinates: (12/47)\n",
+       "  * time                                     (time) datetime64[ns] 2021-06-13...\n",
+       "    id                                       (time) <U54 'S2B_MSIL2A_20210613...\n",
+       "  * band                                     (band) <U3 'B02' 'B03' 'B04'\n",
+       "  * x                                        (x) float64 6.011e+05 ... 6.023e+05\n",
+       "  * y                                        (y) float64 5.66e+06 ... 5.659e+06\n",
+       "    s2:water_percentage                      (time) float64 0.1445 0.1281\n",
+       "    ...                                       ...\n",
+       "    title                                    (band) <U20 'Band 2 - Blue - 10m...\n",
+       "    common_name                              (band) <U5 'blue' 'green' 'red'\n",
+       "    center_wavelength                        (band) float64 0.49 0.56 0.665\n",
+       "    full_width_half_max                      (band) float64 0.098 0.045 0.038\n",
+       "    epsg                                     int64 32632\n",
+       "    cubo:distance_from_center                (y, x) float64 908.2 ... 873.7\n",
+       "Attributes:\n",
+       "    collection:           sentinel-2-l2a\n",
+       "    stac:                 https://planetarycomputer.microsoft.com/api/stac/v1\n",
+       "    epsg:                 32632\n",
+       "    resolution:           20\n",
+       "    edge_size:            64\n",
+       "    edge_size_m:          1280\n",
+       "    central_lat:          51.079225\n",
+       "    central_lon:          10.452173\n",
+       "    central_y:            5659638.0946523\n",
+       "    central_x:            601722.4825156148\n",
+       "    time_coverage_start:  2021-06-01\n",
+       "    time_coverage_end:    2021-07-01
" + ], + "text/plain": [ + "\n", + "dask.array\n", + "Coordinates: (12/47)\n", + " * time (time) datetime64[ns] 2021-06-13...\n", + " id (time) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'sentinel-2-l2a' (time: 2, band: 3, y: 64, x: 64)>\n",
+       "dask.array<fetch_raster_window, shape=(2, 3, 64, 64), dtype=float64, chunksize=(1, 1, 64, 64), chunktype=numpy.ndarray>\n",
+       "Coordinates: (12/47)\n",
+       "  * time                                     (time) datetime64[ns] 2021-06-13...\n",
+       "    id                                       (time) <U54 'S2B_MSIL2A_20210613...\n",
+       "  * band                                     (band) <U3 'B02' 'B03' 'B04'\n",
+       "  * x                                        (x) float64 6.011e+05 ... 6.023e+05\n",
+       "  * y                                        (y) float64 5.66e+06 ... 5.659e+06\n",
+       "    s2:water_percentage                      (time) float64 0.1445 0.1281\n",
+       "    ...                                       ...\n",
+       "    title                                    (band) <U20 'Band 2 - Blue - 10m...\n",
+       "    common_name                              (band) <U5 'blue' 'green' 'red'\n",
+       "    center_wavelength                        (band) float64 0.49 0.56 0.665\n",
+       "    full_width_half_max                      (band) float64 0.098 0.045 0.038\n",
+       "    epsg                                     int64 32632\n",
+       "    cubo:distance_from_center                (y, x) float64 908.2 ... 873.7\n",
+       "Attributes:\n",
+       "    collection:           sentinel-2-l2a\n",
+       "    stac:                 https://planetarycomputer.microsoft.com/api/stac/v1\n",
+       "    epsg:                 32632\n",
+       "    resolution:           20\n",
+       "    edge_size:            64\n",
+       "    edge_size_m:          1280\n",
+       "    central_lat:          51.079225\n",
+       "    central_lon:          10.452173\n",
+       "    central_y:            5659638.0946523\n",
+       "    central_x:            601722.4825156148\n",
+       "    time_coverage_start:  2021-06-01\n",
+       "    time_coverage_end:    2021-07-01
" + ], + "text/plain": [ + "\n", + "dask.array\n", + "Coordinates: (12/47)\n", + " * time (time) datetime64[ns] 2021-06-13...\n", + " id (time) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'sentinel-2-l2a' (time: 2, band: 3, y: 76, x: 76)>\n",
+       "dask.array<fetch_raster_window, shape=(2, 3, 76, 76), dtype=float64, chunksize=(1, 1, 76, 76), chunktype=numpy.ndarray>\n",
+       "Coordinates: (12/47)\n",
+       "  * time                                     (time) datetime64[ns] 2021-06-13...\n",
+       "    id                                       (time) <U54 'S2B_MSIL2A_20210613...\n",
+       "  * band                                     (band) <U3 'B02' 'B03' 'B04'\n",
+       "  * x                                        (x) float64 6.01e+05 ... 6.025e+05\n",
+       "  * y                                        (y) float64 5.66e+06 ... 5.659e+06\n",
+       "    s2:water_percentage                      (time) float64 0.1445 0.1281\n",
+       "    ...                                       ...\n",
+       "    title                                    (band) <U20 'Band 2 - Blue - 10m...\n",
+       "    common_name                              (band) <U5 'blue' 'green' 'red'\n",
+       "    center_wavelength                        (band) float64 0.49 0.56 0.665\n",
+       "    full_width_half_max                      (band) float64 0.098 0.045 0.038\n",
+       "    epsg                                     int64 32632\n",
+       "    cubo:distance_from_center                (y, x) float64 1.078e+03 ... 1.0...\n",
+       "Attributes:\n",
+       "    collection:           sentinel-2-l2a\n",
+       "    stac:                 https://planetarycomputer.microsoft.com/api/stac/v1\n",
+       "    epsg:                 32632\n",
+       "    resolution:           20\n",
+       "    edge_size:            76\n",
+       "    edge_size_m:          1520\n",
+       "    central_lat:          51.079225\n",
+       "    central_lon:          10.452173\n",
+       "    central_y:            5659638.0946523\n",
+       "    central_x:            601722.4825156148\n",
+       "    time_coverage_start:  2021-06-01\n",
+       "    time_coverage_end:    2021-07-01
" + ], + "text/plain": [ + "\n", + "dask.array\n", + "Coordinates: (12/47)\n", + " * time (time) datetime64[ns] 2021-06-13...\n", + " id (time) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'sentinel-2-l2a' (time: 2, band: 3, y: 142, x: 142)>\n",
+       "dask.array<fetch_raster_window, shape=(2, 3, 142, 142), dtype=float64, chunksize=(1, 1, 142, 142), chunktype=numpy.ndarray>\n",
+       "Coordinates: (12/47)\n",
+       "  * time                                     (time) datetime64[ns] 2021-06-13...\n",
+       "    id                                       (time) <U54 'S2B_MSIL2A_20210613...\n",
+       "  * band                                     (band) <U3 'B02' 'B03' 'B04'\n",
+       "  * x                                        (x) float64 6.003e+05 ... 6.031e+05\n",
+       "  * y                                        (y) float64 5.661e+06 ... 5.658e+06\n",
+       "    s2:water_percentage                      (time) float64 0.1445 0.1281\n",
+       "    ...                                       ...\n",
+       "    title                                    (band) <U20 'Band 2 - Blue - 10m...\n",
+       "    common_name                              (band) <U5 'blue' 'green' 'red'\n",
+       "    center_wavelength                        (band) float64 0.49 0.56 0.665\n",
+       "    full_width_half_max                      (band) float64 0.098 0.045 0.038\n",
+       "    epsg                                     int64 32632\n",
+       "    cubo:distance_from_center                (y, x) float64 2.011e+03 ... 1.9...\n",
+       "Attributes:\n",
+       "    collection:           sentinel-2-l2a\n",
+       "    stac:                 https://planetarycomputer.microsoft.com/api/stac/v1\n",
+       "    epsg:                 32632\n",
+       "    resolution:           20\n",
+       "    edge_size:            142\n",
+       "    edge_size_m:          2840\n",
+       "    central_lat:          51.079225\n",
+       "    central_lon:          10.452173\n",
+       "    central_y:            5659638.0946523\n",
+       "    central_x:            601722.4825156148\n",
+       "    time_coverage_start:  2021-06-01\n",
+       "    time_coverage_end:    2021-07-01
" + ], + "text/plain": [ + "\n", + "dask.array\n", + "Coordinates: (12/47)\n", + " * time (time) datetime64[ns] 2021-06-13...\n", + " id (time) =2.0.3", "planetary_computer>=1.0.0", "pystac_client>=0.7.2", + "scipy>=1.12.0", "stackstac>=0.4.4", "xarray>=2023.6.0", ], diff --git a/tests/test_cubo.py b/tests/test_cubo.py index ab8ac12..219409c 100644 --- a/tests/test_cubo.py +++ b/tests/test_cubo.py @@ -35,7 +35,7 @@ def test_stackstac(self): end_date="2021-06-10", edge_size=32, resolution=10, - stackstac_kw=dict(xy_coords='center') + stackstac_kw=dict(xy_coords="center"), ) self.assertIsInstance(da, xr.DataArray)