3DEP search by polygon and clip with PDAL¶

This guide will walk you through the steps to search a STAC Catalog using the STAC API to find 3DEP items within a polygon.

This guide will use Python packages geopandas, leafmap, pdal, and pystac_client.

Important You'll also need GDAL and PDAL installed on your system. If you're using Conda to manage your Python environment, these libraries will be installed automatically. Otherwise, you'll need to install GDAL and PDAL manually before using pip to install the Python extensions.

In [ ]:

# Uncomment and run the following line if working out of Google Colab
# !pip install geopandas
# !pip install leafmap
# !pip install pdal
# !pip install pystac_client

# Uncomment and run the following line if working out of Google Colab # !pip install geopandas # !pip install leafmap # !pip install pdal # !pip install pystac_client

In [ ]:

import json
import sys

import geopandas as gpd
import leafmap
import pdal
from pystac_client import Client

import json import sys import geopandas as gpd import leafmap import pdal from pystac_client import Client

You will need local polygon data in either shapefile or GeoPackage format for the next cell. GeoPandas will read the local file into a DataFrame. Alternatively, you can comment out the first line, uncomment the second line, and provide the URL to a remote polygon dataset, such as a zipped shapefile.

In [ ]:

vector_file = "/Volumes/jj/data/purdue/d2s_workshop_example_data/boundary.shp"
# vector_file = "https://workshop.d2s.org/sample_data/project_boundary.zip"
poly = gpd.read_file(vector_file)

# Now convert the boundary to EPSG:4326
poly = poly.to_crs("EPSG:4326")

vector_file = "/Volumes/jj/data/purdue/d2s_workshop_example_data/boundary.shp" # vector_file = "https://workshop.d2s.org/sample_data/project_boundary.zip" poly = gpd.read_file(vector_file) # Now convert the boundary to EPSG:4326 poly = poly.to_crs("EPSG:4326")

D2S provides a STAC API that hosts a 3DEP collection for you to search. The API is accessible at https://stac-api.d2s.org. For a more user-friendly interface to browse the data exposed by the API, visit https://stac.d2s.org. In the following cells, you'll connect to the STAC API using pystac_client and search the 3DEP collection using your bounding box.

In [ ]:

# Connect to STAC API
client = Client.open("https://stac-api.d2s.org")

# Connect to STAC API client = Client.open("https://stac-api.d2s.org")

In [ ]:

# Get bounding box [xmin, ymin, xmax, ymax] for polygon
bounding_box = poly.total_bounds.tolist()

# Search 3DEP collection
search = client.search(
    max_items=10,
    collections=["3dep"],
    bbox=bounding_box,
)

# Get bounding box [xmin, ymin, xmax, ymax] for polygon bounding_box = poly.total_bounds.tolist() # Search 3DEP collection search = client.search( max_items=10, collections=["3dep"], bbox=bounding_box, )

In [ ]:

# Print STAC Item ID and STAC Browser URL for search results
stac_browser_base_item_url = "https://stac.d2s.org/collections/3dep/items"
items = []
for item in search.items():
    print(f"ID: {item.id}, URL: {stac_browser_base_item_url}/{item.id}")
    # EPSG
    print(f"EPSG: {item.properties['proj:epsg']}")
    # You can also directly access the asset URL from the item
    print(f"URL: {item.assets['ept.json'].href}\n")
    items.append(item)

# Print STAC Item ID and STAC Browser URL for search results stac_browser_base_item_url = "https://stac.d2s.org/collections/3dep/items" items = [] for item in search.items(): print(f"ID: {item.id}, URL: {stac_browser_base_item_url}/{item.id}") # EPSG print(f"EPSG: {item.properties['proj:epsg']}") # You can also directly access the asset URL from the item print(f"URL: {item.assets['ept.json'].href}\n") items.append(item)

Clip and export DTM with PDAL¶

In this final section, you will create and run a PDAL pipeline to clip an EPT dataset from the search results to fit within your polygon boundary. The pipeline will export the clipped point cloud as a .laz file and also generate a DTM for the specified area. For a more in-depth explanation of the pipeline process, refer to the official PDAL documentation: https://pdal.io/en/2.4.3/tutorial/iowa-entwine.html.

