stackstac.mosaic is not producing the desired output.

[ddjustina]

I'm working with the "io-lulc-9-class" collection. My "bbox_of_interest" is returning data from two tiles, as shown bellow.

stack = (
    stackstac.stack(
        items,
        dtype=np.ubyte,
        fill_value=255,
        epsg = 3857,
        bounds_latlon=bbox_of_interest,
        sortby_date=False,
    )
    .assign_coords(
        time=pd.to_datetime([item.properties["start_datetime"] for item in items])
        .tz_convert(None)
        .to_numpy()
    )
    .sortby("time")
)

merged = stack.squeeze().compute()

When I try to create a plot, an error occurs. I believe this error is due to having two data arrays with the same date, each corresponding to different tiles.

g = merged2.plot(col="time")
for ax in g.axes.flat:
    ax.set_axis_off()

ValueError: Coordinates used for faceting cannot contain repeated (nonunique) values.

This is my Dataset time coordinate:
time (time) datetime64[ns]
array(['2017-01-01T00:00:00.000000000', '2017-01-01T00:00:00.000000000',
'2018-01-01T00:00:00.000000000', '2018-01-01T00:00:00.000000000',
'2019-01-01T00:00:00.000000000', '2019-01-01T00:00:00.000000000',
'2020-01-01T00:00:00.000000000', '2020-01-01T00:00:00.000000000',
'2021-01-01T00:00:00.000000000', '2021-01-01T00:00:00.000000000',
'2022-01-01T00:00:00.000000000', '2022-01-01T00:00:00.000000000'],
dtype='datetime64[ns]')

I tried running the mosaic instead, but it is returning a Data Array with no time dimension, while I'm expecting a Data Array with multiple single years from 2017 to 2022.

merged = stackstac.mosaic(stack, dim="time", axis=None, nodata=0).squeeze().compute()

[TomAugspurger]

Can you share a copy-pastable example (including the code to get items)?

I think you want something like .groupby("time.year").agg(stackstac.mosaic, nodata=0), but would need to have an example to double-check.

[ddjustina]

Sorry, here is the code that comes before:

import dask.distributed
from matplotlib.colors import ListedColormap
import pystac_client

import cartopy.crs as ccrs
import numpy as np
import pandas as pd
import geopandas as gpd
import planetary_computer
import rasterio
import rasterio.features
import stackstac

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

bbox_of_interest = [-91.29318352234793, 35.240613512556436, -85.75802862493443, 37.50368465550935]
search = catalog.search(collections=["io-lulc-9-class"], bbox=bbox_of_interest)

items = search.item_collection()
print(f"Returned {len(items)} Items")

# The STAC metadata contains some information we'll want to use when creating
# our merged dataset. Get the EPSG code of the first item and the nodata value.
item = items[0]

# Create a single DataArray from out multiple resutls with the corresponding
# rasters projected to a single CRS. Note that we set the dtype to ubyte, which
# matches our data, since stackstac will use float64 by default.
stack = (
    stackstac.stack(
        items,
        dtype=np.ubyte,
        fill_value=255,
        epsg = 3857,
        bounds_latlon=bbox_of_interest,
        sortby_date=False,
    )
    .assign_coords(
        time=pd.to_datetime([item.properties["start_datetime"] for item in items])
        .tz_convert(None)
        .to_numpy()
    )
    .sortby("time")
)

client = dask.distributed.Client(processes=False)
print(f"/proxy/{client.scheduler_info()['services']['dashboard']}/status")

merged = stack.squeeze().compute()

[ddjustina]

I tried your suggestion, but I get the following:

stack.groupby('time.year').agg(stackstac.mosaic, nodata=0)

AttributeError: 'DataArrayGroupBy' object has no attribute 'agg'

[ddjustina]

it worked with odc.stac :

import dask.distributed
import dask.utils
import numpy as np
import planetary_computer as pc
import xarray as xr
from IPython.display import display
from pystac_client import Client

from odc.stac import configure_rio, stac_load

cfg = {
    "io-lulc-9-class": {
         "*": {"data_type": np.ubyte, "nodata": 0},
    },
    "*": {"warnings": "ignore"},
}

client = dask.distributed.Client()
configure_rio(cloud_defaults=True, client=client)
display(client)

catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")

query = catalog.search(
    collections=["io-lulc-9-class"],
    bbox=bbox_of_interest
)

items = list(query.get_items())
print(f"Found: {len(items):d} datasets")


crs = "EPSG:3857"
resolution = 100

ds = stac_load(
    items,
    crs=crs,
    resolution=resolution,
    groupby="solar_day",
    stac_cfg=cfg,
    bbox=bbox_of_interest,
)

[TomAugspurger]

Ah, it was .apply, not .agg:

by_year = stack.groupby("time.year").apply(stackstac.mosaic, dim="time", nodata=0)
by_year.isel(y=slice(2000,3000), x=slice(2000, 3000)).compute().plot(col="year")

---

odc.stac's default behavior is to group by time and do this kind of mosaicing. This can be controlled with the groupby keyword.