Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

STAC metadata as multidimensional coordinates #230

Draft
wants to merge 31 commits into
base: main
Choose a base branch
from

Conversation

gjoseph92
Copy link
Owner

@gjoseph92 gjoseph92 commented Nov 30, 2023

This refactors the STAC metadata -> xarray coordinates logic to support multi-dimensional coordinates. For example, imagine a field like rescaling (or href, for that matter), that's different for every asset of every item. Previously, we'd drop this field. Now, we can store it as a 2D array indexing the dimensions time, band.

This builds off of @Berhinj's work in #222, which provided a lot of inspiration for this design.

This is a significant rewrite of the metadata logic. The basic idea now is, for each field in STAC metadata:

  1. Make a 2D NumPy array (for the time, band dimensions)
  2. Iterate through all the STAC metadata, and write values into the array
  3. De-duplicate the array: if all values along a dimension are identical, drop that dimension. For example, the type field will be an MxN array of image/tiff; application=geotiff over and over; this is collapsed down into a 0D scalar.

Unfortunately, it does come at a slight performance cost: stackstac.stack on 10,000 Landast-8 items is 13.6% slower, going from 13.2s to 15s. (Of course, there's no performance difference with actually computing the results, which is where the bulk of time is spent anyway.)

There are a number of benefits, however:

  • Robust handling of JSON-y subfields. For example, the roles field, which contains variable-length lists of tags per band like ["data"], ["cloud", "cloud-shadow"], is now stored. Same with complex things like classification:bitfields. More interestingly, we could now easily retain geometry, bbox, etc. xref Store per-item bounding boxes as a coordinate #6
  • Generally, should be more robust to new fields/extensions adopted in STAC, or even inserted by users (they may not be formatted ideally, but almost certainly won't be dropped)
  • raster:bands subfields are now in coordinates, and supporting extensions like eo:bands, raster:bands, etc. is now generalized and easy to extend in the future
  • Finally have tests for coordinates logic!
  • rescale logic could be simplified into a one-liner after the stack is built (stack * stack.coords['raster:bands_scale'] + stack.coords['raster:bands_offset']), rather than plumbing scale/offset factors through many layers. (TBD if this makes a nasty Dask graph though.)

There's maybe even a future where the logic to generate the Dask array takes these coordinates as input, rather than raw STAC dicts or pystac items? Perhaps if we can someday store the geometries as a GeoSeries, or something like that.

Closes #216

types are so messed up, idk what to do with this, or if it matters
- unnesting is slow. can optimize.
    - could probably move into inner loop, so we don't have to iterate over values twice
- deobjectifying is slowest.
makes surprisingly little difference. even with micro-optimization, this is just slow python.

could maybe replace `isinstance` with `type is`; beyond that, idk.
this gets us to just 3x slower (was 4-5x before, with deobjectifying)

string fields are now left as objects, mostly because we'd have to scan them to see if any values were missing before `astype('U')`, because otherwise missing values would just be the string `'None'`
tests aren't quite there
@Berhinj
Copy link

Berhinj commented Dec 6, 2023

@gjoseph92 what is preventing you from merging? Would a deep code review help? I tried going through it quickly, but i's too far in my head and there has been a looot of changes from my initial code but in a good way though.

@gjoseph92
Copy link
Owner Author

@Berhinj simply haven't gotten around to taking another look at it. I wanted to give it a few days before taking another pass at self-review. If you could try out this branch on your use case and confirm it helps, that would be very helpful to get a real-world test!

@Berhinj
Copy link

Berhinj commented Dec 11, 2023

@gjoseph92 just tried it and compared to the results I was getting from my PR, the multicoordinnates aspect is working good!

Issue though, I notice that the "earthsearch:boa_offset_applied" boolean array from sentinel 2 became float...
Same for the "s2:processing_baseline" coords which went from characters to float

@gjoseph92
Copy link
Owner Author

I notice that the "earthsearch:boa_offset_applied" boolean array from sentinel 2 became float

Good point. I should special-case booleans.

@gjoseph92 just tried it and compared to the results I was getting from my PR, the multicoordinnates aspect is working good!

Same for the "s2:processing_baseline" coords which went from characters to float

I can't reproduce this. I'm getting array('05.09', dtype=object), which is still a string. I'm just running through the basic example from the docs for this.

@gjoseph92
Copy link
Owner Author

@Berhinj booleans should be handled correctly now; I'm getting earthsearch:boa_offset_applied as a bool now. Do you have any other feedback here?

@Berhinj
Copy link

Berhinj commented Jan 3, 2024

@gjoseph92 no, that's perfect, thanks

@Berhinj
Copy link

Berhinj commented Jan 19, 2024

any chance this will be merged soon? :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Multidimensional coordinnates is not supported
2 participants