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.

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/v0.3/src/coincident/io/download.py:26: 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...	http://prd-tnm.s3.amazonaws.com/index.html?pre...	http://prd-tnm.s3.amazonaws.com/index.html?pre...	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        http://prd-tnm.s3.amazonaws.com/index.html?pre...
metadata_link         http://prd-tnm.s3.amazonaws.com/index.html?pre...
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/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/geopandas/array.py:403: 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
0	ARIK	Arikaree River NEON	2020-06-01	2020-06-30	https://data.neonscience.org/api/v0/data/DP3.3...	POLYGON ((-102.60902 39.69825, -102.60871 39.7...

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"],
)

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()

# 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)
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");

Actually read in the DEMs

datetime_str = gf_neon.end_datetime.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"
)

/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/rioxarray/_io.py:1146: RuntimeWarning: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order
  if riods.subdatasets:
/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/rioxarray/_io.py:1146: RuntimeWarning: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order
  if riods.subdatasets:

/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/rioxarray/_io.py:1146: RuntimeWarning: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order
  if riods.subdatasets:
/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/rioxarray/_io.py:1146: RuntimeWarning: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order
  if riods.subdatasets:

/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/rioxarray/_io.py:1146: RuntimeWarning: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order
  if riods.subdatasets:
/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/rioxarray/_io.py:1146: RuntimeWarning: TIFFReadDirectoryCheckOrder:Invalid TIFF directory; tags are not sorted in ascending order
  if riods.subdatasets:

CPU times: user 186 ms, sys: 38 ms, total: 224 ms
Wall time: 1.66 s

da_neon_dem

<xarray.DataArray 'elevation' (y: 1146, x: 1318)> Size: 6MB
dask.array<getitem, shape=(1146, 1318), dtype=float32, chunksize=(808, 1000), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 11kB 7.14e+05 7.14e+05 ... 7.153e+05 7.153e+05
  * y            (y) float64 9kB 4.395e+06 4.395e+06 ... 4.396e+06 4.396e+06
    spatial_ref  int64 8B 0

xarray.DataArray

'elevation'

• y: 1146
• x: 1318

• dask.array<chunksize=(338, 40), meta=np.ndarray>

  Array Chunk Bytes 5.76 MiB 3.08 MiB Shape (1146, 1318) (808, 1000) Dask graph 6 chunks in 64 graph layers Data type float32 numpy.ndarray

• Coordinates: (3)
  • 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, ..., 715275.5, 715276.5, 715277.5],
          shape=(1318,))

  • y
    (y)
    float64
    4.395e+06 4.395e+06 ... 4.396e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([4394662.5, 4394663.5, 4394664.5, ..., 4395805.5, 4395806.5, 4395807.5],
          shape=(1146,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],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","32613"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / 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["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],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","32613"]]

    GeoTransform :

    713960.0 1.0 0.0 4394662.0 0.0 1.0

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([713960.5, 713961.5, 713962.5, 713963.5, 713964.5, 713965.5, 713966.5,
           713967.5, 713968.5, 713969.5,
           ...
           715268.5, 715269.5, 715270.5, 715271.5, 715272.5, 715273.5, 715274.5,
           715275.5, 715276.5, 715277.5],
          dtype='float64', name='x', length=1318))

  • y
    PandasIndex

    PandasIndex(Index([4394662.5, 4394663.5, 4394664.5, 4394665.5, 4394666.5, 4394667.5,
           4394668.5, 4394669.5, 4394670.5, 4394671.5,
           ...
           4395798.5, 4395799.5, 4395800.5, 4395801.5, 4395802.5, 4395803.5,
           4395804.5, 4395805.5, 4395806.5, 4395807.5],
          dtype='float64', name='y', length=1146))

• Attributes: (0)

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 553 ms, sys: 245 ms, total: 797 ms
Wall time: 2.62 s

da_usgs_dem

<xarray.DataArray 'elevation' (y: 1146, x: 1317)> Size: 6MB
dask.array<getitem, shape=(1146, 1317), dtype=float32, chunksize=(1146, 1317), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 11kB 7.14e+05 7.14e+05 ... 7.153e+05 7.153e+05
  * y            (y) float64 9kB 4.396e+06 4.396e+06 ... 4.395e+06 4.395e+06
    spatial_ref  int64 8B 0

