rioxarray interoperability¶

lazycogs generates DataArrays that are interoperable with rioxarray. This notebook demonstrates how you can use rioxarray's to_raster function to export a 3D lazycogs array (band, y, x) to a COG.

In [ ]:

Copied!





from pathlib import Path

import rioxarray  # noqa: F401 - registers the .rio accessor
import rustac
from pyproj import Transformer

import lazycogs

dst_crs = "epsg:5070"
dst_bbox = (250_000, 2_600_000, 350_000, 2_700_000)

# transform to epsg:4326 for STAC search
transformer = Transformer.from_crs(dst_crs, "epsg:4326", always_xy=True)
bbox_4326 = transformer.transform_bounds(*dst_bbox)

PARQUET = "data/midwest_summer_2025.parquet"

if not Path(PARQUET).exists():
    Path(PARQUET).parent.mkdir(exist_ok=True)
    await rustac.search_to(
        PARQUET,
        href="https://earth-search.aws.element84.com/v1",
        collections=["sentinel-2-c1-l2a"],
        datetime="2025-06-01/2025-09-30",
        bbox=bbox_4326,
        limit=100,
    )

print(f"Using {PARQUET}")

store = lazycogs.store_for(PARQUET, skip_signature=True)
from pathlib import Path

import rioxarray  # noqa: F401 - registers the .rio accessor
import rustac
from pyproj import Transformer

import lazycogs

dst_crs = "epsg:5070"
dst_bbox = (250_000, 2_600_000, 350_000, 2_700_000)

# transform to epsg:4326 for STAC search
transformer = Transformer.from_crs(dst_crs, "epsg:4326", always_xy=True)
bbox_4326 = transformer.transform_bounds(*dst_bbox)

PARQUET = "data/midwest_summer_2025.parquet"

if not Path(PARQUET).exists():
    Path(PARQUET).parent.mkdir(exist_ok=True)
    await rustac.search_to(
        PARQUET,
        href="https://earth-search.aws.element84.com/v1",
        collections=["sentinel-2-c1-l2a"],
        datetime="2025-06-01/2025-09-30",
        bbox=bbox_4326,
        limit=100,
    )

print(f"Using {PARQUET}")

store = lazycogs.store_for(PARQUET, skip_signature=True)

In [5]:

Copied!





da = lazycogs.open(
    PARQUET,
    crs=dst_crs,
    bbox=dst_bbox,
    resolution=100,
    time_period="P1D",
    bands=["red", "green", "blue"],
    dtype="int16",
    store=store,
)

da
da = lazycogs.open(
    PARQUET,
    crs=dst_crs,
    bbox=dst_bbox,
    resolution=100,
    time_period="P1D",
    bands=["red", "green", "blue"],
    dtype="int16",
    store=store,
)

da

Out[5]:

<xarray.DataArray (band: 3, time: 54, y: 1000, x: 1000)> Size: 324MB
[162000000 values with dtype=int16]
Coordinates:
  * band         (band) <U5 60B 'red' 'green' 'blue'
  * time         (time) datetime64[s] 432B 2025-06-02 2025-06-04 ... 2025-09-27
  * y            (y) float64 8kB 2.7e+06 2.7e+06 2.7e+06 ... 2.6e+06 2.6e+06
  * x            (x) float64 8kB 2.5e+05 2.502e+05 ... 3.498e+05 3.5e+05
    spatial_ref  int64 8B 0
Indexes:
  ┌ x        RasterIndex (crs=EPSG:5070)
  └ y
Attributes:
    grid_mapping:            spatial_ref
    zarr_conventions:        [{'schema_url': 'https://raw.githubusercontent.c...
    spatial:dimensions:      ['y', 'x']
    spatial:bbox:            (250000, 2600000, 350000, 2700000)
    spatial:transform_type:  affine
    spatial:transform:       [100.0, 0.0, 250000.0, 0.0, -100.0, 2700000.0]
    spatial:shape:           [1000, 1000]
    spatial:registration:    pixel
    proj:code:               EPSG:5070
    _FillValue:              0.0

rioxarray can export a 2- or 3D array so we need to eliminate the time dimension before exporting.

Now write the array to a cloud-optimized geotiff file!

In [8]:

Copied!

subset.rio.to_raster(Path("data/cog.tif"), driver="COG")
subset.rio.to_raster(Path("data/cog.tif"), driver="COG")