GPR and Lidar - Snow Observations Cookbook

Author:

This notebook is designed to obtain snow depth, density, and SWE from ground-penetrating radar data. We will use GPR data gathered during the SnowEx Alaska campaigns, with validation from the airborne lidar obtained during the same campaign.

The text and code in this notebook is adapted from Randall Bonnell’s tutorial on lidar and GPR during the 2024 SnowEx Hackweek, found here: https://snowex-2024.hackweek.io/tutorials/gpr_lidar/GPR_Lidar_HackweekTutorial.html

1. GPR Methods for the Retrieval of Snow Depth and SWE¶

What is GPR?¶

GPR transmits a radar signal into the snowpack, which then reflects off objects/interfaces with contrasting dielectric permittivity. The GPR records the amplitude and two-way travel time (twt) of the reflections.
Dielectric permittivity refers to the dielectric properties of the snowpack that define how EM energy trasmits through the medium.
Usually, we are interested in the snow-ground interface, and we measure the snowpack thickness (depth) in two-way travel time (in nanoseconds).
Most analysis-ready GPR products have twt, snow depth, and SWE variables. Some have also been updated to include derived snow density

For this notebook, we will start from the two-way travel time data to derive our snow properties of interest, and compare it to airborne lidar data. The snowexsql package will be needed.

Note that the results derived here are for the user’s reference - most of the GPR products posted on NSIDC already have snow depth and SWE available as variables.

2. SnowEx23 GPR/Lidar Derived Permittivities/Densities in the Boreal Forest, Alaska¶

Deriving Snow Density at Farmer’s Loop/Creamer’s Field¶

For this first step, we will use:

Airborne lidar data collected on 11 March, 2023.
GPR data collected on 7, 11, 13, and 16 March, 2023.

#1.1 Load relevant packages
import os
import numpy as np 
from datetime import date
from scipy.spatial import cKDTree

#packages for figures
import matplotlib.pyplot as plt
from rasterio.plot import show

#geospatial packages
import geopandas as gpd #for vector data
import xarray as xr
import pandas as pd
from shapely.geometry import box, Point

import rasterio as rio
import rioxarray as rxa

from snowexsql.lambda_client import SnowExLambdaClient

# Initialize client
client = SnowExLambdaClient()

# Get all measurement classes dynamically
classes = client.get_measurement_classes()
LayerMeasurements = classes['LayerMeasurements']
PointMeasurements = classes['PointMeasurements']
RasterMeasurements = classes['RasterMeasurements']

print("🔍 Testing Lambda connection...")
connection_test = client.test_connection()
print(f"✅ Connected: {connection_test.get('connected', False)}")
if connection_test.get('connected'):
    print(f"📊 Database: {connection_test.get('version', 'Unknown version')}")
else:
    print("❌ Connection failed")

🔍 Testing Lambda connection...

✅ Connected: True
📊 Database: PostgreSQL 16.10 on x86_64-conda-linux-gnu, compiled by x86_64-conda-linux-gnu-cc (conda-forge gcc 14.3.0-4) 14.3.0, 64-bit

Part 1: Load the GPR data from the SnowEx data base¶

# Query all point measurements within 10000 m of Creamer's Field center
# Creamer's Field coordinates: 64.872265°N, 147.695426°W (EPSG:4326)
# Converting to a Point geometry for the query

# Center point of Creamer's Field
creamers_field_center = Point(-147.695426, 64.872265)

# Query with 1000m buffer around the center point
df = PointMeasurements.from_area(
    pt=creamers_field_center,
    crs=4326,  # WGS84 coordinate system
    buffer=10000,  # 10000 meters
    date = date(2023, 3, 11),
    type='two_way_travel',
    observer='Randall Bonnell',
    limit=25000,  
    verbose=False
)

# Let's look at the distribution of gpr two-way travel times and estimated snow depths
#Estimate snow depths from twt by assuming a velocity of 0.25 m/ns --> Is this an appropriate velocity estimate?
df['Depth_estimated'] = (df['value']/2)*0.25

ax1 = df.plot.hist(column=["value"], edgecolor='black', title='two-way travel time (ns)')
ax2 = df.plot.hist(column=["Depth_estimated"], edgecolor='black', title='Snow depth (m)')

Load Lidar-Derived Canopy Height and Snow Depth¶

To compare against the GPR data, we will load airborne lidar data products, specifically canopy height and snow depths. To facilitate our analysis, we will create a bounding box for the lidar results using our GPR data.