xarray.DataArray

'elevation'

• y: 1146
• x: 1317

• dask.array<chunksize=(1146, 1317), meta=np.ndarray>

  Array Chunk Bytes 5.76 MiB 5.76 MiB Shape (1146, 1317) (1146, 1317) Dask graph 1 chunks in 8 graph layers Data type float32 numpy.ndarray

• Coordinates: (3)
  • 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([713961.499973, 713962.499973, 713963.499973, ..., 715275.499973,
           715276.499973, 715277.499973], shape=(1317,))

  • 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.500046, 4395806.500046, 4395805.500046, ..., 4394664.500046,
           4394663.500046, 4394662.500046], shape=(1146,))

  • 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.999973323 1.0 0.0 4395808.000045877 0.0 -1.0

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([713961.499973323, 713962.499973323, 713963.499973323, 713964.499973323,
           713965.499973323, 713966.499973323, 713967.499973323, 713968.499973323,
           713969.499973323, 713970.499973323,
           ...
           715268.499973323, 715269.499973323, 715270.499973323, 715271.499973323,
           715272.499973323, 715273.499973323, 715274.499973323, 715275.499973323,
           715276.499973323, 715277.499973323],
          dtype='float64', name='x', length=1317))

  • y
    PandasIndex

    PandasIndex(Index([4395807.500045877, 4395806.500045877, 4395805.500045877,
           4395804.500045877, 4395803.500045877, 4395802.500045877,
           4395801.500045877, 4395800.500045877, 4395799.500045877,
           4395798.500045877,
           ...
           4394671.500045877, 4394670.500045877, 4394669.500045877,
           4394668.500045877, 4394667.500045877, 4394666.500045877,
           4394665.500045877, 4394664.500045877, 4394663.500045877,
           4394662.500045877],
          dtype='float64', name='y', length=1146))

• Attributes: (0)

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

Download

Note

coincident.io.download.download_neon_dem needs the NEON site’s end_datetime, not start_datetime (separated by 1 month) to work

gf_neon

	id	title	start_datetime	end_datetime	product_url	geometry
0	ARIK	Arikaree River NEON	2020-06-01	2020-06-30	https://data.neonscience.org/api/v0/data/DP3.3...	POLYGON ((-102.60902 39.69825, -102.60871 39.7...

local_output_dir = "/tmp"

coincident.io.download.download_neon_dem(
    aoi=aoi,
    datetime_str=gf_neon.end_datetime.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.

%%time
gf_neon_lpc_tiles = coincident.io.download.fetch_neon_lpc_tiles(
    aoi=aoi, datetime_str=gf_neon.start_datetime.item(), site_id=gf_neon.id.item()
)

CPU times: user 55.7 ms, sys: 2.93 ms, total: 58.7 ms
Wall time: 328 ms

gf_neon_lpc_tiles.head()

	name	url	geometry
0	NEON_D10_ARIK_DP1_715000_4395000_classified_po...	https://storage.googleapis.com/neon-aop-produc...	POLYGON ((-102.48149 39.67754, -102.48116 39.6...
1	NEON_D10_ARIK_DP1_714000_4394000_classified_po...	https://storage.googleapis.com/neon-aop-produc...	POLYGON ((-102.49347 39.66879, -102.49314 39.6...
2	NEON_D10_ARIK_DP1_713000_4395000_classified_po...	https://storage.googleapis.com/neon-aop-produc...	POLYGON ((-102.50479 39.67804, -102.50447 39.6...
3	NEON_D10_ARIK_DP1_713000_4394000_classified_po...	https://storage.googleapis.com/neon-aop-produc...	POLYGON ((-102.50511 39.66904, -102.50479 39.6...
4	NEON_D10_ARIK_DP1_714000_4395000_classified_po...	https://storage.googleapis.com/neon-aop-produc...	POLYGON ((-102.49314 39.67779, -102.49281 39.6...

%%time
gf_usgs_lpc_tiles = coincident.io.download.fetch_usgs_lpc_tiles(
    aoi=aoi, project=usgs_project
)

CPU times: user 8.7 ms, sys: 1.24 ms, total: 9.94 ms
Wall time: 1.2 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)

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

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'}]}

