Sprint 2, Module A

Summary

We selected an ALS (Airborne Laser Scanning) scan from Shaver Lake and generated its corresponding Canopy Height Model (CHM) and treelist. ALS is a remote sensing technique that uses LiDAR (Light Detection and Ranging) sensors mounted on aircraft to collect high-resolution, three-dimensional data about the Earth’s surface and vegetation. As the aircraft flies over a forest, the LiDAR sensor emits laser pulses toward the ground and measures how long it takes for the light to bounce back. This produces a dense cloud of elevation points that represent both the terrain and the vertical structure of vegetation, including individual trees. Using this data, we constructed a CHM to visualize tree heights and then derived a treelist identifying the location and height of individual trees. We compared these metrics to those from the Independence Lake site to identify structural differences in forest canopy and tree distributions across regions.

Part 1: Creating a CHM

# Task 1 - Your code here
# import packages
import pc_rasterize as pcr
import glob
import numpy as np
import requests
import pandas as pd
from dask.distributed import Client, LocalCluster, Lock
import rasterio
import matplotlib.pyplot as plt
from skimage.feature import peak_local_max
from pyproj import Transformer

# download ALS data for Shaver Lake
url = "https://wifire-data.sdsc.edu/nc/s/nnSMqWfZAN6Cz6m/download?path=%2Fnew_data%2FRaw%20ALS%20and%20Detected%20Trees%2FALS%20data&files=USGS_LPC_CA_SouthernSierra_2020_B20_w0298n4108.laz&downloadStartSecret=aik5h9wexsu"
local_als = "./USGS_LPC_CA_SouthernSierra_2020_B20_w0298n4108_Dream-Team!.laz"

response = requests.get(url, stream=True)
with open(local_als, "wb") as file:
    for chunk in response.iter_content(chunk_size=8192):
        file.write(chunk)

print(f"Downloaded: {local_als}")

Downloaded: ./USGS_LPC_CA_SouthernSierra_2020_B20_w0298n4108_Dream-Team!.laz

pcr.get_file_quickinfo('./USGS_LPC_CA_SouthernSierra_2020_B20_w0298n4108_Dream-Team!.laz')

LidarInfo(n_points=38817490, dims=('X', 'Y', 'Z', 'Intensity', 'ReturnNumber', 'NumberOfReturns', 'ScanDirectionFlag', 'EdgeOfFlightLine', 'Classification', 'Synthetic', 'KeyPoint', 'Withheld', 'Overlap', 'ScanAngleRank', 'UserData', 'PointSourceId', 'GpsTime', 'ScanChannel', 'Amplitude', 'Reflectance', 'Deviation'), crs=CRS.from_wkt('COMPD_CS["NAD83(2011) / UTM zone 11N + NAVD88 height - Geoid18 (m)",PROJCS["NAD83(2011) / UTM zone 11N",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spatial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","6318"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["meter",1],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","6340"]],VERT_CS["NAVD88 height - Geoid18 (m)",VERT_DATUM["North American Vertical Datum 1988",2005],UNIT["meter",1],AXIS["Up",UP],AUTHORITY["EPSG","5703"]]]'), vert_units='meter', bounds_2d=(298000, 4108000, 298999.99, 4108999.99), bounds_3d=(298000, 4108000, 1661.46, 298999.99, 4108999.99, 1832.1), bbox=<POLYGON ((298999.99 4108000, 298999.99 4108999.99, 298000 4108999.99, 29800...>)

# create pipeline for new file 
als = glob.glob("./USGS_LPC_CA_SouthernSierra_2020_B20_w0298n4108_Dream-Team!.laz") # We load our file

# Create a GeoBox grid specification with a 100m buffer around data
geobox = pcr.build_geobox(als, resolution=1., crs="5070")

# Build a lazy CHM raster
chm = pcr.rasterize(
    als,
    geobox,
    cell_func="max",
    # Set custom dask chunk-size
    chunksize=(1000, 1000),
    nodata=np.nan,
    pdal_filters=[
    {
        "type":"filters.range",
        "limits":"Classification[1:5]" # Keep points classified as 'vegetation' - https://desktop.arcgis.com/en/arcmap/latest/manage-data/las-dataset/lidar-point-classification.htm
    },
    {
      "type":"filters.outlier" # Remove outliers
    },
    {
        "type":"filters.hag_nn" # Compute Height Above Ground (HAG) using the nearest-neighbor method. 
    },
    {
        "type":"filters.ferry",
        "dimensions":"HeightAboveGround=>Z" # Move the computed HAG values into the `Z` coordinate for easier processing.  
    },
    {
      "type":"filters.outlier" # Remove outliers again
    },
    {
        "type": "filters.expression",
        "expression": "Z < 75" # Filters out points higher than 75m (unrealistic height in our context) 
    }
    ]
)
# A quick look to our xarray
chm

<xarray.DataArray 'rasterized_pc' (band: 1, y: 1203, x: 1189)> Size: 6MB
dask.array<broadcast_to, shape=(1, 1203, 1189), dtype=float32, chunksize=(1, 1000, 1000), chunktype=numpy.ndarray>
Coordinates:
  * band         (band) int64 8B 1
  * y            (y) float64 10kB 1.812e+06 1.812e+06 ... 1.811e+06 1.811e+06
  * x            (x) float64 10kB -2.029e+06 -2.029e+06 ... -2.028e+06
    spatial_ref  int64 8B 0
Attributes:
    _FillValue:  nan

xarray.DataArray

'rasterized_pc'

• band: 1
• y: 1203
• x: 1189

• dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>

  Array Chunk Bytes 5.46 MiB 3.81 MiB Shape (1, 1203, 1189) (1, 1000, 1000) Dask graph 4 chunks in 5 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • y
    (y)
    float64
    1.812e+06 1.812e+06 ... 1.811e+06

    array([1812214.5, 1812213.5, 1812212.5, ..., 1811014.5, 1811013.5, 1811012.5])

  • x
    (x)
    float64
    -2.029e+06 ... -2.028e+06

    array([-2028689.5, -2028688.5, -2028687.5, ..., -2027503.5, -2027502.5,
           -2027501.5])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["NAD83 / Conus Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5070"]]

    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 / Conus Albers

    grid_mapping_name :

    albers_conical_equal_area

    standard_parallel :

    (29.5, 45.5)

    latitude_of_projection_origin :

    23.0

    longitude_of_central_meridian :

    -96.0

    false_easting :

    0.0

    false_northing :

    0.0

    spatial_ref :

    PROJCS["NAD83 / Conus Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5070"]]

    array(0)