#Extract x/y limits from GPR data --> these will be used when loading the lidar snow depths
bounds = df.total_bounds

# Create a bounding box
gpr_limits = box(*bounds)

We will load the lidar data directly from an S3 bucket hosted on AWS:

# S3 base URL for tutorial data
S3_BASE_URL = "https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/"

# Load Lidar Canopy Heights

# #Read in the canopy heights raster from Farmer's Loop/Creamer's Field Alaska
flcf_ds = rxa.open_rasterio(
    S3_BASE_URL + "SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.tif",
    mask_and_scale=True
).rio.clip_box(
    minx=gpr_limits.bounds[0],
    miny=gpr_limits.bounds[1],
    maxx=gpr_limits.bounds[2],
    maxy=gpr_limits.bounds[3],
    crs=df.crs  # gpr_limits is in EPSG:4326, this reprojects on the fly
)

# Reproject GPR points to match the lidar CRS
df_reprojected = df.to_crs(flcf_ds.rio.crs)

fig, ax = plt.subplots(figsize=(10, 8))
flcf_ds.squeeze().plot(ax=ax, cmap='Greens', vmin=0, vmax=5, 
                               cbar_kwargs={'label': 'Canopy Height (m)'})
ax.set_title('Canopy Height (m) with GPR Points')
df_reprojected.plot(ax=ax, color='red', markersize=2)
plt.show()

Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.hdr: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.HDR: 403

Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.IMD: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.imd: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.RPB: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.rpb: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.PVL: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.pvl: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_rpc.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_RPC.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_metadata.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_METADATA.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_MTL.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_MTL.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/METADATA.DIM: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/metadata.dim: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_metadata.xml: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0_METADATA.XML: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.RPC: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.rpc: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/summary.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SUMMARY.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDR_Lidar_FLCF_CH_0.5M_20221024_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDR_Lidar_FLCF_CH_0.5M_20221024_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDRX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDRX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPC_Lidar_FLCF_CH_0.5M_20221024_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPC_Lidar_FLCF_CH_0.5M_20221024_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPCX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPCX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.tif.msk: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_CH_0.5M_20221024_V01.0.tif.MSK: 403

# Load Lidar Snow depths 

# #Read in the snow depth raster from Farmer's Loop/Creamer's Field Alaska
flcf_ds = rxa.open_rasterio(
    S3_BASE_URL + "SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.tif",
    mask_and_scale=True
).rio.clip_box(
    minx=gpr_limits.bounds[0],
    miny=gpr_limits.bounds[1],
    maxx=gpr_limits.bounds[2],
    maxy=gpr_limits.bounds[3],
    crs=df.crs  # gpr_limits is in EPSG:4326, this reprojects on the fly
)

# Reproject GPR points to match the lidar CRS
df_reprojected = df.to_crs(flcf_ds.rio.crs)

fig, ax = plt.subplots(figsize=(10, 8))
flcf_ds.squeeze().plot(ax=ax, cmap='Blues', vmin=0, vmax=1.5, 
                               cbar_kwargs={'label': 'Snow Depth (m)'})
ax.set_title('Snow Depth (m) with GPR Points')
df_reprojected.plot(ax=ax, color='red', markersize=2)
plt.show()

Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.hdr: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.HDR: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.IMD: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.imd: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.RPB: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.rpb: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.PVL: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.pvl: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_rpc.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_RPC.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_metadata.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_METADATA.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_MTL.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_MTL.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_metadata.xml: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0_METADATA.XML: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.RPC: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.rpc: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDR_Lidar_FLCF_SD_0.5M_20230311_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDR_Lidar_FLCF_SD_0.5M_20230311_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDRX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/HDRX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPC_Lidar_FLCF_SD_0.5M_20230311_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPC_Lidar_FLCF_SD_0.5M_20230311_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPCX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.txt: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/RPCX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.TXT: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.tif.msk: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/lidar/SNEX23_Lidar_FLCF_SD_0.5M_20230311_V01.0.tif.MSK: 403

Match GPR data with lidar data¶

To match the GPR and lidar data, we can either rasterize the GPR data or vectorize the lidar data. For simplicity, we will vectorize the lidar data and perform a nearest neighbor search.

The GPR data is ~0.1 m resolution, whereas the lidar data is ~0.5 m resolution. So, we can expect ~5 GPR data points per lidar pixel.

# #2.1 Let's learn a bit about the resolution of the lidar rasters