NCALM

Search

aoi = gpd.read_file(
    "https://raw.githubusercontent.com/unitedstates/districts/refs/heads/gh-pages/states/WA/shape.geojson"
)
aoi.plot();

gf_ncalm = coincident.search.search(
    dataset="ncalm", intersects=aoi, datetime=["2018-09-19"]
)

gf_ncalm

	id	name	title	start_datetime	end_datetime	geometry
0	OTLAS.072019.6339.1	WA18_Wall	High-Resolution Mapping of Goat Rock Volcano, WA	2018-09-19	2018-09-20	POLYGON ((-121.46701 46.48376, -121.45914 46.4...

Now, let’s subset to a small AOI for convenience sake

buffer_size = 0.01
centroid = gf_ncalm.centroid
mini_aoi = gpd.GeoDataFrame(
    geometry=[
        box(
            centroid.x - buffer_size,
            centroid.y - buffer_size,
            centroid.x + buffer_size,
            centroid.y + buffer_size,
        )
    ],
    crs="EPSG:4326",
)

/tmp/ipykernel_1059/991173689.py:2: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  centroid = gf_ncalm.centroid
/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/shapely/geometry/polygon.py:91: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead
  return [float(c) for c in o]

m = gf_ncalm.explore()
mini_aoi.explore(m=m, color="red")

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

Note

The NCALM DEMs are not tiled, so reading in the data and downloading takes a longer time

Subset data

%%time
da_ncalm_dtm = coincident.io.xarray.load_ncalm_dem(
    aoi=mini_aoi, product="dtm", dataset_id=gf_ncalm["name"].item()
)

CPU times: user 6.29 s, sys: 2.25 s, total: 8.54 s
Wall time: 22.3 s

da_ncalm_dtm

<xarray.DataArray 'elevation' (y: 2254, x: 1579)> Size: 14MB
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]],
      shape=(2254, 1579), dtype=float32)
Coordinates:
  * x            (x) float64 13kB 6.226e+05 6.226e+05 ... 6.242e+05 6.242e+05
  * y            (y) float64 18kB 5.152e+06 5.152e+06 ... 5.15e+06 5.15e+06
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

xarray.DataArray

'elevation'

• y: 2254
• x: 1579

• nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan

  array([[nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         ...,
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan]],
        shape=(2254, 1579), dtype=float32)

• Coordinates: (3)
  • x
    (x)
    float64
    6.226e+05 6.226e+05 ... 6.242e+05

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([622599.000001, 622600.000001, 622601.000001, ..., 624175.000001,
           624176.000001, 624177.000001], shape=(1579,))

  • y
    (y)
    float64
    5.152e+06 5.152e+06 ... 5.15e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([5152454.000016, 5152453.000016, 5152452.000016, ..., 5150203.000016,
           5150202.000016, 5150201.000016], shape=(2254,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 10N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -123.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]

    GeoTransform :

    622598.5000010493 1.0000000000086315 0.0 5152454.500016331 0.0 -1.0000000000086808

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([622599.0000010493, 622600.0000010493, 622601.0000010493,
           622602.0000010493, 622603.0000010494, 622604.0000010494,
           622605.0000010494, 622606.0000010494, 622607.0000010494,
           622608.0000010494,
           ...
           624168.0000010629, 624169.0000010629, 624170.0000010629,
           624171.0000010629, 624172.0000010629, 624173.0000010629,
           624174.0000010629, 624175.0000010629, 624176.0000010629,
           624177.0000010629],
          dtype='float64', name='x', length=1579))

  • y
    PandasIndex

    PandasIndex(Index([5152454.000016331, 5152453.000016331, 5152452.000016331,
           5152451.000016331, 5152450.000016331, 5152449.000016331,
           5152448.000016331, 5152447.000016331, 5152446.000016331,
            5152445.00001633,
           ...
           5150210.000016311, 5150209.000016311, 5150208.000016311,
           5150207.000016311, 5150206.000016311, 5150205.000016311,
           5150204.000016311, 5150203.000016311, 5150202.000016311,
           5150201.000016311],
          dtype='float64', name='y', length=2254))

