Trouble merging rasters derived from USGS 3DEP LPC collects

[ATL-Flaneur]

I'm trying to switch a workflow for creating digital surface models (DSMs) from downloading USGS 3DEP LPC tiles, converting them to DSM tiles, and stitching the tiles together with GDAL to downloading a 3DEP lidar DSM from Planetary Computer.

I've been pulling from example code to get the following:

import planetary_computer
import stackstac
import pystac_client
import pystac
import rasterio

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]

data = stackstac.stack(
        signed_items,
        epsg=4326,
        bounds=FT_Myers)

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

This fails with :

Exception has occurred: RuntimeError
Error opening 'https://usgslidareuwest.blob.core.windows.net/usgs-3dep-cogs/usgs-cogs/USGS_LPC_FL_Southwest_A_2018_LAS_2019/dsm/USGS_LPC_FL_Southwest_A_2018_LAS_2019-dsm-2m-2-4.tif?st=2022-06-30T20%3A34%3A03Z&se=2022-07-02T20%3A34%3A03Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2022-07-01T03%3A23%3A46Z&ske=2022-07-08T03%3A23%3A46Z&sks=b&skv=2021-06-08&sig=AqW6uaqYru2hK4KU10vHHTis60u8VvsL2nuWhT3bOgA%3D': AttributeError("module 'rasterio' has no attribute 'parse_path'")

The above exception was the direct cause of the following exception:

  File "/Users/ayank/Documents/spectrum/propagation/tools/geo_download/src/dsm_download.py", line 79, in <module>
    mosaic=stackstac.mosaic(data).persist()

What am I getting wrong? It's very possible that I'm thinking about the nature of the returned data incorrectly. The 3DEP lidar-derived products don't appear to be a single seamless raster but a set of per-collect rasters.

I'm surprised that the returned data doesn't have a CRS and I need to specify one in the stack() call.

[TomAugspurger]

Can you share what version of rasterio you're using. The error AttributeError("module 'rasterio' has no attribute 'parse_path'"), is likely related to rasterio/rasterio#2365, if you're using a new enough rasterio.

I'm surprised that the returned data doesn't have a CRS and I need to specify one in the stack() call.

We're using proj:projjson for this collection. @hobu do you happen to remember why we chose that? I vaguely remember terms like "compund CRS"s.

stackstac currently just supports CRSs defined by EPSG code. I'm not sure if odc-stac will work with projjson, but it might be worth a shot.

The 3DEP lidar-derived products don't appear to be a single seamless raster but a set of per-collect rasters.

Yes, I believe that's correct. The 3dep:usgs_id on the item properties should point back to the original USGS collection.

[ATL-Flaneur]

Can you share what version of rasterio you're using.

I'm on rasterio 1.3b3 to get Python 3.10 support. Best I can tell, using 1.2.10 would get complicated since it appears to top out at Python 3.9. Up until now, I've relied on the system package manager (apt of Ubuntu Linux, MacPorts on macOS) for Python packages up until now; I'm learning pip along with Planetary Computer. Do I want to be setting up some prior version of Python and rasterio to run in parallel using virtualenv?

We're using proj:projjson for this collection.

Makes sense—the CRS of USGS LPC collects varies between collects. I've run across some odd ones…

[ATL-Flaneur]

I set up a parallel Python 3.9 install on my machine and installed rasterio 1.2.10. The stackstac.stack() operation no longer throws an exception and returns a data object that I can make method calls against. I see a thread pool fire up when the interpreter hits stackstac.mosaic() , but the call quickly executes and the Python interpreter exits. So something is going wrong in there.

I've replaced items = search.get_items() with items = list(search.get_items()) to persist the iterator.

[TomAugspurger]

Seems like something in rasterio / GDAL is segfaulting. I think you're hitting rasterio/rasterio#2415.

Thanks to rasterio/rasterio#2423 (review), that deprecation didn't quite go cleanly. So you're in a bit of a bind between rasterio versions and stackstack right now.

$ cat segfault.py 
# file segfault.py
import pystac
import planetary_computer
import stackstac

item = planetary_computer.sign(
    pystac.read_file("https://planetarycomputer.microsoft.com/api/stac/v1/collections/3dep-lidar-dsm/items/USGS_LPC_FL_Southwest_A_2018_LAS_2019-dsm-2m-4-4")
)

data = stackstac.stack(
    item,
    epsg=4326,
)

data.compute(scheduler="single-threaded")

Running that

$ python -X faulthandler segfault.py 
Fatal Python error: Segmentation fault

Current thread 0x00007f4643adf740 (most recent call first):
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/rasterio/crs.py", line 156 in to_epsg
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/stackstac/rio_reader.py", line 348 in _open
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/stackstac/rio_reader.py", line 381 in dataset
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/stackstac/rio_reader.py", line 385 in read
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/stackstac/to_dask.py", line 185 in fetch_raster_window
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/core.py", line 119 in _execute_task
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/core.py", line 149 in get
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/optimization.py", line 969 in __call__
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/core.py", line 119 in _execute_task
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/local.py", line 219 in execute_task
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/local.py", line 233 in <listcomp>
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/local.py", line 233 in batch_execute_tasks
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/local.py", line 537 in submit
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/local.py", line 489 in fire_tasks
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/local.py", line 494 in get_async
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/local.py", line 552 in get_sync
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/dask/base.py", line 573 in compute
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/xarray/core/dataset.py", line 861 in load
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/xarray/core/dataarray.py", line 921 in load
  File "/srv/conda/envs/notebook/lib/python3.8/site-packages/xarray/core/dataarray.py", line 947 in compute
  File "segfault.py", line 15 in <module>
Segmentation fault (core dumped)

Seems like rasterio.crs is having trouble with that STAC item's CRS info. That line in stackstac seems to segafult when doing rasterio.DatasetReader.crs Indeed,

cat segfault2.py 

import planetary_computer
import stackstac
import pystac
import rasterio
item = planetary_computer.sign(
    pystac.read_file("https://planetarycomputer.microsoft.com/api/stac/v1/collections/3dep-lidar-dsm/items/USGS_LPC_FL_Southwest_A_2018_LAS_2019-dsm-2m-4-4")
)
ds = rasterio.open(item.assets["data"].href)
print(ds.read_crs())

will trigger the segfault.

[ATL-Flaneur]

Seems like rasterio.crs is having trouble with that STAC item's CRS info.

Would immediately warping all the rasters to a simple, common CRS help?

I'm seeing that rasterio 1.3 is now out but I'm having a devil of a time install it—running MacPorts and pip then trying to juggle Python 3.9 and 3.10 is proving a nightmare.

[ATL-Flaneur]

So you're in a bit of a bind between rasterio versions and stackstack right now.

It required a grade & pave of my MacPorts install, but I got a base Python 3.10 system installed with it then used pip to install rasterio 1.3, which supports Python 3.10. Your segfault2.py now passes a smoke test and prints a CRS:

COMPD_CS["NAD83 / UTM zone 17N + NAVD88 height",PROJCS["NAD83 / UTM zone 17N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26917"]],VERT_CS["NAVD88 height",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","5703"]]]

Will continue working on this in the morning.

[gjoseph92]

stackstac 0.4.2 has just been released on pypi and supports rasterio 1.3 (though older versions will work as well, just throw a deprecation warning).

[ATL-Flaneur]

The code runs as expected with Python 3.10 after updating stackstac and rasterio.