search

stac-utils / xpystac

For extending xarray.open_dataset to accept pystac objects
MIT License
32 stars 2 forks source link

bands is read as a dimension without a coordinate #36

Closed rbavery closed 8 months ago

rbavery commented 9 months ago 
import pystac_client
import planetary_computer
from shapely.geometry import Point
import stackstac
import xarray as xr

area_of_interest = Point((-120.963122, 37.025915)) # wright solar farm lon lat
range_old = "2015-01-01/2018-01-01"
range_new = "2020-01-01/2021-01-01"
max_cloud_cover = 15

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)
search_new = catalog.search(
    collections=["sentinel-2-l2a"], 
    intersects=area_of_interest, 
    datetime=range_new, 
    query={"eo:cloud_cover": {"lt": max_cloud_cover}, 's2:mgrs_tile':{'eq':'10SFF'}}, 
    max_items=4 
)
items_new = search_new.item_collection()

using stackstac directly shows alls dimensions are coordinates

new_range_xarr = stackstac.stack(items_new, assets=["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B09", "B10", "B11", "B12"])
new_range_xarr

using xpystac, only band is left without a coordinate

new_range_xarr = xr.open_dataset(items_new, engine="stac", assets=["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B09", "B10", "B11", "B12"])
new_range_xarr
jsignell commented 9 months ago

Yes! I think this is another flavor of https://github.com/pydata/xarray/issues/6038

One big difference between stackstac and xpystac is xpystac returns at dataset rather than a dataarray and to do that conversion it uses DataArray.to_dataset. I think the behavior should be that coords that depend only on band should be promoted to attrs on the data_vars. I'm opening up a PR to that effect.