stackstac
stackstac copied to clipboard
Support `proj:wkt2` if `proj:epsg` is null
As an example, see this Landsat Albers scene from the Landsat Collection 2 STAC API:
https://ibhoyw8md9.execute-api.us-west-2.amazonaws.com/prod/collections/landsat-c2l2alb-sr/items/LC08_L2SR_082024_20210806_20210811_02_A1_SR
I hacked stackstac to accept any CRS string in the epsg keyword rather than requiring an EPSG code which is basically just handing it along here:
https://github.com/gjoseph92/stackstac/blob/c431e091690d6b8af7eaba544af13bf8ebf5ebe6/stackstac/prepare.py#L515
Would be good to follow the same logic here that the new GDAL STAC Item driver does: https://github.com/OSGeo/gdal/pull/4138
which allows for:

Rather than drop the epsg code may be better to add deprecation warning, and add a new crs keyword. When parsing the STAC it can then construct the right CRS based on what is available for the Item/Asset.
@gjoseph92 I'm happy to issue a PR for this, lmk if I should, don't want to duplicate work if you had already been working on this.
I was intentional to only use epsg at first, because I really didn't want to get into the mess of CRS representations right out of the gate. I was hoping the rest of the community might come to consensus around how to represent CRSs in the xarray data model (see https://github.com/gjoseph92/stackstac/issues/50), but that still hasn't happened. The rioxarray model is complete but cumbersome and hard to understand for users. A plain crs field that can hold any type (EPSG code, proj string, wkt, wkt2, projjson, etc.) is probably where we'll end up, though I'm a little sad about it. I think xarray flexible indexes will be the real answer.
I think there are a number of places where we'd have to change logic to switch epsg to crs. I think if you start in geom_utils.py, update signatures to take crs instead of epsg, and then refactor everything calling those functions, it might work.
I'm open to a PR doing this, but I'd also like to hear your thoughts on the design. CRSs are a mess that I'm hesitant to wade into without a good plan.
As I've said in other places, my extremely strong preference would be for some other library (geoxarray?) to define the data model and provide all the utility functions (like our current geom_utils.py) for working with it, then have stackstac, rioxarray, etc. just deal with the IO and generate xarrays that conform with that data model.
This is an important issue for us. Some reference systems that are frequently used are not defined by an EPSG number such as all those defined with an ESRI: code (https://spatialreference.org/ref/esri/). EPSG doesn't have, for example, an equal area projection for North or South America.
Yes, does ESRI:102008 work? Never had cause to use this one.
I think only EPSG codes are recognized by stackstac, so I wouldn't know how to enter an ESRI: code. Just entering 102008 doesn't work.