In [ ]:

# Choose the item to clip
item = items[1]

# Reproject polygon to match first item's coordinate system
poly = poly.to_crs(f"EPSG:{item.properties['proj:epsg']}")
bounding_box = poly.total_bounds.tolist()

# Get Asset URL for first item
asset_url = item.assets["ept.json"].href
print(bounding_box)

# Choose the item to clip item = items[1] # Reproject polygon to match first item's coordinate system poly = poly.to_crs(f"EPSG:{item.properties['proj:epsg']}") bounding_box = poly.total_bounds.tolist() # Get Asset URL for first item asset_url = item.assets["ept.json"].href print(bounding_box)

In [ ]:

# Provide different filepath and name if desired
out_laz_filename = "./clip.laz"
out_tif_filename = "./clip.tif"

# Pipeline
json_dict = {
    "pipeline": [
        {
            "bounds": f"([{bounding_box[0]}, {bounding_box[2]}], [{bounding_box[1]}, {bounding_box[3]}])",
            "filename": asset_url,
            "type": "readers.ept",
            "tag": "readdata"
        },
        {
            "limits": "Classification![7:7]",
            "type": "filters.range",
            "tag": "nonoise"
        },
        {
            "assignment": "Classification[:]=0",
            "type": "filters.assign",
            "tag": "wipeclasses"
        },
        {
            "out_srs": "EPSG:26916",
            "type": "filters.reprojection",
            "tag": "reprojectUTM"
        },
        {
            "tag": "groundify",
            "type": "filters.smrf"
        },
        {
            "limits": "Classification[2:2]",
            "type": "filters.range",
            "tag": "classify"
        },
        {
            "filename": out_laz_filename,
            "inputs": [ "classify" ],
            "type": "writers.las",
            "tag": "writerslas"
        },
        {
            "filename": out_tif_filename,
            "gdalopts": "tiled=yes,     compress=deflate",
            "inputs": [ "writerslas" ],
            "nodata": -9999,
            "output_type": "idw",
            "resolution": 1,
            "type": "writers.gdal",
            "window_size": 6
        }
    ]
}

# Provide different filepath and name if desired out_laz_filename = "./clip.laz" out_tif_filename = "./clip.tif" # Pipeline json_dict = { "pipeline": [ { "bounds": f"([{bounding_box[0]}, {bounding_box[2]}], [{bounding_box[1]}, {bounding_box[3]}])", "filename": asset_url, "type": "readers.ept", "tag": "readdata" }, { "limits": "Classification![7:7]", "type": "filters.range", "tag": "nonoise" }, { "assignment": "Classification[:]=0", "type": "filters.assign", "tag": "wipeclasses" }, { "out_srs": "EPSG:26916", "type": "filters.reprojection", "tag": "reprojectUTM" }, { "tag": "groundify", "type": "filters.smrf" }, { "limits": "Classification[2:2]", "type": "filters.range", "tag": "classify" }, { "filename": out_laz_filename, "inputs": [ "classify" ], "type": "writers.las", "tag": "writerslas" }, { "filename": out_tif_filename, "gdalopts": "tiled=yes, compress=deflate", "inputs": [ "writerslas" ], "nodata": -9999, "output_type": "idw", "resolution": 1, "type": "writers.gdal", "window_size": 6 } ] }

In [ ]:

# Execute pipeline
pipeline = pdal.Pipeline(json.dumps(json_dict))
count = pipeline.execute()

# Execute pipeline pipeline = pdal.Pipeline(json.dumps(json_dict)) count = pipeline.execute()

In [ ]:

# Display DTM on map
m = leafmap.Map()
m.add_raster("clip.tif", colormap="gray")
m.add_geojson(json.loads(poly.to_json(to_wgs84=True)), style={"color": "yellow", "weight": 3, "fill": False})
m

# Display DTM on map m = leafmap.Map() m.add_raster("clip.tif", colormap="gray") m.add_geojson(json.loads(poly.to_json(to_wgs84=True)), style={"color": "yellow", "weight": 3, "fill": False}) m