• Indexes: (3)
  • y
    PandasIndex

    PandasIndex(Index([1812214.5, 1812213.5, 1812212.5, 1812211.5, 1812210.5, 1812209.5,
           1812208.5, 1812207.5, 1812206.5, 1812205.5,
           ...
           1811021.5, 1811020.5, 1811019.5, 1811018.5, 1811017.5, 1811016.5,
           1811015.5, 1811014.5, 1811013.5, 1811012.5],
          dtype='float64', name='y', length=1203))

  • x
    PandasIndex

    PandasIndex(Index([-2028689.5, -2028688.5, -2028687.5, -2028686.5, -2028685.5, -2028684.5,
           -2028683.5, -2028682.5, -2028681.5, -2028680.5,
           ...
           -2027510.5, -2027509.5, -2027508.5, -2027507.5, -2027506.5, -2027505.5,
           -2027504.5, -2027503.5, -2027502.5, -2027501.5],
          dtype='float64', name='x', length=1189))

  • band
    PandasIndex

    PandasIndex(Index([1], dtype='int64', name='band'))

• Attributes: (1)

  _FillValue :

  nan

# perform CHM computation
with LocalCluster(n_workers=8, threads_per_worker=1, memory_limit='24GB') as cluster, Client(cluster) as client:
    chm_squeezed = chm.squeeze()  # We remove an extra dimension with no data
    chm_squeezed.rio.to_raster("./chm_Dream-Team!.tiff", tiled=True, lock=Lock("rio"))
print("Successfuly computed the CHM!")

Successfuly computed the CHM!

Part 2: Visualizing the CHM & Creating a Treelist

# Task 2 - Your code here
chm_path = "./chm_Dream-Team!.tiff"
with rasterio.open(chm_path) as src:
    chm = src.read(1)  
    transform = src.transform  # Affine transform for georeferencing
    crs = src.crs  # Get raster Coordinate Reference System

plt.figure(figsize=(10, 6))
plt.imshow(chm, cmap="viridis", interpolation="nearest")
plt.colorbar(label="Canopy Height (m)")
plt.title("Canopy Height Model @ Shaver Lake")
plt.show()

# Normalize CHM 
chm_norm = (chm - np.nanmin(chm)) / (np.nanmax(chm) - np.nanmin(chm))

# Locate tree tops using local maxima filtering
tree_tops = peak_local_max(chm_norm, min_distance=3, threshold_abs=0.1, exclude_border=False)

# Convert pixel coordinates to geospatial coordinates and extract heights
tree_points = []
for row, col in tree_tops:
    x, y = transform * (int(col), int(row))  
    height = chm[row, col] 
    tree_points.append((height, y, x)) 

df = pd.DataFrame(tree_points, columns=["HT", "y", "x"])

# Let´s look at the df
df.head()

	HT	y	x
0	63.060001	1811859.0	-2027711.0
1	62.360001	1811791.0	-2027928.0
2	61.439999	1811907.0	-2028066.0
3	61.439999	1811864.0	-2027709.0
4	61.380001	1811788.0	-2027934.0

# Let's take a look at the distribution of tree heights
plt.figure(figsize=(10, 5))
plt.hist(df["HT"], bins=30, edgecolor="black", alpha=0.7)
plt.xlabel("Tree Height (HT)")
plt.ylabel("Frequency")
plt.title("Distribution of Tree Heights (HT) at Shaver Lake")
plt.grid(axis="y", linestyle="--", alpha=0.7)
plt.show()

# Convert to Lat/Lon system
if crs and crs.to_epsg() != 4326:
    transformer = Transformer.from_crs(crs.to_epsg(), "EPSG:4326", always_xy=True)
    df["x"], df["y"] = transformer.transform(df["x"].values, df["y"].values) 

# Save CSV
df.to_csv("./ttops_Dream-Team!.csv", index=False)
print(f"TreeTops saved with {len(df)} points")

TreeTops saved with 6002 points

Analysis

The canopy height models and tree height distributions at Independence Lake and Shaver Lake reveal several similarities and differences. Both sites exhibit a right-skewed distribution of tree heights, with the majority of trees falling in the mid-height range and fewer tall trees, as seen in their respective histograms. Visually, the canopy height models show clustered vegetation with distinct patches of taller trees and some areas of low or no canopy. However, the Independence Lake graph has a higher maximum tree height (over 70 meters) compared to Shaver Lake (just above 60 meters), and its distribution is more centered around 30 meters, suggesting a generally taller and more uniform canopy. In contrast, Shaver Lake has a broader spread of lower tree heights, with a higher frequency of trees between 15 and 25 meters, indicating a potentially younger or more disturbed forest. Additionally, Independence Lake appears to have denser regions of tall canopy cover, while Shaver Lake shows more scattered patterns with greater presence of low-canopy or open areas across the CHM map. Overall, Independence Lake seems to host a more mature forest structure, while Shaver Lake reflects a younger or more variable landscape.