Unable to access values in stacked data

[ATL-Flaneur]

I'm missing something in trying to build a single raster from Planetary Computer data. This is my first experience working with STAC data—I'm familiar with slicing and dicing local raster data with GDAL.

This is my trial code to plot DSM data derived from 3DEP lidar collects.

#!/usr/bin/env python3

import os.path
import planetary_computer
import stackstac
import pystac_client
import pystac
import rasterio
import rasterio.merge
import rasterio.plot

api_key = "" # Your API key goes here.
planetary_computer.settings.set_subscription_key(api_key)

URL = "https://planetarycomputer.microsoft.com/api/stac/v1/"
catalog = pystac_client.Client.open(URL)

collections = [
    "3dep-lidar-dsm"
]

FT_Myers = [-82.063, 26.541, -81.922, 26.698]
AOI = FT_Myers

time_range = "2000/2023"

search = catalog.search(
    collections=collections,
    bbox = AOI,
    datetime = time_range)

items = search.get_items()
signed_items = [planetary_computer.sign(item).to_dict() for item in items]

# Stack in WGS84 at 1 arc sec resolution.

data = stackstac.stack(
        signed_items,
        resolution=1/3600,
        epsg=4326,
        resampling=rasterio.enums.Resampling.average,
        bounds=FT_Myers
    ).drop("band").squeeze()

mosaic=stackstac.mosaic(data).persist()

mosaic.plot()

This results in a plot with expected axes but only a few dots:

.

The data type on mosaic and mosaic.data looks reasonable:

$ type(mosaic)
<class 'xarray.core.dataarray.DataArray'>
$ type(mosaic.data)
<class 'dask.array.core.Array'>
$ mosaic.data
dask.array<mosaic-aggregate, shape=(566, 508), dtype=float64, chunksize=(566, 508), chunktype=numpy.ndarray>

But doesn't seem like the data values are getting populated:

print(mosaic.data[0:9,0:9])
dask.array<getitem, shape=(9, 9), dtype=float64, chunksize=(9, 9), chunktype=numpy.ndarray>

I must be missing something on getting the data sent locally to my machine—some kind of lazy evaluation that isn't being triggered 🤔

[ATL-Flaneur]

The following code is giving me a plot that matches expectations for part of Ft. Myers, FL:

#!/usr/bin/env python3

import odc.stac
import os.path
import planetary_computer
import stackstac
import pystac_client
import pystac

from odc.geo.geobox import GeoBox

api_key = "" # Your API key goes here.
planetary_computer.settings.set_subscription_key(api_key)

URL = "https://planetarycomputer.microsoft.com/api/stac/v1/"
catalog = pystac_client.Client.open(URL)

collections = [
    "3dep-lidar-hag"
]

# AOI format is lower left and upper right corners.
# (x_min, y_min,x_max,y_max)

AOI = (-82.07, 26.54, -81.92, 26.70)

gb = GeoBox.from_bbox(AOI, "epsg:4326", resolution=1/3600)

time_range = "2010/2023"

search = catalog.search(
    collections=collections,
    bbox = AOI,
    datetime = time_range)

items = list(search.get_items())
signed_items = [planetary_computer.sign(item) for item in items]

gb = GeoBox.from_bbox(AOI, "epsg:4326", resolution=1/(3600*3))

data = odc.stac.load(
        signed_items,
        # assets=["data"],
        # chunksize=1024,
        geobox=gb
    ).squeeze()

# Enable "robust" so that color mapping doesn't get hijacked by missing values and outliers.
plt = data.data.plot.imshow(cmap="plasma", robust = True)

I'm seeing a lot of artifacts in the data, which can be seen along the bottom of the resulting plot. Maybe there's something wrong in how I'm accessing the data?

Denver (AOI = (-105.09, 39.73, -104.94, 39.83)) is pretty ugly. This is admittedly HAG data based on an old lidar collect instead of the newer 3DEP ones. Again, this seems bad enough that the problem is likely on my end.

[TomAugspurger]

Thanks for continuing to pursue this, and apologies for not responding earlier.

Looks like some of the parameters in the PDAL pipeline need to be adjusted, and we'll possibly need to reprocess some data.

[ATL-Flaneur]

I traded emails with @hobu and @mmcfarland on Friday and confirmed the presence of unfiltered noise in the AOIs I was downloading. Any estimate on a timeframe for reprocessing the raster data?

[hobu]

No timeframe, but we're working on it.