Save stackstac output to file

[edornd]

Hi! Sorry in advance for the rather stupid question, but I couldn't find a "clean" solution so far.
In short, what is the recommended/canonical way to store the results of a stackstac.stack call to file?

For instance, the cloudless mosaic tutorial produces a RGB xarray, however simply calling to_netcdf or any other saving mechanism throws errors due to some "typeless" numpy arrays (proj:transform, proj:bbox, etc.):

ValueError:  unable to infer dtype on variable 'proj:transform'; xarray cannot serialize arbitrary Python objects

Where proj:transform is indeed a bit weird: array({0.0, 300000.0, 10.0, 6100020.0, -10.0}, dtype=object)

I managed to store the result as a NetCDF by dropping these variables, however I'm not a fan of dropping the georeferencing parts.
I'd like to be able to store it as COG, or even HDF5/NetCDF, but none of it work at the moment without some tinkering, I bet I'm missing something very obvious.
Given I'm new to Planetary and these tools, my questions are: is this supposed to happen? If so, how should I store results?

Thank you very much!

[TomAugspurger]

I think this is still an active area of development, but we should have better solutions in the next year or so.

Right now, the easiest output to use might be Zarr, but I don't know if that's compatible with your downstream use-case.

If you're writing to NetCDF or COG, you might want to look at how rioxarray manages coordinate information: https://corteva.github.io/rioxarray/html/getting_started/crs_management.html. It tends to write output that's compatible with GDAL.

https://github.com/TomAugspurger/xcog is a little proof of concept for writing out the chunked DataArray to a bunch of COGs (one per chunk) with STAC items per chunk. It's not really ready for use, but might be an avenue worth exploring.

[edornd]

To be fair, COG is not really mandatory for my use case: I'm building an image dataset for computer vision for internal use, so the format is not really enforced. It just seemed more "natural" to use an image format, that's all.

In any case, I'm having a similar issue saving with to_zarr (same stackstac-related issue with the proj attributes).
So far, the only stupid workaround I managed to come up with is to drop those vars:

data = data.drop_vars(["proj:shape", "proj:transform", "proj:bbox"])
# now I can save it in whatever format, netCDF and Zarr included
data.to_netcdf("data.nc")

However, as far as I know I basically lose the spatial information by doing so.
Another thing I have in mind is to convert/move these values (haven't figured it out yet).

It doesn't feel "clean" and a bit of a hack, but if there's no other proper way of doing that at the moment, I'm fine with it (as long as it works!).

[TomAugspurger]

Those properites like proj:shape​ are coming from the STAC metadata. They aren't used by any tools downstream of stackstac (like xarray, NetCDF, zarr, rioxarray) to interpret geospatial data, so they should be fine to drop. The rioxarray docs I linked to earlier has more information on which properties are used.

…

________________________________ From: Edoardo Arnaudo ***@***.***> Sent: Wednesday, November 1, 2023 10:32 AM To: microsoft/PlanetaryComputer ***@***.***> Cc: Tom Augspurger ***@***.***>; Comment ***@***.***> Subject: Re: [microsoft/PlanetaryComputer] Save stackstac output to file (Discussion #288) To be fair, COG is not really mandatory for my use case: I'm building an image dataset for computer vision for internal use, so the format is not really enforced. It just seemed more "natural" to use an image format, that's all. In any case, I'm having a similar issue saving with to_zarr (same stackstac-related issue with the proj attributes). So far, the only stupid workaround I managed to come up with is to drop those vars: data = data.drop_vars(["proj:shape", "proj:transform", "proj:bbox"]) # now I can save it in whatever format, netCDF and Zarr included data.to_netcdf("data.nc") However, as far as I know I basically lose the spatial information by doing so. Another thing I have in mind is to convert/move these values (haven't figured it out yet). It doesn't feel "clean" and a bit of a hack, but if there's no other proper way of doing that at the moment, I'm fine with it (as long as it works!). — Reply to this email directly, view it on GitHub<#288 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAKAOIROLXIH6FUM7Z4GWADYCJTSLAVCNFSM6AAAAAA6ZJAEISVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TINBXGAZTG>. You are receiving this because you commented.

[edornd]

Alright, thanks for the help, I'm starting to get the hang of it!
Looking around your examples and the docs, I managed to find a solution:

import rioxarray
import xarray as xr

# stack the available items, as per example
data = stackstac.stack(
    items,
    assets=["B04", "B03", "B02", "B08"],  # red, green, blue
    chunksize=512,
    bounds=rasterio.features.bounds(aoi.geometry),
    epsg=crs,
).where(...).assign_coords(...)

# build the mosaic
# keep_attrs=True is the important part, otherwise it doesn't get projected right (at least from my tests)
data.persist()
median = data.median(dim="time", keep_attrs=True).compute()

# then it's as easy as calling to_raster
data.rio.to_raster("test2.tif", driver="COG")

And that was it! I now have correct GeoTIFF files who are also (hopefully) COG-compliant. I skipped the rioxarray bits, my bad.
Thanks again for the useful inputs!

[rbavery]

this would be a great addition to the stackstac docs! didn't know stackstac was interoperable with rioxarray and the docs don't currently mention it: https://stackstac.readthedocs.io/en/latest/search.html?q=rioxarray&check_keywords=yes&area=default