ValueError: conflicting sizes for dimension with single item / asset

StackSTAC

Turn a list of STAC items into a 4D xarray DataArray (dims: time, band, y, x), including reprojection to a common grid. The array is a lazy Dask array, so loading and processing the data in parallel—locally or on a cluster—is just a compute() call away.

For more information and examples, please see the documentation.

import stackstac
import pystac_client

URL = "https://earth-search.aws.element84.com/v1"
catalog = pystac_client.Client.open(URL)

stac_items = catalog.search(
    intersects=dict(type="Point", coordinates=[-105.78, 35.79]),
    collections=["sentinel-2-l2a"],
    datetime="2020-04-01/2020-05-01"
).get_all_items()

stack = stackstac.stack(stac_items)
print(stack)

<xarray.DataArray 'stackstac-fefccf3d6b2f9922dc658c114e79865b' (time: 13, band: 17, y: 10980, x: 10980)>
dask.array<fetch_raster_window, shape=(13, 17, 10980, 10980), dtype=float64, chunksize=(1, 1, 1024, 1024), chunktype=numpy.ndarray>
Coordinates: (12/24)
  * time                        (time) datetime64[ns] 2020-04-01T18:04:04 ......
    id                          (time) <U24 'S2B_13SDV_20200401_0_L2A' ... 'S...
  * band                        (band) <U8 'overview' 'visual' ... 'WVP' 'SCL'
  * x                           (x) float64 4e+05 4e+05 ... 5.097e+05 5.098e+05
  * y                           (y) float64 4e+06 4e+06 ... 3.89e+06 3.89e+06
    view:off_nadir              int64 0
    ...                          ...
    instruments                 <U3 'msi'
    created                     (time) <U24 '2020-09-05T06:23:47.836Z' ... '2...
    sentinel:sequence           <U1 '0'
    sentinel:grid_square        <U2 'DV'
    title                       (band) object None ... 'Scene Classification ...
    epsg                        int64 32613