• Attributes: (3)

  AREA_OR_POINT :

  Area

  scale_factor :

  1.0

  add_offset :

  0.0

da_ncalm_dtm.coarsen(x=20, y=20, boundary="trim").mean().plot.imshow();

Download

%%time
coincident.io.download.download_ncalm_dem(
    aoi=mini_aoi,
    dataset_id=gf_ncalm["name"].item(),
    product="dtm",
    output_dir=local_output_dir,
)

CPU times: user 6.35 s, sys: 1.75 s, total: 8.1 s
Wall time: 21.3 s

%%time
gf_ncalm_lpc_tiles = coincident.io.download.fetch_ncalm_lpc_tiles(
    aoi=mini_aoi, dataset_name=gf_ncalm.name.item()
)

CPU times: user 69.4 ms, sys: 9.78 ms, total: 79.2 ms
Wall time: 467 ms

gf_ncalm_lpc_tiles.head()

	name	url	geometry
0	WA18_Wall/622000_5150000.laz	https://opentopography.s3.sdsc.edu/pc-bulk/WA1...	POLYGON ((-121.39723 46.49234, -121.39697 46.5...
1	WA18_Wall/622000_5151000.laz	https://opentopography.s3.sdsc.edu/pc-bulk/WA1...	POLYGON ((-121.39697 46.50133, -121.3967 46.51...
2	WA18_Wall/622000_5152000.laz	https://opentopography.s3.sdsc.edu/pc-bulk/WA1...	POLYGON ((-121.3967 46.51033, -121.39644 46.51...
3	WA18_Wall/623000_5150000.laz	https://opentopography.s3.sdsc.edu/pc-bulk/WA1...	POLYGON ((-121.3842 46.49216, -121.38394 46.50...
4	WA18_Wall/623000_5151000.laz	https://opentopography.s3.sdsc.edu/pc-bulk/WA1...	POLYGON ((-121.38394 46.50115, -121.38367 46.5...

%%time
gf_ncalm_lpc_tiles = coincident.io.download.fetch_ncalm_lpc_tiles(
    aoi=mini_aoi, dataset_name="WA18_Wall"
)

CPU times: user 72.3 ms, sys: 14.1 ms, total: 86.3 ms
Wall time: 463 ms

gf_ncalm.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

m = gf_ncalm.explore()
gf_ncalm_lpc_tiles.explore(m=m, color="red")
mini_aoi.explore(m=m, color="black")

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

NOAA

Search

aoi = gpd.read_file(
    "https://raw.githubusercontent.com/unitedstates/districts/refs/heads/gh-pages/states/FL/shape.geojson"
)
gf_noaa = coincident.search.search(
    dataset="noaa", intersects=aoi, datetime=["2022-10-27"]
)

---------------------------------------------------------------------------
HTTPError                                 Traceback (most recent call last)
Cell In[43], line 1
----> 1 aoi = gpd.read_file(
      2     "https://raw.githubusercontent.com/unitedstates/districts/refs/heads/gh-pages/states/FL/shape.geojson"
      3 )
      4 gf_noaa = coincident.search.search(
      5     dataset="noaa", intersects=aoi, datetime=["2022-10-27"]
      6 )

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/site-packages/geopandas/io/file.py:288, in _read_file(filename, bbox, mask, columns, rows, engine, **kwargs)
    282 from_bytes = False
    283 if _is_url(filename):
    284     # if it is a url that supports random access -> pass through to
    285     # pyogrio/fiona as is (to support downloading only part of the file)
    286     # otherwise still download manually because pyogrio/fiona don't support
    287     # all types of urls (https://github.com/geopandas/geopandas/issues/2908)
--> 288     with urllib.request.urlopen(filename) as response:
    289         if not response.headers.get("Accept-Ranges") == "bytes":
    290             filename = response.read()

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/urllib/request.py:215, in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    213 else:
    214     opener = _opener