height, width = flcf_ds.squeeze().shape  # Find the height and width of the array

# Get resolution directly from rioxarray
res = flcf_ds.rio.resolution()
print("The x resolution of the snow depth raster is:", abs(res[0]))
print("The y resolution of the snow depth raster is:", abs(res[1]))

# Extract x/y coordinate vectors from the DataArray directly
x_lidar = flcf_ds.x.values
y_lidar = flcf_ds.y.values

# Create meshgrid for vectorizing later
x_lidar_2d, y_lidar_2d = np.meshgrid(x_lidar, y_lidar)

The x resolution of the snow depth raster is: 0.5
The y resolution of the snow depth raster is: 0.5

# # 2.2 Matching GPR to the lidar grid

# #Two conceptual paths forward: rasterize the GPR data, or convert lidar data to points

# Let's vectorize the raster data using the meshgrid created above
x_lidar_vec = x_lidar_2d.flatten()
y_lidar_vec = y_lidar_2d.flatten()
flcf_ds_vec = flcf_ds.squeeze().values.flatten()  # .values converts DataArray to numpy array

# Pull vectors from geo dataframe - reproject GPR to match lidar CRS first
df_reprojected = df.to_crs(flcf_ds.rio.crs)
gpr_arr = np.stack([df_reprojected.geometry.x, df_reprojected.geometry.y, df_reprojected['value']], axis=1)
gpr_x = gpr_arr[:, 0]
gpr_y = gpr_arr[:, 1]
gpr_twt = gpr_arr[:, 2].reshape(len(gpr_arr[:, 2]), 1)

print("Number of lidar points:", len(x_lidar_vec))
print("Number of GPR points:", len(gpr_x))

Number of lidar points: 487925
Number of GPR points: 20213

For the nearest-neighbor approach, we will be using a K-D tree to efficiently find adjacent points.

# 2.3 Create sets of coordinates for the nearest neighbors search
coordinates_set1 = np.column_stack((x_lidar_vec, y_lidar_vec))
coordinates_set2 = np.column_stack((gpr_x, gpr_y))

# Build KDTree from the GPR coordinates
tree = cKDTree(coordinates_set2)

# Define the radius (in meters) - half the lidar pixel size
radius = 0.25

# Function to find the median of travel times within a radius
def find_median_travel_time_within_radius(point, tree, coordinates_set2, gpr_twt, radius):
    indices = tree.query_ball_point(point, radius)
    if indices:
        # Retrieve travel times for the nearest neighbors
        neighbor_twt = gpr_twt[indices]
        median_twt = np.median(neighbor_twt)
        return median_twt
    else:
        return np.nan  # Return NaN if no neighbors are within the radius

# Find medians for each lidar point
medians = np.array([find_median_travel_time_within_radius(point, tree, coordinates_set2, gpr_twt, radius) for point in coordinates_set1])
print("Median twt array size:", medians.shape)
print("Number of non-nan values:", np.sum(~np.isnan(medians)))

Median twt array size: (487925,)
Number of non-nan values: 3788

The GPR data is not as spatially continuous as the lidar data, so twt_median has a lot of NaN values. To reduce memory usage, let’s remove those.

#2.4 Before we get to the math part, let's clear out the nan's from all important vectors:
#Create mask for gpr medians that are nan's
mask = np.isnan(medians)

#Remove entries from the lidar snow depth, x, and y vectors that align with the nan twt values
flcf_ds_vec_clean = flcf_ds_vec[~mask]
coordinates_set1_clean=coordinates_set1[~mask]

#Lastly, remove entries from the twt medians
medians_clean = medians[~mask]

#Let's check the new size of the twt array
print(medians_clean.shape)

(3788,)

Now that the data is matched, and we have removed NaN values, we can now calculate relative permittivity and snow density.

# Calculate relative permittivity

c=0.2998 # The speed of light in a vacuum
e_s = ((c * medians_clean) / (2 * flcf_ds_vec_clean)) ** 2

#And then calculate density
rho_s = ((np.sqrt(e_s) - 1) / 0.845) * 1000

Examine the derived densities¶

Now that we’ve gone through all of this trouble, let’s take a look at the derived snow densities.

# Plot the snow density
plt.figure()
plt.scatter(coordinates_set1_clean[:,0], coordinates_set1_clean[:,1], s=10, c=rho_s, cmap='viridis', clim=(0, 500), edgecolor=None)

