Question on tile merging logic

[ljstrnadiii]

I am considering using this as an intermediate solution to write out a region or collection to an xarray zarr datastore.

Something like: given a collection of items, merge all tiles into one large xarray lazily, and write to zarr with specific chunk over dask.

Why? I have access patterns like random slices of much smaller regions than the native zoom of the underlying raster data.

Questions:

1. is this reasonable?
2. can someone point to the merge logic? e.g. how are lat/lon coordinates handled where two adjacent tiles meet?

Sweet project! Excited to give it a try soon.

[gjoseph92]

Are you talking about making a mosaic? Sounds like you're basically talking about flattening the stack of items, whether through a compositing method like mean(), median(), etc. or mosaic (picking the first valid pixel). The example from the docs on making composites may be helpful.

given a collection of items, merge all tiles into one large xarray lazily, and write to zarr with specific chunk over dask

This basically describes what stackstac.stack does, so a more thorough description here would help.

Note that doing this over a large spatial area is challenging, because it produces more tasks than dask can handle, see this thread: microsoft/PlanetaryComputer#12. There isn't a recommended way to do this yet.

[ljstrnadiii]

Yes, I think I mean mosaic. I mean to say given 3 rgb geotiffs that are contiguous in space (neighboring in this sense) of shape (100,100,3) stackstac will merge into one xarray dataset of shape (300,100,3) depending on what direction they are neighboring geotiffs from what I understand. Is that correct?

If so, how does the final dataset height, width and coordinates get estimated and how do the geotiffs get aligned?

[gjoseph92]

Well fist off, stackstac can't handle multi-band assets right now #62. Each band would have to be a separate GeoTIFF. So let's say shape (100, 100, 1) right now (actually (1, 100, 100), because we put the bands as the outer axis).

If you have three GeoTIFFs laid out like this, each with shape (1, 100, 100):

+--++--++--+
|  ||  ||  |
+--++--++--+

Then yes, stackstac.stack will produce an xarray DataArray like this with shape (3, 1, 100, 300) (this is items, bands, y, x)

+--++--++--+
|          |
+--++--++--+

how does the final dataset height, width and coordinates get estimated and how do the geotiffs get aligned?

This is all calculated by stackstac from the STAC metadata. We look at the proj:bbox, proj:shape, proj:transform, geometry, and proj:epsg fields to figure out the resolution and bounding-box of each item. Then we calculate a bounding-box that covers all the items, and pick the finest resolution of all the items. Basically, we try to pick a good "common grid" for all the items automatically. But you can override any of these choices by passing epsg=, resolution=, bounds/bounds_latlon= to stackstac.stack.

Then, when we read in each item, we reproject it to that common resolution and CRS, and clip it to those bounds. So all the geotiffs are inherently aligned, because they're being warped to that common grid as necessary.

[ljstrnadiii]

Got it. Thanks for sharing!

What is the priority on supporting multiple bands? That is pretty much a deal breaker for me I believe :/

From what I understand I would have to open each geotiff and save out a new one for each band, correct?

[gjoseph92]

That's one of the top priority items. However, working on stackstac isn't part of my job, so I don't have much time for it these days. So since it's not something I can add in an hour or two, it's less likely to happen that soon.

From what I understand I would have to open each geotiff and save out a new one for each band, correct?

Yes.