Quickstart#

This library provides functions to hone in on coincident datasets and visualize their coverage.

The current recommended workflow is to:

  1. start with a lidar dataset

  2. reduce to region with coincident maxar stereo within an acceptable temporal range

  3. optionally reduce further with additional datasets such as icesat-2 and gedi altimetry

This notebook provides an example starting from USGS 3DEP LiDAR in Colorado, USA

import geopandas as gpd
import coincident

print(coincident.__version__)

/home/docs/checkouts/readthedocs.org/user_builds/coincident/checkouts/104/src/coincident/io/download.py:27: 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

0.5.dev4+gc0bf3bd.d20250804

Define an area of interest#

Searches often start with an spatial subset, in this notebook we know we are interested in datasets in Colorado. We recommend restricting searches to areas at the ‘State’ scale rather than Country or Global scales

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

aoi.explore(color="black", style_kwds=dict(fill=False))
Make this Notebook Trusted to load map: File -> Trust Notebook

Note

this polygon for the state of Colorado is simple and has only 4 vertices in the corners, but it’s good practice good to check the number of vertices you have before searching. The simpler your search polygons are the faster your searches will be!

print("Original number of vertices:", aoi.count_coordinates().iloc[0])
aoi = aoi.simplify(0.01)
print("Simplified number of vertices:", aoi.count_coordinates().iloc[0])

Original number of vertices: 260
Simplified number of vertices: 7

Uniform search method#

the coincident package provides a search() method that has the same syntax regardless of which dataset you are searching. Behind the scenes, polygons intersecting your area of interest are efficiently located and returned as a geodataframe.

For 3DEP LiDAR, we start by searching bounding boxes for each ‘workunit’. Once we identify a workunit, we load a precise polygon delineating the extent of LiDAR measurements.

The search function is essentially a pystac_client.ItemSearch with an extra argument dataset.

coincident.datasets.aliases

['maxar',
 '3dep',
 'icesat-2',
 'gedi',
 'gliht',
 'ablvis2_1',
 'aflvis2_1',
 'cop30',
 'worldcover',
 'tdx',
 'noaa',
 'ncalm',
 'neon']