--> 215 return opener.open(url, data, timeout)

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/urllib/request.py:521, in OpenerDirector.open(self, fullurl, data, timeout)
    519 for processor in self.process_response.get(protocol, []):
    520     meth = getattr(processor, meth_name)
--> 521     response = meth(req, response)
    523 return response

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/urllib/request.py:630, in HTTPErrorProcessor.http_response(self, request, response)
    627 # According to RFC 2616, "2xx" code indicates that the client's
    628 # request was successfully received, understood, and accepted.
    629 if not (200 <= code < 300):
--> 630     response = self.parent.error(
    631         'http', request, response, code, msg, hdrs)
    633 return response

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/urllib/request.py:559, in OpenerDirector.error(self, proto, *args)
    557 if http_err:
    558     args = (dict, 'default', 'http_error_default') + orig_args
--> 559     return self._call_chain(*args)

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/urllib/request.py:492, in OpenerDirector._call_chain(self, chain, kind, meth_name, *args)
    490 for handler in handlers:
    491     func = getattr(handler, meth_name)
--> 492     result = func(*args)
    493     if result is not None:
    494         return result

File ~/checkouts/readthedocs.org/user_builds/coincident/checkouts/v0.3/.pixi/envs/dev/lib/python3.12/urllib/request.py:639, in HTTPDefaultErrorHandler.http_error_default(self, req, fp, code, msg, hdrs)
    638 def http_error_default(self, req, fp, code, msg, hdrs):
--> 639     raise HTTPError(req.full_url, code, msg, hdrs, fp)

HTTPError: HTTP Error 429: Too Many Requests

buffer_size = 0.02
centroid = gf_noaa.centroid
mini_aoi = gpd.GeoDataFrame(
    geometry=[
        box(
            centroid.x - buffer_size,
            centroid.y - buffer_size,
            centroid.x + buffer_size,
            centroid.y + buffer_size,
        )
    ],
    crs="EPSG:4326",
)

m = gf_noaa.explore()
mini_aoi.explore(m=m, color="red")

the name and id being identical is expected

gf_noaa

Note

Our coincident.search.search(dataset=”noaa”) returns the dataset ids “Lidar Datasets at NOAA Digital Coast” whereas coincident.io.xarray.load_noaa_dem() requires the ids from the “Imagery and Elevation Raster Datasets at NOAA Digital Coast” dataset. The corresponding elevation raster dataset id is the same as the lidar dataset id + 1. e.g. “Great Bay NERR UAS Lidar” has id 10175 for lidar data and id 10176 for dem data

noaa_dem_id = int(gf_noaa.id.item()) + 1
print(f"NOAA LiDAR id: {gf_noaa.id.item()}  NOAA DEM id: {noaa_dem_id}")

Warning

The larger the NOAA flight, the longer the below function takes regardless of your input AOI

Subset

%%time
da_noaa_dem = coincident.io.xarray.load_noaa_dem(mini_aoi, noaa_dem_id)

da_noaa_dem

da_noaa_dem.coarsen(x=50, y=50, boundary="trim").mean().plot.imshow();

buffer_size = 0.008
centroid = gf_noaa.centroid
mini_aoi = gpd.GeoDataFrame(
    geometry=[
        box(
            centroid.x - buffer_size,
            centroid.y - buffer_size,
            centroid.x + buffer_size,
            centroid.y + buffer_size,
        )
    ],
    crs="EPSG:4326",
)
m = gf_noaa.explore()
mini_aoi.explore(m=m, color="red")

Download

%%time
coincident.io.download.download_noaa_dem(
    aoi=mini_aoi, dataset_id=noaa_dem_id, output_dir=local_output_dir
)

list(mini_aoi.bounds.values)

gf_noaa.id.item()

%%time
gf_noaa_lpc_tiles = coincident.io.download.fetch_noaa_lpc_tiles(
    aoi=mini_aoi, dataset_id=gf_noaa.id.item()
)

gf_noaa_lpc_tiles

m = gf_noaa.explore()
gf_noaa_lpc_tiles.explore(m=m, color="red")