Can't leverage custom projection like modis ?

[Berhinj]

Planetary computer provides a very neat stac catalog to acces geotifs from MODIS. I would have loved to leverage stackstac instead of the suggested odc.stac library to access such data but it requires a proj:epsg properties for the stac item that modis doesn't per nature. Though the stac contains a proj:wkt2 which contains the information about modis custom projection.

I guess I'm probably missing something. How should it be managed?

import planetary_computer
import pystac_client
import rich.table
import rioxarray as rxr
import xarray as xr

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

# Boise, Idaho
latitude = 43.6
longitude = -116.2
buffer = 1
bbox = [longitude - buffer, latitude - buffer, longitude + buffer, latitude + buffer]

import stackstac
search = catalog.search(
        collections=["modis-11A1-061"],
        bbox=bbox,
        datetime="2022-01-01/2022-01-02",
    )
_ = search.matched()

xs = search.item_collection()
stack = stackstac.stack(
        xs,
    )

Note, the wkt2 projection is supported by rasterio.

rxr.open_rasterio(xs[0].assets["LST_Day_1km"].href).rio.crs

[gjoseph92]

See discussion in #67 (comment). The short answer is that it's not supported right now.

[Berhinj]

@gjoseph92 if I spend the time to implement it and make a PR - would it be something that could potentially be of interest in main (if done right) ?

[gjoseph92]

@Berhinj absolutely. The main work is less the implementation, but figuring out how we want to represent other CRSs—both as kwarg-input to stack, and more importantly, as metadata in xarray.

[gjoseph92]

@Berhinj absolutely. The main work is less the implementation, but figuring out how we want to represent other CRSs—both as kwarg-input to stack, and more importantly, as metadata in xarray.

[Berhinj]

if I'm not mistaken (I'm probably missing the point), currently the stack method supports geographic projections via the 'epsg' parameter (only accepting integers).

Per defaults, rioxarray already supports a wider range of projection system (WKT and geotransform) so we don't really have to reinvent the wheel on how to represent crs as metadata in xarray, right?

To do that rioxarray is relying on rasterio.crs.CRS objects which we could easily create from a wide variety of input (including epsg, pyproj object, and wkt) with the rasterio.crs.CRS.from_user_input().

How realistic would it be to shift from working with an epsg as an integer as currently implemented to a rasterio.crs.CRS object which we'd create with rasterio.crs.CRS.from_user_input()? I'd be happy to investigate but having your opinion on it first would be valuable.

From stack input perspective, it'd means 'epsg' should probably renamed with the more generic 'crs' then, right?