How to map GeoZarr to datasets with many multiple CRSs

[rbavery]

How should a STAC Collection containing many geospatial rasters that span many CRSs map to GeoZarr?

For example if I have Sentinel-2 or Landsat data spanning UTM zones, I don't want to have to reproject every scene to a common CRS in order to take advantage of chunking and sharding.

As I'm loading scenes from GeoZarr , making an ML prediction, and georeferencing the result, I also want to keep the original CRS so that I can keep the area preserving projection.

It seems like this has been discussed a bit here

It may be helpful to specify in which case a dataset has children dataset and how to index the hierarchy (Combined with STAC object, the GeoZarr Dataset might be restricted to a single variable and only contains multiscales as children alternative variables)

From #34

and @benbovy had a helpful proposal of how this might be achieved in Xarray. If this is to big of a change to Xarray maybe it would be possible to achieve this index of indexes in ran extension to Xarray? rioxarray?

pydata/xarray#9543 (comment)

curious if others have thoughts on this. I don't have a well formulated design yet, still in the stage of reading and digesting the discussions up to today.

[mdsumner]

My opinion is it shouldn't, if you want a standardized index set one up with the footprint of every source to a common crs, use that to index and filter to the sources that will contribute to a final "render" in that common crs, or to any other grid spec required.

This is how standard tools already work (the GDAL API for example), and how specific implementations that load from stac work: gdalcubes, opendatacube, stackstac, the GDAL STACIT driver and (more generally) the GDAL GTI driver.

