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/stable/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/stable/.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... | POLYGON ((-102.60871 39.70725, -102.6084 39.71... | 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"],
)
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.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");
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 199 ms, sys: 31.6 ms, total: 231 ms
Wall time: 11.9 s
/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/stable/src/coincident/io/xarray.py:480: FutureWarning: In a future version of xarray the default value for join will change from join='outer' to join='exact'. This change will result in the following ValueError: cannot be aligned with join='exact' because index/labels/sizes are not equal along these coordinates (dimensions): 'x' ('x',) The recommendation is to set join explicitly for this case.
merged = xr.concat(dem_dataarrays, dim="band").max(dim="band").rename("elevation")
/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/stable/src/coincident/io/xarray.py:480: FutureWarning: In a future version of xarray the default value for join will change from join='outer' to join='exact'. This change will result in the following ValueError: cannot be aligned with join='exact' because index/labels/sizes are not equal along these coordinates (dimensions): 'y' ('y',) The recommendation is to set join explicitly for this case.
merged = xr.concat(dem_dataarrays, dim="band").max(dim="band").rename("elevation")
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:
* y (y) float64 9kB 4.395e+06 4.395e+06 ... 4.396e+06 4.396e+06
* x (x) float64 11kB 7.14e+05 7.14e+05 ... 7.153e+05 7.153e+05
spatial_ref int64 8B 0
Attributes:
AREA_OR_POINT: Area
scale_factor: 1.0
add_offset: 0.0usgs_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 79.7 ms, sys: 7.79 ms, total: 87.4 ms
Wall time: 4.06 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: 100da_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... | POLYGON ((-102.60871 39.70725, -102.6084 39.71... | 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 44.6 ms, sys: 4.61 ms, total: 49.2 ms
Wall time: 184 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_713000_4395000_classified_po... | https://storage.googleapis.com/neon-aop-produc... | POLYGON ((-102.50479 39.67804, -102.50447 39.6... |
| 2 | NEON_D10_ARIK_DP1_713000_4394000_classified_po... | https://storage.googleapis.com/neon-aop-produc... | POLYGON ((-102.50511 39.66904, -102.50479 39.6... |
| 3 | NEON_D10_ARIK_DP1_714000_4395000_classified_po... | https://storage.googleapis.com/neon-aop-produc... | POLYGON ((-102.49314 39.67779, -102.49281 39.6... |
| 4 | NEON_D10_ARIK_DP1_715000_4394000_classified_po... | https://storage.googleapis.com/neon-aop-produc... | POLYGON ((-102.48182 39.66854, -102.48149 39.6... |
%%time
gf_usgs_lpc_tiles = coincident.io.download.fetch_lpc_tiles(
aoi=aoi, dataset_id=usgs_project, provider="USGS"
)
CPU times: user 19.3 ms, sys: 1.04 ms, total: 20.3 ms
Wall time: 1.64 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
Now, we can download the laz files
coincident.io.download.download_files(
gf_neon_lpc_tiles["url"], output_dir=local_output_dir
)
['/tmp/NEON_D10_ARIK_DP1_715000_4395000_classified_point_cloud_colorized.laz',
'/tmp/NEON_D10_ARIK_DP1_713000_4395000_classified_point_cloud_colorized.laz',
'/tmp/NEON_D10_ARIK_DP1_713000_4394000_classified_point_cloud_colorized.laz',
'/tmp/NEON_D10_ARIK_DP1_714000_4395000_classified_point_cloud_colorized.laz',
'/tmp/NEON_D10_ARIK_DP1_715000_4394000_classified_point_cloud_colorized.laz',
'/tmp/NEON_D10_ARIK_DP1_714000_4394000_classified_point_cloud_colorized.laz']
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'}]}