Attributes:
    spec:        RasterSpec(epsg=32613, bounds=(399960.0, 3890220.0, 509760.0...
    crs:         epsg:32613
    transform:   | 10.00, 0.00, 399960.00|\n| 0.00,-10.00, 4000020.00|\n| 0.0...
    resolution:  10.0

Once in xarray form, many operations become easy. For example, we can compute a low-cloud weekly mean-NDVI timeseries:

lowcloud = stack[stack["eo:cloud_cover"] < 40]
nir, red = lowcloud.sel(band="B08"), lowcloud.sel(band="B04")
ndvi = (nir - red) / (nir + red)
weekly_ndvi = ndvi.resample(time="1w").mean(dim=("time", "x", "y")).rename("NDVI")
# Call `weekly_ndvi.compute()` to process ~25GiB of raster data in parallel. Might want a dask cluster for that!

Installation

pip install stackstac

Windows notes:

It's a good idea to use conda to handle installing rasterio on Windows. There's considerably more pain involved with GDAL-type installations using pip. Then pip install stackstac.

Things `stackstac` does for you:

Figure out the geospatial parameters from the STAC metadata (if possible): a coordinate reference system, resolution, and bounding box.
Transfer the STAC metadata into xarray coordinates for easy indexing, filtering, and provenance of metadata.
Efficiently generate a Dask graph for loading the data in parallel.
Mediate between Dask's parallelism and GDAL's aversion to it, allowing for fast, multi-threaded reads when possible, and at least preventing segfaults when not.
Mask nodata and rescale by STAC item asset scales/offsets.
Display data in interactive maps in a notebook, computed on the fly by Dask.

Limitations:

Raster data only! We are currently ignoring other types of assets you might find in a STAC (XML/JSON metadata, point clouds, video, etc.).
Single-band raster data only! Each band has to be a separate STAC asset—a separate red, green, and blue asset on each Item is great, but a single RGB asset containing a 3-band GeoTIFF is not supported yet.
COGs work best. "Normal" GeoTIFFs that aren't internally tiled, or don't have overviews, will see much worse performance. Sidecar files (like .msk files) are ignored for performance. JPEG2000 will probably work, but probably be slow unless you buy kakadu. Formats make a big difference.
BYOBlocksize. STAC doesn't offer any metadata about the internal tiling scheme of the data. Knowing it can make IO more efficient, but actually reading the data to figure it out is slow. So it's on you to set this parameter. (But if you don't, things should be fine for any reasonable COG.)
Doesn't make geospatial data any easier to work with in xarray. Common operations (picking bands, clipping to bounds, etc.) are tedious to type out. Real geospatial operations (shapestats on a GeoDataFrame, reprojection, etc.) aren't supported at all. rioxarray might help with some of these, but it has limited support for Dask, so be careful you don't kick off a huge computation accidentally.
I haven't even written tests yet! Don't use this in production. Or do, I guess. Up to you.

Roadmap:

Short-term:

Write tests and add CI (including typechecking)
Support multi-band assets
Easier access to s3://-style URIs (right now, you'll need to pass in gdal_env=stackstac.DEFAULT_GDAL_ENV.updated(always=dict(session=rio.session.AWSSession(...))))
Utility to guess blocksize (open a few assets)
Support item assets to provide more useful metadata with collections that use it (like S2 on AWS)
Rewrite dask graph generation once the Blockwise IO API settles

Long term (if anyone uses this thing):

Support other readers (aiocogeo?) that may perform better than GDAL for specific formats
Interactive mapping with xarray_leaflet, made performant with some Dask graph-rewriting tricks to do the initial IO at coarser resolution for lower zoom levels (otherwize zooming out could process terabytes of data)
Improve ergonomics of xarray for raster data (in collaboration with rioxarray)
Implement core geospatial routines (warp, vectorize, vector stats, GeoPandas/spatialpandas interop) in Dask

	{% set docname = env.doc2path(env.docname, base=False) %}
	.. note::
	You can view & download the original notebook
	`on Github <https://github.com/gjoseph92/stackstac/blob/main/docs/{{
	"../" + docname if docname.startswith("examples") else docname
	}}>`_.

	Or, `click here <https://cloud.coiled.io/gjoseph92/jobs/stackstac>`_
	to run these notebooks on Coiled with access to Dask clusters.
	"""

	# TODO enable binder once Coiled supports websocket clusters over 443.
	# (Binder blocks outbound traffic on all ports besides 22, 80, and 443, so we can't connect to Coiled on 8786.)
	# nbsphinx_prolog = """
	# {% set docname = env.doc2path(env.docname, base=False) %}
	# .. note::
	# You can run this notebook interactively here: \|Binder\|, or view & download the original
	# `on Github <https://github.com/gjoseph92/stackstac/blob/main/docs/{{ docname }}>`_.

	# .. \|Binder\| image:: https://mybinder.org/badge_logo.svg
	# :target: https://mybinder.org/v2/gh/gjoseph92/stackstac/main?urlpath=lab/tree/docs/{{ docname }}
	# """

	VSI_CACHE=True
	# ^ cache HTTP requests for opening datasets. This is critical for `ThreadLocalRioDataset`,
	# which re-opens the same URL many times---having the request cached makes subsequent `open`s
	# in different threads snappy.
	),
	read=dict(
	VSI_CACHE=False
	# ^ don't cache HTTP requests for actual data. We don't expect to re-request data,
	# so this would just blow out the HTTP cache that we rely on to make repeated `open`s fast
	# (see above)

	# There's no bbox, nor shape and transform. The only info we have is `item.bbox` in lat-lon.
	else:
	if item_bbox_proj is None:
	try:
	bbox_lonlat = item["bbox"]
	except KeyError:
	asset_bbox_proj = None
	else:
	# TODO handle error
	asset_bbox_proj = reproject_bounds(
	bbox_lonlat, 4326, out_epsg
	)
	item_bbox_proj = asset_bbox_proj
	# ^ so we can reuse for other assets

	"time": pd.to_datetime(
	[item["properties"]["datetime"] for item in items],
	infer_datetime_format=True,
	errors="coerce",
	),

gjoseph92 / stackstac Goto Github PK

stackstac's Introduction

StackSTAC

Installation

Things stackstac does for you:

Limitations:

Roadmap:

stackstac's People

Contributors

Stargazers

Watchers

Forkers

stackstac's Issues

Recommend Projects

Recommend Topics

Recommend Org

Things `stackstac` does for you: