Reading and Downloading DEM and LPC Products

The coincident.io.xarray and coincident.io.download modules support the option to load DEMs into memory via odc-stac or download DEMs into local directories.

There is also support for Lidar Point Cloud (LPC) spatial filtering for aerial lidar catalogs, where the user can return a GeoDataFrame with the respective .laz tile filename, download url, and geometry (epsg 4326) for each tile intersecting an input aoi. These laz files can then be downloaded locally with coincident.io.download.download_files()

There is specific support for USGS 3DEP EPT readers where the user can return a PDAL pipeline configured with the EPT URL, the AOI’s bounds, and polygon WKT, all in the EPT’s spatial reference system.

Note

Coincident does not support the processing of lidar point cloud products. Please see the lidar_tools repository for information on processing the returned GeoDataFrame with lidar point cloud products.

import coincident
import geopandas as gpd
from shapely.geometry import box

/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/src/coincident/io/download.py:25: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  from tqdm.autonotebook import tqdm

3DEP and NEON overlapping Flights

Note

For all of these functions, you will need identification metadata from the coincident.search.search functions for each respective catalog

Search

workunit = "CO_CentralEasternPlains_1_2020"
df_wesm = coincident.search.wesm.read_wesm_csv()
gf_usgs = coincident.search.wesm.load_by_fid(
    df_wesm[df_wesm.workunit == workunit].index
)

gf_usgs

	workunit	workunit_id	project	project_id	start_datetime	end_datetime	ql	spec	p_method	dem_gsd_meters	...	seamless_category	seamless_reason	lpc_link	sourcedem_link	metadata_link	geometry	collection	datetime	dayofyear	duration
0	CO_CentralEasternPlains_1_2020	192973	CO_CentralEasternPlains_2020_D20	192976	2020-05-09	2020-06-10	QL 2	USGS Lidar Base Specification 2.1	linear-mode lidar	1.0	...	Meets	Meets 3DEP seamless DEM requirements	https://rockyweb.usgs.gov/vdelivery/Datasets/S...	https://prd-tnm.s3.amazonaws.com/index.html?pr...	https://prd-tnm.s3.amazonaws.com/index.html?pr...	MULTIPOLYGON (((-103.7138 39.8835, -103.7138 3...	3DEP	2020-05-25	146	32

1 rows × 33 columns

# We will examine the 'sourcedem' 1m
# NOTE: it's important to note the *vertical* CRS of the data
gf_usgs.iloc[0]

workunit                                 CO_CentralEasternPlains_1_2020
workunit_id                                                      192973
project                                CO_CentralEasternPlains_2020_D20
project_id                                                       192976
start_datetime                                      2020-05-09 00:00:00
end_datetime                                        2020-06-10 00:00:00
ql                                                                 QL 2
spec                                  USGS Lidar Base Specification 2.1
p_method                                              linear-mode lidar
dem_gsd_meters                                                      1.0
horiz_crs                                                          6342
vert_crs                                                           5703
geoid                                                           GEOID18
lpc_pub_date                                        2020-12-11 00:00:00
lpc_update                                                          NaT
lpc_category                                                      Meets
lpc_reason                                  Meets 3DEP LPC requirements
sourcedem_pub_date                                  2020-12-10 00:00:00
sourcedem_update                                                    NaT
sourcedem_category                                                Meets
sourcedem_reason                     Meets 3DEP source DEM requirements
onemeter_category                                                 Meets
onemeter_reason                         Meets 3DEP 1-m DEM requirements
seamless_category                                                 Meets
seamless_reason                    Meets 3DEP seamless DEM requirements
lpc_link              https://rockyweb.usgs.gov/vdelivery/Datasets/S...
sourcedem_link        https://prd-tnm.s3.amazonaws.com/index.html?pr...
metadata_link         https://prd-tnm.s3.amazonaws.com/index.html?pr...
geometry              MULTIPOLYGON (((-103.7138 39.8835, -103.7138 3...
collection                                                         3DEP
datetime                                            2020-05-25 00:00:00
dayofyear                                                           146
duration                                                             32
Name: 0, dtype: object

Now, we’ll explore an overlapping NEON flight

gf_neon = coincident.search.search(
    dataset="neon", intersects=gf_usgs, datetime=["2020"]
)

/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/.pixi/envs/docs/lib/python3.14/site-packages/geopandas/array.py:411: UserWarning: Geometry is in a geographic CRS. Results from 'sjoin_nearest' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  warnings.warn(

gf_neon.head()

	id	title	start_datetime	end_datetime	product_url	geometry	collection
0	ARIK	Arikaree River NEON	2020-06-01	2020-06-30	https://data.neonscience.org/api/v0/data/DP3.3...	POINT (-102.44715 39.75821)	NEON

m = gf_usgs.explore(color="black")
gf_neon.explore(m=m, column="id", popups=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

Subset data

We will evaluate small subset of this overlap for deomstrative purposes.

Let’s subset based on some contextual LULC data

gf_wc = coincident.search.search(
    dataset="worldcover",
    intersects=gf_neon,
    datetime=["2020"],
)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[7], line 1
----> 1 gf_wc = coincident.search.search(
      2     dataset="worldcover",
      3     intersects=gf_neon,
      4     datetime=["2020"],

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/src/coincident/search/main.py:85, in search(dataset, intersects, datetime, **kwargs)
     77 # NOTE: not very robust, explode() demotes MultiPolygons to single Polygon (seems many GeoJSONs have this)
     78 # ANd 'exterior' not available for Multipolygons, just
     79 # NOTE: force_2d as some STAC searches fail with 3D polygons
     80 # https://github.com/uw-cryo/coincident/issues/101#issuecomment-3104277451
     81 # shapely_geometry = intersects.geometry.force_2d().explode().iloc[0]
     82 # Since coincident PCD fixtures are MultiPolygon
     83 shapely_geometry = intersects.geometry.force_2d().convex_hull.iloc[0]
---> 85 if not shapely_geometry.exterior.is_ccw:
     86     shapely_geometry = (
     87         shapely_geometry.reverse()
     88     )  # Apparently NASA CMR enforces polygon CCW order
     90 aoi = _pystac_client._format_intersects(shapely_geometry)  # to JSON geometry

AttributeError: 'Point' object has no attribute 'exterior'

dswc = coincident.io.xarray.to_dataset(
    gf_wc,
    bands=["map"],
    aoi=gf_neon,
    mask=True,
    resolution=0.00027,  # ~30m
)
dswc = dswc.rename(map="landcover")
dswc = dswc.compute()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 2
      1 dswc = coincident.io.xarray.to_dataset(
----> 2     gf_wc,
      3     bands=["map"],
      4     aoi=gf_neon,
      5     mask=True,

NameError: name 'gf_wc' is not defined

# arbitrary bbox that will be our subset area (all cropland)
bbox_geometry = box(-102.505, 39.675, -102.49, 39.685)
aoi = gpd.GeoDataFrame(geometry=[bbox_geometry], crs="EPSG:4326")

ax = coincident.plot.plot_esa_worldcover(dswc.landcover)
aoi.plot(ax=ax, facecolor="none", edgecolor="black", linestyle="--", linewidth=2)
from matplotlib.lines import Line2D

custom_line = Line2D([0], [0], color="black", linestyle="--", lw=2)
ax.legend([custom_line], ["Area of Interest"], loc="upper right", fontsize=10)
ax.set_title("ESA WorldCover");

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 ax = coincident.plot.plot_esa_worldcover(dswc.landcover)
      2 aoi.plot(ax=ax, facecolor="none", edgecolor="black", linestyle="--", linewidth=2)
      3 from matplotlib.lines import Line2D
      4 

NameError: name 'dswc' is not defined

Actually read in the DEMs

datetime_str = gf_neon.end_datetime.dt.strftime("%Y-%m-%d").item()
site_id = gf_neon.id.item()
datetime_str, site_id

('2020-06-30', 'ARIK')

%%time
da_neon_dem = coincident.io.xarray.load_neon_dem(
    aoi, datetime_str=datetime_str, site_id=site_id, product="dsm"
)

CPU times: user 45.3 ms, sys: 2.04 ms, total: 47.4 ms
Wall time: 158 ms

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[12], line 1
----> 1 get_ipython().run_cell_magic('time', '', 'da_neon_dem = coincident.io.xarray.load_neon_dem(\n    aoi, datetime_str=datetime_str, site_id=site_id, product="dsm"\n)\n')

File <timed exec>:1
----> 1 'Could not get source, probably due dynamically evaluated source code.'

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/src/coincident/io/xarray.py:437, in load_neon_dem(aoi, datetime_str, site_id, product, res, clip)
    434 aoi_geom_utm = aoi_utm.union_all()
    436 # 3: Query the NEON API
--> 437 data_json = query_neon_data_api(site_id, month_str)
    438 if (
    439     data_json.get("data", {}).get("release") == "PROVISIONAL"
    440     and len(data_json.get("data", {}).get("files", [])) == 0
    441 ):
    442     msg_provisional = (
    443         f"Data for {site_id} in {month_str} has PROVISIONAL status with no files."
    444     )

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/src/coincident/search/neon_api.py:297, in query_neon_data_api(site_id, month_str, product_code)
    293 if response.status_code != 200:
    294     msg_neon_fail = (
    295         f"NEON API request failed with status code {response.status_code}"
    296     )
--> 297     raise RuntimeError(msg_neon_fail)
    298 json_data: dict[str, Any] = response.json()  # Explicit cast for MyPy
    299 return json_data

RuntimeError: NEON API request failed with status code 403

da_neon_dem

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 1
----> 1 da_neon_dem

NameError: name 'da_neon_dem' is not defined

usgs_project = gf_usgs["project"].item()
usgs_project

'CO_CentralEasternPlains_2020_D20'

%%time
da_usgs_dem = coincident.io.xarray.load_usgs_dem(aoi, usgs_project)

CPU times: user 126 ms, sys: 12.6 ms, total: 139 ms
Wall time: 4.78 s

da_usgs_dem

<xarray.DataArray 'elevation' (y: 1147, x: 1319)> Size: 6MB
dask.array<getitem, shape=(1147, 1319), dtype=float32, chunksize=(1147, 1319), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 9kB 4.396e+06 4.396e+06 ... 4.395e+06 4.395e+06
  * x            (x) float64 11kB 7.14e+05 7.14e+05 ... 7.153e+05 7.153e+05
    band         int64 8B 1
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:             Area
    STATISTICS_MAXIMUM:        1290.7430419922
    STATISTICS_MEAN:           1256.105713952
    STATISTICS_MINIMUM:        1201.7416992188
    STATISTICS_STDDEV:         11.399835977004
    STATISTICS_VALID_PERCENT:  100

xarray.DataArray

'elevation'

• y: 1147
• x: 1319

• dask.array<chunksize=(1147, 1319), meta=np.ndarray>

  Array Chunk Bytes 5.77 MiB 5.77 MiB Shape (1147, 1319) (1147, 1319) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • y
    (y)
    float64
    4.396e+06 4.396e+06 ... 4.395e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([4395807.5, 4395806.5, 4395805.5, ..., 4394663.5, 4394662.5, 4394661.5],
          shape=(1147,))

  • x
    (x)
    float64
    7.14e+05 7.14e+05 ... 7.153e+05

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([713960.5, 713961.5, 713962.5, ..., 715276.5, 715277.5, 715278.5],
          shape=(1319,))

  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["NAD83 / UTM zone 13N",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",-105],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","26913"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    NAD83

    horizontal_datum_name :

    North American Datum 1983

    projected_crs_name :

    NAD83 / UTM zone 13N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -105.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["NAD83 / UTM zone 13N",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",-105],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","26913"]]

    GeoTransform :

    713960.0 1.0 0.0 4395808.0 0.0 -1.0

    array(0)

• Attributes: (6)

  AREA_OR_POINT :

  Area

  STATISTICS_MAXIMUM :

  1290.7430419922

  STATISTICS_MEAN :

  1256.105713952

  STATISTICS_MINIMUM :

  1201.7416992188

  STATISTICS_STDDEV :

  11.399835977004

  STATISTICS_VALID_PERCENT :

  100

da_usgs_dem.coarsen(x=5, y=5, boundary="trim").mean().plot.imshow();

Download

Note

coincident.io.download.download_neon_dem needs the NEON site’s start_datetime OR end_datetime to work

gf_neon

	id	title	start_datetime	end_datetime	product_url	geometry	collection
0	ARIK	Arikaree River NEON	2020-06-01	2020-06-30	https://data.neonscience.org/api/v0/data/DP3.3...	POINT (-102.44715 39.75821)	NEON

local_output_dir = "/tmp"

coincident.io.download.download_neon_dem(
    aoi=aoi,
    datetime_str=gf_neon.end_datetime.dt.strftime("%Y-%m-%d").item(),
    site_id=gf_neon.id.item(),
    product="dsm",
    output_dir=local_output_dir,
)

# USGS_1M_13_x71y440_CO_CentralEasternPlains_2020_D20.tif:  236MB
coincident.io.download.download_usgs_dem(
    aoi=aoi,
    project=usgs_project,
    output_dir=local_output_dir,
    save_parquet=True,  # save a STAC-like geoparquet of the tiles you download
)

Finally, you can grab the LPC tile metadata. For USGS 3DEP data, you can also return a PDAL pipeline based on the available EPT data. This PDAL pipeline will be returned as a JSON file where the user can add their own custom parameters (additional filters, writers, etc.) to this pipeline dictionary before executing it with PDAL.

Note

coincident.io.download.fetch_lpc_tiles needs the NEON site’s end_datetime OR start_datetime to work if you are in fact accessing NEON data

%%time
gf_neon_lpc_tiles = coincident.io.download.fetch_lpc_tiles(
    aoi=aoi,
    dataset_id=gf_neon.id.item(),
    provider="NEON",
    datetime_str=gf_neon.end_datetime.dt.strftime("%Y-%m-%d").item(),
)

CPU times: user 46.8 ms, sys: 924 μs, total: 47.7 ms
Wall time: 152 ms

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[20], line 1
----> 1 get_ipython().run_cell_magic('time', '', 'gf_neon_lpc_tiles = coincident.io.download.fetch_lpc_tiles(\n    aoi=aoi,\n    dataset_id=gf_neon.id.item(),\n    provider="NEON",\n    datetime_str=gf_neon.end_datetime.dt.strftime("%Y-%m-%d").item(),\n)\n')

File <timed exec>:1
----> 1 'Could not get source, probably due dynamically evaluated source code.'

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/src/coincident/io/download.py:1192, in fetch_lpc_tiles(aoi, dataset_id, provider, datetime_str, output_dir)
   1190         msg_no_neon_date = "datetime_str is required for NEON provider"
   1191         raise ValueError(msg_no_neon_date)
-> 1192     return _fetch_neon_lpc_tiles(
   1193         aoi, datetime_str=datetime_str, site_id=dataset_id, output_dir=output_dir
   1194     )
   1195 if prov == "ncalm":
   1196     return _fetch_ncalm_lpc_tiles(
   1197         aoi, dataset_name=dataset_id, output_dir=output_dir
   1198     )

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/src/coincident/io/download.py:648, in _fetch_neon_lpc_tiles(aoi, datetime_str, site_id, output_dir)
    645 aoi_u = aoi_utm.union_all()
    647 # 3. Query NEON API for LPC product
--> 648 data = query_neon_data_api(site_id, month_str, product_code="DP1.30003.001")
    649 files = data["data"]["files"]
    651 # 4. Filter by filename and AOI intersection

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/latest/src/coincident/search/neon_api.py:297, in query_neon_data_api(site_id, month_str, product_code)
    293 if response.status_code != 200:
    294     msg_neon_fail = (
    295         f"NEON API request failed with status code {response.status_code}"
    296     )
--> 297     raise RuntimeError(msg_neon_fail)
    298 json_data: dict[str, Any] = response.json()  # Explicit cast for MyPy
    299 return json_data

RuntimeError: NEON API request failed with status code 403

gf_neon_lpc_tiles.head()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[21], line 1
----> 1 gf_neon_lpc_tiles.head()

NameError: name 'gf_neon_lpc_tiles' is not defined

%%time
gf_usgs_lpc_tiles = coincident.io.download.fetch_lpc_tiles(
    aoi=aoi, dataset_id=usgs_project, provider="USGS"
)

CPU times: user 23.9 ms, sys: 1.96 ms, total: 25.9 ms
Wall time: 4.08 s

gf_usgs_lpc_tiles.head()

	name	url	geometry
0	5fded821d34e30b9123e230c	https://rockyweb.usgs.gov/vdelivery/Datasets/S...	POLYGON ((-102.50479 39.66904, -102.50479 39.6...
1	5fded821d34e30b9123e230e	https://rockyweb.usgs.gov/vdelivery/Datasets/S...	POLYGON ((-102.50447 39.67804, -102.50447 39.6...
2	5fded82fd34e30b9123e235a	https://rockyweb.usgs.gov/vdelivery/Datasets/S...	POLYGON ((-102.49314 39.66879, -102.49314 39.6...
3	5fded830d34e30b9123e235c	https://rockyweb.usgs.gov/vdelivery/Datasets/S...	POLYGON ((-102.49281 39.67779, -102.49281 39.6...
4	5fded837d34e30b9123e23a8	https://rockyweb.usgs.gov/vdelivery/Datasets/S...	POLYGON ((-102.48149 39.66854, -102.48149 39.6...

m = gf_usgs_lpc_tiles.explore(color="black")
gf_neon_lpc_tiles.explore(m=m)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[24], line 2
      1 m = gf_usgs_lpc_tiles.explore(color="black")
----> 2 gf_neon_lpc_tiles.explore(m=m)

NameError: name 'gf_neon_lpc_tiles' is not defined

Now, we can download the laz files

coincident.io.download.download_files(
    gf_neon_lpc_tiles["url"], output_dir=local_output_dir
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[25], line 2
      1 coincident.io.download.download_files(
----> 2     gf_neon_lpc_tiles["url"], output_dir=local_output_dir
      3 )

NameError: name 'gf_neon_lpc_tiles' is not defined

pdal_pipeline = coincident.io.download.build_usgs_ept_pipeline(
    aoi, workunit=gf_usgs.workunit.item(), output_dir=local_output_dir
)

pdal_pipeline

{'pipeline': [{'type': 'readers.ept',
   'filename': 'https://s3-us-west-2.amazonaws.com/usgs-lidar-public/CO_CentralEasternPlains_1_2020/ept.json',
   'bounds': '(([-11410804.4037645068, -11409134.6114026085], [4818825.9525320893, 4820272.3693662276]))'},
  {'type': 'filters.crop',
   'polygon': 'POLYGON ((-11409134.611402608 4818825.952532089, -11409134.611402608 4820272.369366228, -11410804.403764507 4820272.369366228, -11410804.403764507 4818825.952532089, -11409134.611402608 4818825.952532089))'},
  {'type': 'writers.las',
   'filename': 'CO_CentralEasternPlains_1_2020_EPT_subset_pipeline.laz',
   'compression': 'laszip'}]}