gjoseph92 / stackstac

Turn a STAC catalog into a dask-based xarray
https://stackstac.readthedocs.io
MIT License
245 stars 49 forks source link

fetch_raster_window should use `fill_value` #113

Closed TomAugspurger closed 2 years ago

TomAugspurger commented 2 years ago

At https://github.com/gjoseph92/stackstac/blob/4029a681bc5c5bb4ff69c6d0767c241a17bb9384/stackstac/to_dask.py#L178, fetch_raster_window will fill with NaN, which can result in data.dtype not matching data.compute().dtype. It should use the user-provided fill_value instead.

A non-minimal reproducer, from https://github.com/microsoft/PlanetaryComputer/discussions/17

import numpy as np
import xarray as xr
import rasterio.features
import stackstac
import pystac_client
import pystac
import planetary_computer
from rasterio.enums import Resampling
import rioxarray
import dask
from dask_gateway import GatewayCluster
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")

aoi = {
    "type": "Polygon",
    "coordinates": [
        [
            [-100, 30],
            [-100, 66],
            [-47, 66],
            [-47, 30],
            [-100, 30]
        ]
    ],
}

search = catalog.search(
    #collections=["cop-dem-glo-90"], intersects=aoi
    collections=["alos-dem"], intersects=aoi
)
items = list(search.get_items())
signed_items = [planetary_computer.sign(item).to_dict() for item in items]

a = stackstac.stack(
        signed_items,
        assets=["data"],  
        chunksize=512,
        resolution=1000,
        epsg=32198,
        resampling=Resampling.average,
        dtype="int16",
        fill_value=0,
        # bounds=[-2009488, -715776, 1401061, 2597757],
      ).drop('band').squeeze()
# mosaic = stackstac.mosaic(data)
# mosaic
assert a.dtype == a[0, :10, :10].compute().dtype

I'll make a PR.