Look at how odc-geo maps pixel index coordinates from various crs to the tile pixel index of a target source, I suppose you could store references that do that, VirtualiZarr-like, but it's so easily calculated in the fly for an on-demand grid spec I really don't see the value. (Ok maybe I've talked myself into illustrating something ... this was brought up by Deepak in terms of getting the index and weights from GDAL, so maybe that's what's really being asked here, to virtualize and undercut current libs). 😆

here's the mailing list pointers to how one might get "index and weights" from the warper: https://lists.osgeo.org/pipermail/gdal-dev/2024-April/058906.html

In opendatacube I think it's a slightly different motivation, because they are trying to cull-filter the mapping between scenes and which tile/s they contribute to , so you can line up sources and target tiles (aka 'chunks' ...) . So yeah, maybe I've framed this topic in a way I can negotiate now (!)

[benbovy]

From GDAL's GTI driver docs:

The tiles are listed as features of any GDAL supported vector format.

This would translate much more naturally to a GeoPandas' GeoDataFrame of Xarray DataArray objects than a plain data cube (see here for how this could work). The GeoDataFrame geometry column would store the footprint (bbox) of each scene in a common CRS so we can easily retrieve the scenes by spatial filtering. Basically the same approach than the one @mdsumner describes and that is indeed used by other tools (opendatacube, STAC tools, etc.), although those tools rely on other mechanisms than a GeoDataFrame for indexing the scenes.

This representation of individual scenes + geospatial index has some limitations, e.g., we cannot easily slice along other dimensions (time, band, etc.) as we could do with a plain datacube, but I guess that in most applications it is better anyway to do some spatial filtering first. It is also easier to reason about especially when the scenes are heterogeneous (i.e., different shapes and/or CRS), sparsely distributed and/or partially overlapping. In any case, this is much easier to implement than what I suggest in pydata/xarray#9543 (comment).

In terms of storage, with this approach we would still want to store the scenes each as a separate GeoZarr dataset. We could then store their footprint as a Geoparquet file or even into a basic relational database.

IIUC odc-geo's GridSpec is an alternative way of indexing the scenes in the specific (but pretty common?) case where those are arranged like a regular spatial grid. I'm wondering if in this case we couldn't have a Zarr "meta" dataset storing the information of each tile, i.e., their CRSs (e.g., using EPSG codes), their geo-transform (N+1 dimension array where the last dimension represents the 6 affine parameters) and even their footprint (assuming that GeoZarr would support vector geometries #28). That might be a little far-fetched, though, since we don't really need the full power of Zarr (cloud-optimized / n-arbitrary dimensions) to store this kind of information (pretty cheap, 2-dimensions only).

[benbovy]

Also xref geoarrow/geoarrow#24, which is relevant for this use case I think.

[RichardScottOZ]

What about xvec?

[benbovy]

Xvec and vector data cubes in general work best with geometries having attributes that also depend on other, non-geospatial dimensions like time, etc. Here this is a bit different: each geometry materializes the extent of a raster.

Similarly to a GeoDataFrame of DataArrays, maybe we could imagine an xvec DataArray of DataArrays? I'm not sure at all if/how this would work, though...
Would there be any advantage of using xarray DataArrays over simpler structures like a GeoDataFrame or odc-geo's GridSpec here?

Perhaps Xvec may be useful for loading and handling the Zarr "meta" dataset (cf. my suggestion above)? Off the top of my head something similar to how we use stackstac to query and turn a STAC collection into a DatArray, but here exclusively based on (Geo)Zarr and Xarray (Xvec).

[mdsumner]

From GDAL's GTI driver docs:

The tiles are listed as features of any GDAL supported vector format.

This would translate much more naturally to a GeoPandas' GeoDataFrame of Xarray DataArray objects than a plain data cube (see here for how this could work).

It *is a * vector layer when you treat it as one (it's an abstraction), when you treat it as a "GTI:" it is a raster (a 2D datacube with 1 or more vars), the entire point of the driver - and exactly like loading from stac into xarray via the various rasterio wrappers that exist now. I should have made that clearer, it doesn't present as a data frame when you access it via GTI, and the grid-resolving tasks are not invoked at all until you trigger a read-tasks. I'm not suggesting to crush tabular data into a Zarr, just that the idea of having multiple grid-specs somehow lined up in an array doesn't make much sense ... unless you are trying to index the details of how each footprint and underlying tiling maps into the final array (?). I think there is some utility in that.

[benbovy]

Thanks for the clarification @mdsumner ! So if I understand well GDAl GTI (collection of rasters -> virtual raster) is quite similar to stackstac (STAC items -> lazy xarray DataArray), right?

To be clear, I'm not suggesting either to crush tabular data into Zarr. I guess I was trying to say that a table of rasters (data frame) vs. a mosaic of rasters (single array) have both pros and cons, but that might be slightly off-topic.

I also agree that a multi-crs array doesn't make much sense and that something like a tile index would be useful. Two questions:

1. Is there any interest in adding the tile index itself to the GeoZarr specs? Or should we leave that to other specs (e.g., STAC + GeoZarr)?

2. If the answer to 1 is yes, how the spec would look like for the tile index? What do we need? Can / should we put both the tile index and all the individual tiles in the same Zarr dataset (leveraging Zarr's group hierarchy)?

[mdsumner]

So if I understand well GDAl GTI (collection of rasters -> virtual raster) is quite similar to stackstac (STAC items -> lazy xarray DataArray), right?

yes - I was a too reactive on the "crush" thing, I went and made an example and while it's amazingly simple to craft, it's quite abstract and needs some explaining. I think it's just not worth mapping disparate sources to an actual array, the array has a spec (rank, coordinates) and that's completely abstract, you can pour any sources into that so I guess it comes down to how to think about the catalogue of possible sources and when the mapping starts to occur.

(and I think that mapping is interesting, like for a given target grid specification, which sources are actually relevant, but maybe filtering spatially is best done after temporal and data type filtering etc)

your Qs I don't think so - I think the catalogue topic is very dynamic and will change and consolidate a lot in fairly short time. I'm concerned that STAC is treated specially, because why not have any old grid sources in a lighter collection - GTI (and the GeoParquet example, already compatible with GTI shows you can do that, STAC becomes a bespoke case folded in.

here's a very very simple GTI example I crafted, two disparate polar crs grids to a shared global one: https://gist.github.com/mdsumner/724f1db0bd0d211de6e3775486cd3df0

[maxrjones]

From previous conversations, I believe @ljstrnadiii may have some experience with the DataFrame of DataArray's approach to this problem. Len, am I remembering that right and do you have any insights to share for this conversation?

[ljstrnadiii]

Yeah, @maxrjones I used a similar schema as GTI for the purpose of indexing tiles for different collections (features), date times, and crs. Columns in my geopandas df are geometry (bbox in the geopandas df's crs coordinates), path to cog, collection name, datetime, crs.

My goal is to reproject and resample heterogeneous tiles (w.r.t. crs) into a single Xarray data array with the final dims: band, time, y, x. However, gdalwarp wasn't happy when I provided a list of tiles in different crs. This could be a user error thing lol. So, I first query relevant tiles from parquet/geopandas, then groupby crs, mosaic with vrt, warp to target crs, and finally mosaic warped vrts. I then do that per collection and per datetime.

I am actually really excited to try out GTI soon for larger mosaics because the current approach seems to take time and I think requires http request per tif to discover metadata to build the vrt. I can build GTI data while I'm ingesting the tifs and avoid the lookup and provide it per tif as recommended in the GTI docs.

My guess is that GTI will help a ton with that and seems like it supports heterogeneous tiles crs and resolution!

In regards to the post, I don't see how chunking is relevant across crs areas without reprojecting. I'm failing to follow here, but I don't have a good developer-level understanding of geozarr, GDAL, or Xarray.

I suppose I'd try GTI to reproject resample to one mosaic Xarray and then apply the model for predictions, etc. and reproject back to UTM per zone. That depends on information loss and is honestly a question I've had: does going from a particular UTM zone to 3857 (for example) and back cause data loss/error? It will certainly generate redundant data. Or maybe make one Xarray data array per UTM zone and chunk within utm zones mosaiked by GTI?

I have tried to use a dimension in the past for tile id for example, but that was awkward since each tile has different heights/widths and therefore coordinates indexed by tile id so concatenating over tile id was not ideal not to mention edge effects, etc.

[mdsumner]

However, gdalwarp wasn't happy when I provided a list of tiles in different

That definitely works, keen to follow up on this. (You can see in my gist above that it works). Certainly it has improved a lot in recent years, I used to encounter crashes a lot ~8 years ago but it's super robust now. I usually use it via the API directly, the high level Warp() function in C++ , or via python or R.

Super keen to explore GTI with you, GDAL is also very responsive and developing fast, and some parts don't get tested enough by the community (or tested too late via downstream packages).

🙏

[RichardScottOZ]

I have a problem like Len's currently - not tiles but global heterogeneous surveys at different resolutions and shapes and CRS - so 'domain' is a possible index - gravity, magnetics, etc.

Might not just be errors @ljstrnadiii but artefacts too?

[rbavery]

Wow thanks all for the quick and insightful comments! @mdsumner great recommendation on GTI that short example you provided is compelling. Looking forward to trying it out.

@ljstrnadiii you laid out a lot of options I've thought of trying and more I hadn't thought of, thanks!

In regards to the post, I don't see how chunking is relevant across crs areas without reprojecting.

My line of thinking was roughly what @benbovy discusses above with more clarity wrt to the Zarr "meta" dataset that tracks CRSs, geo-transform, and footprint. Each Zarr in a group would be for a single CRS. Reprojections would happen only on the fly if needed by the operation, like with GTI. Writing this out, it sounds like GTI with STAC geoparquet referencing COGs could be a great solution.

The only thing I'm wondering though is that Zarr v3 seems to have better support than other array formats for sharding, making it easy to rechunk, and compressing chunks with different methods. Or maybe (likely) this is all available in GDAL and I just need to chat with @mdsumner :)

I suppose I'd try GTI to reproject resample to one mosaic Xarray and then apply the model for predictions, etc. and reproject back to UTM per zone. That depends on information loss and is honestly a question I've had: does going from a particular UTM zone to 3857 (for example) and back cause data loss/error?

There's definitely some data corruption at higher latitudes but I'm very unsure on the extent of it and I think it is dependent on the size of the reprojected regions and latitude since mercator distorts more toward the poles. This makes me think it would be nice to tell GTI to reproject adjacent UTM scenes not to Web Mercator (or another virtual mosaic CRS) but something more locally accurate for the two scenes involved. Like one of the UTM projections of one of the scenes, determined by which UTM zone has the majority of pixel area.

Anyway I think the two questions posed by Benoit are good to discuss more! I'm curious to hear more from @mdsumner on his take that GeoZarr should not have a tile index.

the catalogue topic is very dynamic and will change and consolidate a lot in fairly short time.

I had the sense that GeoZarr is where this sort of OGC standardization on cataloguing should take place.

[RichardScottOZ]

However, gdalwarp wasn't happy when I provided a list of tiles in different

That definitely works, keen to follow up on this. (You can see in my gist above that it works). Certainly it has improved a lot in recent years, I used to encounter crashes a lot ~8 years ago but it's super robust now. I usually use it via the API directly, the high level Warp() function in C++ , or via python or R.

Super keen to explore GTI with you, GDAL is also very responsive and developing fast, and some parts don't get tested enough by the community (or tested too late via downstream packages).

🙏

How far have you seen this taken Michael? In terms of lots of CRS, size etc?