plt.colorbar()
plt.title('Snow Density (kg m-3)')
plt.show()

We can also check the histogram distribution of snow density.

# Define bin edges
bin_edges = np.arange(np.min(rho_s), np.max(rho_s), 25)  # Create bin edges from min(x) to max(x) with step size 25

# Create the histogram
plt.figure()  # Create a new figure
plt.hist(rho_s, bins=bin_edges, edgecolor=None)  # Plot histogram with specified bin edges

plt.title('Snow Density Histogram')

# Show the plot
plt.show()

There’s some pretty unrealistic values in the data, so let’s zoom in to more realistic densities.

# Re-define bin edges
bin_edges = np.arange(0, 500, 25)  # Create bin edges from min(x) to max(x) with step size 25

# Create the histogram
plt.figure()  # Create a new figure
plt.hist(rho_s, bins=bin_edges, edgecolor='black')  # Plot histogram with specified bin edges

plt.title('Snow Density Histogram')

# Show the plot
plt.show()

Even at this range, there is a fair amount of random error. This may be caused by the measurement accuracy of the lidar or the GPR, the depth of the snow, and/or geolocation uncertainty.

A user could further improve these densities by upsampling the lidar to match the GPR’s geolocation uncertainty (3 m in this case), remove erroneous values in the snow depth or permittivity data, or run a spatial averaging filter.

References by topic¶

Lidar Datasets¶

Larsen (2024). Larsen (2024)

Relevant GPR LWC Studies¶

Webb et al. (2018). Webb et al. (2018)
Webb et al. (2020). Webb et al. (2020)
Bonnell et al. (2021). Bonnell et al. (2021)
Webb et al. (2022). Webb et al. (2022)

Relevant GPR Density Studies¶

Yildiz et al. (2021). Yildiz et al. (2021)
McGrath et al. (2022). McGrath et al. (2022)
Bonnell et al. (2023). Bonnell et al. (2023)

References¶

Larsen, C. (2024). SnowEx23 Airborne Lidar-Derived 0.25M Snow Depth and Canopy Height, Version 1. NASA National Snow. 10.5067/BV4D8RRU1H7U
Webb, R. W., Jennings, K. S., Fend, M., & Molotch, N. P. (2018). Combining Ground‐Penetrating Radar With Terrestrial LiDAR Scanning to Estimate the Spatial Distribution of Liquid Water Content in Seasonal Snowpacks. Water Resources Research, 54(12). 10.1029/2018wr022680
Webb, R. W., Wigmore, O., Jennings, K., Fend, M., & Molotch, N. P. (2020). Hydrologic connectivity at the hillslope scale through intra‐snowpack flow paths during snowmelt. Hydrological Processes, 34(7), 1616–1629. 10.1002/hyp.13686
Bonnell, R., McGrath, D., Williams, K., Webb, R., Fassnacht, S. R., & Marshall, H.-P. (2021). Spatiotemporal Variations in Liquid Water Content in a Seasonal Snowpack: Implications for Radar Remote Sensing. Remote Sensing, 13(21), 4223. 10.3390/rs13214223
Webb, R. W., Musselman, K. N., Ciafone, S., Hale, K. E., & Molotch, N. P. (2022). Extending the vadose zone: Characterizing the role of snow for liquid water storage and transmission in streamflow generation. Hydrological Processes, 36(3). 10.1002/hyp.14541
Yildiz, S., Akyurek, Z., & Binley, A. (2021). Quantifying snow water equivalent using terrestrial ground penetrating radar and unmanned aerial vehicle photogrammetry. Hydrological Processes, 35(5). 10.1002/hyp.14190
McGrath, D., Bonnell, R., Zeller, L., Olsen-Mikitowicz, A., Bump, E., Webb, R., & Marshall, H.-P. (2022). A Time Series of Snow Density and Snow Water Equivalent Observations Derived From the Integration of GPR and UAV SfM Observations. Frontiers in Remote Sensing, 3. 10.3389/frsen.2022.886747
Bonnell, R., McGrath, D., Hedrick, A. R., Trujillo, E., Meehan, T. G., Williams, K., Marshall, H., Sexstone, G., Fulton, J., Ronayne, M. J., Fassnacht, S. R., Webb, R. W., & Hale, K. E. (2023). Snowpack relative permittivity and density derived from near‐coincident lidar and ground‐penetrating radar. Hydrological Processes, 37(10). 10.1002/hyp.14996