gf = coincident.search.search(
    dataset="3dep",
    intersects=aoi,
    datetime=["2018", "2024"],
)
gf.explore(column="workunit", popup=True)
Make this Notebook Trusted to load map: File -> Trust Notebook
# From this search we identify a specific lidar acquisition of interest
gf = gf[gf.workunit == "CO_WestCentral_2019"]
gf
minx maxx miny maxy geometry workunit workunit_id project project_id start_datetime ... onemeter_reason seamless_category seamless_reason lpc_link sourcedem_link metadata_link collection datetime dayofyear duration
fid
2446 -107.989311 -106.135696 38.134274 39.389549 POLYGON ((-106.1357 38.13427, -106.1357 39.389... CO_WestCentral_2019 175984 CO_WestCentral_2019_A19 175987 2019-08-21 ... Meets 3DEP 1-m DEM requirements 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... 3DEP 2019-09-04 12:00:00 247 29

1 rows × 37 columns

Each Workunit has a unique ‘Feature ID’, or ‘FID’ that can be used to efficiently retrieve the full-resultion detailed MultiPolygon footprint from the USGS’s WESM.gpkg file in AWS

gf.index.values

array([2446])

USGS 3DEP Lidar#

Some datasets have additional functions available to load auxiliary data. For example, we can load original high resolution polygon. In addition to loading “swath” polygons to understand exact days when acquisitions were made.

Warning

Swath polygons are not available for all LiDAR workunits. They are more likely to exist for acquisitions after 2018.

gf_wesm = coincident.search.wesm.load_by_fid(fids=gf.index.values)

gf_wesm.explore(column="workunit", popup=True)
Make this Notebook Trusted to load map: File -> Trust Notebook
gf_wesm[["workunit", "start_datetime", "end_datetime", "duration"]]  # duration in days
workunit start_datetime end_datetime duration
0 CO_WestCentral_2019 2019-08-21 2019-09-19 29

Within this ‘workunit’ data was collected over 29 days, in order to see when data was collected within this polygon we need to load a corresponding ‘swath’ polygon. coincident provides a helper function for this, which loads the swath polygon for a given workunit

Note

Swath polygons have detailed timing information for individual LiDAR flight lines composing a given ‘workunit’. Not all workunits have swath polygons. They tend to be available for data collected after 2019.

# NOTE: be patient here, reading shapefiles from s3 can be slow depending on your bandwidth
gf_swath = coincident.search.wesm.get_swath_polygons(gf.workunit.iloc[0])

print("Number of swath polygons:", len(gf_swath))

Number of swath polygons: 314

# Plot them, simplify first for faster plotting
gf_swath["geometry"] = gf_swath.simplify(0.01)
gf_swath.explore(column="dayofyear", cmap="plasma")
Make this Notebook Trusted to load map: File -> Trust Notebook

As you can see in the above figure, the actual day of observation for any point on the ground can vary in complex ways based on the flight paths taken during the LiDAR collection period.

GEDI#

Search for Coincident GEDI L2A.

# We've refined our search polygon again, so use a new AOI
aoi = gpd.GeoSeries(gf_i.union_all().convex_hull, crs="EPSG:4326")

gf_gedi = coincident.search.search(
    dataset="gedi",
    intersects=aoi,
    datetime=date_range,
)

print(len(gf_gedi))
gf_gedi.head()

9
assets bbox collection geometry id links stac_extensions stac_version type datetime end_datetime start_datetime storage:schemes dayofyear
0 {'002/GEDI02_A_2019221181154_O03725_02_T03077_... {'xmin': -130.7824147, 'ymin': 14.5375478, 'xm... GEDI02_A_002 POLYGON ((-130.74077 14.53755, -128.39101 17.5... GEDI02_A_2019221181154_O03725_02_T03077_02_003... [{'href': 'https://cmr.earthdata.nasa.gov/stac... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-08-09 18:11:54+00:00 2019-08-09 19:44:41+00:00 2019-08-09 18:11:54+00:00 {'aws': {'bucket': 'lp-prod-protected', 'platf... 221
1 {'002/GEDI02_A_2019228215125_O03836_03_T03938_... {'xmin': -154.5108058, 'ymin': 0.0764856, 'xma... GEDI02_A_002 POLYGON ((-154.51081 51.7766, -148.53387 51.60... GEDI02_A_2019228215125_O03836_03_T03938_02_003... [{'href': 'https://cmr.earthdata.nasa.gov/stac... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-08-16 21:51:25+00:00 2019-08-16 23:24:15+00:00 2019-08-16 21:51:25+00:00 {'aws': {'bucket': 'lp-prod-protected', 'platf... 228
2 {'002/GEDI02_A_2019232140204_O03893_02_T00583_... {'xmin': -131.649277, 'ymin': 15.7666963, 'xma... GEDI02_A_002 POLYGON ((-131.60759 15.7667, -129.21939 18.73... GEDI02_A_2019232140204_O03893_02_T00583_02_003... [{'href': 'https://cmr.earthdata.nasa.gov/stac... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-08-20 14:02:04+00:00 2019-08-20 15:34:53+00:00 2019-08-20 14:02:04+00:00 {'aws': {'bucket': 'lp-prod-protected', 'platf... 232
3 {'002/GEDI02_A_2019247142840_O04126_03_T02362_... {'xmin': -154.8002265, 'ymin': -0.0163617, 'xm... GEDI02_A_002 POLYGON ((-154.80023 51.77634, -148.81565 51.5... GEDI02_A_2019247142840_O04126_03_T02362_02_003... [{'href': 'https://cmr.earthdata.nasa.gov/stac... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-09-04 14:28:40+00:00 2019-09-04 16:01:29+00:00 2019-09-04 14:28:40+00:00 {'aws': {'bucket': 'lp-prod-protected', 'platf... 247
4 {'002/GEDI02_A_2019251125026_O04187_03_T01704_... {'xmin': -153.6270991, 'ymin': 0.0029707, 'xma... GEDI02_A_002 POLYGON ((-153.6271 51.77622, -147.64024 51.60... GEDI02_A_2019251125026_O04187_03_T01704_02_003... [{'href': 'https://cmr.earthdata.nasa.gov/stac... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-09-08 12:50:26+00:00 2019-09-08 14:23:15+00:00 2019-09-08 12:50:26+00:00 {'aws': {'bucket': 'lp-prod-protected', 'platf... 251 
gf_gedi.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook

Altimeter granules span a large geographic area! So let’s again cut results down to the area of intersection.

# Normalize colormap across both dataframes
vmin, vmax = gpd.pd.concat([gf_i.dayofyear, gf_gedi.dayofyear]).agg(["min", "max"])
cmap = "plasma"

# NOTE: here we just take the boundary of the union of all maxar+lidar regions to avoid many overlay geometries
union = gpd.GeoDataFrame(geometry=[gf_i.union_all()], crs="EPSG:4326")
m = gf_i.explore(column="dayofyear", cmap=cmap, vmin=vmin, vmax=vmax)
gf_gedi = union.overlay(gf_gedi, how="intersection")
gf_gedi.explore(m=m, column="dayofyear", cmap=cmap, vmin=vmin, vmax=vmax, legend=False)
m
Make this Notebook Trusted to load map: File -> Trust Notebook

Immediately we see some potential for close-in time acquisitions.

Note

As an additional step, you might want to add another filtering step to remove tiny polygons

# For each Stereo pair, describe number of altimeter passes and day offset
stereo_pair = gf_i.iloc[[0]]
gf_dt = stereo_pair.overlay(gf_gedi, how="intersection")
m = stereo_pair.explore(
    column="dayofyear", vmin=vmin, vmax=vmax, cmap=cmap, legend=False
)
gf_dt.explore(m=m, column="dayofyear_2", vmin=vmin, vmax=vmax, cmap=cmap)
Make this Notebook Trusted to load map: File -> Trust Notebook
print("Maxar Stereo Acquisition DOY - GEDI DOY:")
print(stereo_pair.stereo_pair_identifiers.values[0])
(gf_dt.dayofyear_1 - gf_dt.dayofyear_2).describe()

Maxar Stereo Acquisition DOY - GEDI DOY:
['11a11283-6661-4f91-9bcb-5dab3c6a5d02-inv']

count     4.000000
mean     23.000000
std      17.492856
min       3.000000
25%      14.250000
50%      22.000000
75%      30.750000
max      45.000000
dtype: float64

# Consider just the minimum offset for each stereo pair
# NOTE: probably a fancier way to do this using a multiindex/groupby
print("Maxar Stereo Acquisition DOY - GEDI DOY")
for i in range(len(gf_i)):
    stereo_pair = gf_i.iloc[[i]]
    gf_dt = stereo_pair.overlay(gf_gedi, how="intersection")
    min_dt = (gf_dt.dayofyear_1 - gf_dt.dayofyear_2).min()
    print(stereo_pair.stereo_pair_identifiers.values[0], min_dt)

Maxar Stereo Acquisition DOY - GEDI DOY
['11a11283-6661-4f91-9bcb-5dab3c6a5d02-inv'] 3
['21d1d4b2-47a4-4c4d-a26f-5c688199b518-inv'] -36
['28b60c45-3164-4cd0-9b43-67dac4e8c7b4-inv'] -31
['41fef87b-0d7c-46c2-a598-bcd178dfb9e0-inv'] -29
['6002024a-ca85-4f3c-8c7f-ce849c175a52-inv'] 3

Note

In this case, we have some altimeter acquisitions that were just 3 days after a maxar stereo acquisition!

ICESat-2#

Search for Coincident ICESat-2 Altimetry (ATL03)

gf_is2 = coincident.search.search(
    dataset="icesat-2",
    intersects=aoi,
    datetime=date_range,
)
print(len(gf_is2))
gf_is2

5
assets bbox collection geometry id links stac_extensions stac_version type datetime end_datetime start_datetime storage:schemes dayofyear
0 {'browse': {'href': 'https://data.nsidc.earthd... {'xmin': -109.90889, 'ymin': 26.94277, 'xmax':... ATL03_006 POLYGON ((-109.68811 59.53895, -109.90889 59.5... ATL03_20190807082015_06150402_006_02.h5 [{'href': 'https://cmr.earthdata.nasa.gov/clou... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-08-07 08:20:14.997000+00:00 2019-08-07 08:28:45.648000+00:00 2019-08-07 08:20:14.997000+00:00 {'aws': {'bucket': 'nsidc-cumulus-prod-protect... 219
1 {'browse': {'href': 'https://data.nsidc.earthd... {'xmin': -108.09761, 'ymin': 26.96544, 'xmax':... ATL03_006 POLYGON ((-103.41614 59.54107, -103.63704 59.5... ATL03_20190903191210_10340406_006_02.h5 [{'href': 'https://cmr.earthdata.nasa.gov/clou... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-09-03 19:12:09.748000+00:00 2019-09-03 19:20:40.343000+00:00 2019-09-03 19:12:09.748000+00:00 {'aws': {'bucket': 'nsidc-cumulus-prod-protect... 246
2 {'browse': {'href': 'https://data.nsidc.earthd... {'xmin': -109.64978, 'ymin': 26.94274, 'xmax':... ATL03_006 POLYGON ((-109.42923 59.53777, -109.64978 59.5... ATL03_20190905065620_10570402_006_02.h5 [{'href': 'https://cmr.earthdata.nasa.gov/clou... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-09-05 06:56:19.240000+00:00 2019-09-05 07:04:49.503000+00:00 2019-09-05 06:56:19.240000+00:00 {'aws': {'bucket': 'nsidc-cumulus-prod-protect... 248
3 {'browse': {'href': 'https://data.nsidc.earthd... {'xmin': -110.29975, 'ymin': 26.95232, 'xmax':... ATL03_006 POLYGON ((-110.07845 59.54577, -110.29975 59.5... ATL03_20190909064801_11180402_006_02.h5 [{'href': 'https://cmr.earthdata.nasa.gov/clou... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-09-09 06:48:00.991000+00:00 2019-09-09 06:56:31.581000+00:00 2019-09-09 06:48:00.991000+00:00 {'aws': {'bucket': 'nsidc-cumulus-prod-protect... 252
4 {'browse': {'href': 'https://data.nsidc.earthd... {'xmin': -107.83927, 'ymin': 26.96543, 'xmax':... ATL03_006 POLYGON ((-103.15539 59.54093, -103.37634 59.5... ATL03_20191002174817_00890506_006_02.h5 [{'href': 'https://cmr.earthdata.nasa.gov/clou... [https://stac-extensions.github.io/storage/v2.... 1.1.0 Feature 2019-10-02 17:48:16.589000+00:00 2019-10-02 17:56:47.166000+00:00 2019-10-02 17:48:16.589000+00:00 {'aws': {'bucket': 'nsidc-cumulus-prod-protect... 275 
# Normalize colormap across both dataframes
vmin, vmax = gpd.pd.concat([gf_i.dayofyear, gf_is2.dayofyear]).agg(["min", "max"])
cmap = "plasma"

# NOTE: here we just take the boundary of the union of all maxar+lidar regions to avoid many overlay geometries
union = gpd.GeoDataFrame(geometry=[gf_i.union_all()], crs="EPSG:4326")
m = gf_i.explore(
    column="dayofyear", cmap=cmap, vmin=vmin, vmax=vmax, style_kwds=dict(color="black")
)
gf_is2 = union.overlay(gf_is2, how="intersection")
gf_is2.explore(m=m, column="dayofyear", cmap=cmap, vmin=vmin, vmax=vmax, legend=False)
m
Make this Notebook Trusted to load map: File -> Trust Notebook
# For each Stereo pair, describe number of altimeter passes and day offset
stereo_pair = gf_i.iloc[[0]]
gf_dt = stereo_pair.overlay(gf_is2, how="intersection")
m = stereo_pair.explore(
    column="dayofyear",
    vmin=vmin,
    vmax=vmax,
    cmap=cmap,
    legend=False,
    style_kwds=dict(color="black"),
)
gf_dt.explore(m=m, column="dayofyear_2", vmin=vmin, vmax=vmax, cmap=cmap)
Make this Notebook Trusted to load map: File -> Trust Notebook
print("Maxar Stereo Acquisition DOY - ICESat-2 DOY")
for i in range(len(gf_i)):
    stereo_pair = gf_i.iloc[[i]]
    gf_dt = stereo_pair.overlay(gf_is2, how="intersection")
    min_dt = (gf_dt.dayofyear_1 - gf_dt.dayofyear_2).min()
    print(stereo_pair.stereo_pair_identifiers.values[0], min_dt)

Maxar Stereo Acquisition DOY - ICESat-2 DOY

['11a11283-6661-4f91-9bcb-5dab3c6a5d02-inv'] -2
['21d1d4b2-47a4-4c4d-a26f-5c688199b518-inv'] -18
['28b60c45-3164-4cd0-9b43-67dac4e8c7b4-inv'] -36
['41fef87b-0d7c-46c2-a598-bcd178dfb9e0-inv'] -30
['6002024a-ca85-4f3c-8c7f-ce849c175a52-inv'] 27

Summary#

  • Any of these Maxar stereo pairs could be worth ordering. But if the goal is to have coincident acquisitions as close as possible in time ‘11a11283-6661-4f91-9bcb-5dab3c6a5d02-inv’ seems like a good bet:

    • There is an overlapping ICESat-2 track 2 days later

    • There is an overlapping GEDI track 3 days earlier

    • 3DEP Lidar was acquired 27 to 36 days earlier

# Can save this dataframe, or just note STAC ids to work with later ('102001008EC5AC00','102001008BE9BB00')
pair_id = "11a11283-6661-4f91-9bcb-5dab3c6a5d02-inv"
full_frame = gf_maxar[gf_maxar.stereo_pair_id == pair_id]
subset = gf_i[gf_i.stereo_pair_id == pair_id]

# Save the original STAC metadata and polygon
full_frame.to_parquet("/tmp/maxar-mono-images.parquet")
# Also save the subsetted region for this stereo pair
subset.to_parquet("/tmp/stereo-subset.parquet")