.ipynb .pdf

repository open issue suggest edit

Binder

Contents

Why would I use UAVSAR for snow?

L-band SAR penetrates through the snowpack. However when it crosses into the snowpack from the air it refracts at an angle, similar to light entering water. This refraction leads to a phase shift relative to an image with no or less snow. Using this difference in phase between two images we can calculate the change in snow height between flights using:

\[ \Delta d = - \frac{\Delta \phi \lambda}{4 \pi} \frac{1}{\cos^{ } \alpha - \sqrt{\epsilon_{s} - \sin^{2} \alpha}} \]

Where \(\Delta\) d is the change in snow height, \(\Delta \phi\) is the phase shift between two SAR images, \(\lambda\) is the radar wavelength, \(\alpha\) is the incidence angle, and \(\epsilon_{s}\) is the dielectric constant of snow which is dependent on the density and liquid water content.

Fig. 5 Conceptual diagram of radar refraction across the air-snow interface.

Download uavsar data from GitHub

import os
# clone tifs to tmp directory
os.chdir('/tmp/')
if not os.path.exists('/tmp/uavsar-tutorial-data-2022'):
    !git clone --quiet https://github.com/snowex-hackweek/uavsar-tutorial-data-2022.git

# list files downloaded
data_dir = '/tmp/uavsar-tutorial-data-2022/'
os.chdir(data_dir)
# ! ls -l

Import Libraries

# Database imports
from snowexsql.db import get_db
from snowexsql.data import PointData, ImageData, LayerData, SiteData
from snowexsql.conversions import query_to_geopandas

# Uavsar_pytools imports
from uavsar_pytools.snow_depth_inversion import depth_from_phase, phase_from_depth

# Other imports
from os.path import join
from datetime import date
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from shapely.geometry import box
import rasterio as rio
import rioxarray as rxa
import contextily as cx
import holoviews as hv
import rioxarray as rxa
from bokeh.plotting import show
import datashader as ds
from datashader.mpl_ext import dsshow
hv.extension('bokeh', logo=False)
%config InlineBackend.figure_format='retina'

Set variables

# Directory of the uavsar tiffs
data_dir = '/tmp/uavsar-tutorial-data-2022/'

# Mesa Lake Snotel Coordinates
snotel_coords = (-108.05, 39.05)

Phase Change between February 1st and 13th UAVSAR Image Pairs

You learned in the first section how to access and download UAVSAR imagery. For this section the data has already been downloaded, converted to GeoTiffs and cropped down to an area of interest that overlaps the main field sites of Grand Mesa. Lets take a look at the coherence and unwrapped phase between these two flights. If you don’t remember what these two represent check out the previous section of this tutorial.

# Create figures and subplots
fig, axes = plt.subplots(2, 1, figsize = (12,8))

# Select colormap for each image type
vis_dic = {'cor': 'Blues', 'unw':'magma'}

# Loop through coherence and unwrapphed phase images
for i, type in enumerate(vis_dic.keys()):
    # select correct axis
    ax = axes[i]
    # open image with rioxarray
    img = rxa.open_rasterio(join(data_dir, f'{type}.tif'))
    # calculate visualization parameters
    vmin, vmax = img.quantile([0.1,0.9])
    # plot images
    img.plot(ax = ax, vmin = vmin, vmax = vmax, cmap = vis_dic[type], zorder = 1, alpha = 0.7)
    # zoom out a bit
    ax.set_xlim(-108.28,-108)
    ax.set_ylim(38.98, 39.08)
    # add topo basemap
    cx.add_basemap(ax, crs=img.rio.crs, alpha = 0.8, source = cx.providers.USGS.USTopo) #cx.providers.USGS.USTopo)
    # turn off labels
    ax.xaxis.label.set_visible(False)
    ax.yaxis.label.set_visible(False)
# set titles
axes[0].set_title('Coherence')
axes[1].set_title('Unwrapped Phase Change')

plt.show()

fig, ax = plt.subplots(figsize = (12,8))

# Plot the snotel location
ax.scatter(x = snotel_coords[0], y = snotel_coords[1], marker = 'x', color = 'black')

# Plot bounding box of uavsar
uavsar_bounds = rxa.open_rasterio(join(data_dir, f'cor.tif')).rio.bounds()
x,y = box(*uavsar_bounds).exterior.xy
ax.plot(x,y, color = 'blue')

# Set overview bounds
ax.set_xlim(-108.4,-107.75)
ax.set_ylim(38.75, 39.3)

# Add background map
# Unfortunately, the below line does not render as of 20240301
# cx.add_basemap(ax, crs=img.rio.crs, source = cx.providers.Stamen.TerrainLabels)
cx.add_basemap(ax, crs=img.rio.crs, alpha = 0.8, source = cx.providers.USGS.USImageryTopo)
plt.title('Overview Map')
plt.show()

Using the SnowEx SQL Database to collect snow depth and lidar datasets

Lets explore how many overlapping depth observations we have between these two days.

# This is what you will use for all of hackweek to access the db
db_name = 'snow:hackweek@db.snowexdata.org/snowex'

# Using the function get_db, we receive 2 ways to interact with the database
engine, session = get_db(db_name)

# Its convenient to store a query like the following 
qry = session.query(PointData)

# Filter to snow depths
qry = qry.filter(PointData.type == 'depth')
qry = qry.filter(PointData.site_name == 'Grand Mesa')
qry = qry.filter(PointData.instrument != 'Mala 800 MHz GPR')

# Then filter on it first date. We are gonna get one day either side of our flight date
qry_feb1 = qry.filter(PointData.date >= date(2020, 1, 31))
qry_feb1 = qry_feb1.filter(PointData.date <= date(2020, 2, 2))
df_feb_1 = query_to_geopandas(qry_feb1, engine)

# Get depths from second flight date
qry_feb12 = qry.filter(PointData.date >= date(2020, 2, 11))
qry_feb12 = qry_feb12.filter(PointData.date <= date(2020, 2, 13))
df_feb_12 = query_to_geopandas(qry_feb12, engine)

# Get depths that were captured on both days
df_both = df_feb_1.overlay(df_feb_12, how = 'intersection')

# Convert crs to match our uavsar images
df_both = df_both.to_crs(epsg = 4326)

# Calculate the snow depth change for each point
df_both['sd_diff'] = df_both.value_2 - df_both.value_1

fig, ax = plt.subplots(figsize = (12,4))

# Plot depth measurements
df_both.plot(ax = ax, column = 'sd_diff', legend = True, legend_kwds = {'label': 'Snow Depth Change [cm]'}, cmap = 'magma')

# Plot the snotel location
snotel_coords = (-108.05, 39.05)
ax.scatter(x = snotel_coords[0], y = snotel_coords[1], marker = 'x', color = 'black')

# Plot bounding box of uavsar
img = rxa.open_rasterio(join(data_dir, f'cor.tif'))
uavsar_bounds = img.rio.bounds()
x,y = box(*uavsar_bounds).exterior.xy
ax.plot(x,y, color = 'blue')

# Set same bounds as uavsar image plot
ax.set_xlim(-108.28,-108)
ax.set_ylim(38.98, 39.08)

# Add background map
cx.add_basemap(ax, crs=img.rio.crs, alpha = 0.8, source = cx.providers.USGS.USImageryTopo)
plt.title('Database Snow Depth Measurements')
plt.show()

Getting the remaining parameters

Incidence Angle

We can recall the formula to calculate snow depth change from incidence angle, phase change, and the snow permittivity.

\[ \Delta d = - \frac{\Delta \phi \lambda}{4 \pi} \frac{1}{\cos^{ } \alpha - \sqrt{\epsilon_{s} - \sin^{2} \alpha}} \]

We have two of these variables already: incidence angle and phase change.

# Create figures and subplots
fig, axes = plt.subplots(2, 1, figsize = (12,8))

# Select colormap for each image type
vis_dic = {'inc': 'Greys', 'unw':'magma'}

# Loop through each image type
for i, type in enumerate(vis_dic.keys()):
    ax = axes[i]
    # Open image with rioxarray
    img = rxa.open_rasterio(join(data_dir, f'{type}.tif'))
    # convert incidence angle from radians to degrees
    if type == 'inc':
        img = np.rad2deg(img)
    # this is a great convenience feature to calculate good visualization levels
    vmin, vmax = img.quantile([0.1,0.9])
    # plot the image
    im = img.plot(ax = ax, vmin = vmin, vmax = vmax, cmap = vis_dic[type], zorder = 1, alpha = 0.7)
    # Zoom out a big
    ax.set_xlim(-108.28,-108)
    ax.set_ylim(38.98, 39.08)
    # Add a topo basemap
    cx.add_basemap(ax, crs=img.rio.crs, alpha = 0.8, source = cx.providers.USGS.USTopo)
    # Remove unnecessary 'x' 'y' labels
    ax.xaxis.label.set_visible(False)
    ax.yaxis.label.set_visible(False)

# Add titles
axes[0].set_title('Incidence Angle')
axes[1].set_title('Unwrapped Phase Change')
plt.show()

Getting Permittivity

We have two ways of getting the \(e_{s}\), or the real part of the snow’s dielectric permittivity. One is by estimating from the snow density. For dry snow we can estimate the permittivity using the density. There are a number of equations for calculating this value, but we will use the equation from Guneriussen et al. 2001:

\[ e_{s} = 1 + 0.0016 \rho + 1.8 1\mathrm{e}{-9} \rho^{3} \]

where \(e_{s}\) is the real part of the snow’s dielectric permittivity and \(\rho\) is the density of the new snow accumulated between the two images in \(\frac{kg}{m^{3}}\).

The other method is to use the directly measured values for permittivity from the field and averaging the top layer.

# Its convenient to store a query like the following 
qry = session.query(LayerData)

# Then filter on it first date. We are gonna get one day either side of second flight date
qry = qry.filter(LayerData.date >= date(2020, 1, 31))
qry = qry.filter(LayerData.date <= date(2020, 2, 2))
qry = qry.filter(LayerData.site_name == 'Grand Mesa')
# Filter to snow density
qry_p = qry.filter(LayerData.type == 'density')
# Change the qry to a geopandas dataframe
df = query_to_geopandas(qry_p, engine)
# create a list to hold the density values
p_values = []
# Loop through each snowpit (each unique site-id is a snowpit) 
for id in np.unique(df.site_id):
    sub = df[df.site_id == id]
    # get the density for the top layer identified in each snowpit
    p = float(sub.sort_values(by = 'depth', ascending = False).iloc[0]['value'])
    # add it our list
    p_values.append(p)
# calculate the mean density of the top layer for each snowpit
mean_new_density = np.nanmean(p_values)
# Use our equation above to estimate our new snow permittivity
es_estimate = 1 + 0.0016*mean_new_density + 1.8e-09*mean_new_density**3

## We can also use snowpits where permittivity was directly observed to compare to
# our density estimates
qry = qry.filter(LayerData.type == 'permittivity')
df = query_to_geopandas(qry, engine)
es_values = []
for id in np.unique(df.site_id):
    sub = df[df.site_id == id]
    es_str = sub.sort_values(by = 'depth', ascending = False).iloc[0]['value']
    if es_str != None:
        es = float(es_str)
        if es != None:
            es_values.append(es)
es_measured = np.nanmean(es_values)

print(f'New snow measured permittivity: {es_measured}. Permittivity from density: {es_estimate}')

New snow measured permittivity: 1.2042714285714287. Permittivity from density: 1.2761086947906

Now we have a new snow permittivity (either from density or directly measured) and we can use that along with our unwrapped phase to calculate the Uavsar snow depth change.

Take a moment to code up the formula for snow depth change from phase and incidence angle:

\[ \Delta d = - \frac{\Delta \phi \lambda}{4 \pi} \frac{1}{\cos^{ } \alpha - \sqrt{\epsilon_{s} - \sin^{2} \alpha}} \]

Where \(\Delta\) d is the change in snow height, \(\Delta \phi\) is the phase shift between two SAR images, \(\lambda\) is the radar wavelength, \(\alpha\) is the incidence angle, and \(\epsilon_{s}\) is the dielectric constant of snow which is dependent on the density and liquid water content.

# Open up our two rasters (unwrapped phase and incidence angle)
unw = rxa.open_rasterio(join(data_dir, f'unw.tif'))
inc = rxa.open_rasterio(join(data_dir, f'inc.tif'))

# This uses the pytool's function to directly give you snow depth change
# feel free to rerun with this to check your results
# https://github.com/SnowEx/uavsar_pytools/blob/main/uavsar_pytools/snow_depth_inversion.py
sd_change = depth_from_phase(unw, inc, density = mean_new_density)

# convert to centimeters from meters
sd_change = sd_change*100

No permittivity data provided -- calculating permittivity from snow density using method guneriussen2001.

# Now we can plot the results!
f, ax = plt.subplots(figsize = (12,8))

# Plot our uavsar snow depth change
sd_change.plot(ax = ax, cmap = 'Blues', vmin = -10, vmax = 10)
# plot black shadow for field observations
df_both.plot(ax = ax, color = 'black', markersize = 90)
# plot field observed snow depth difference
df_both.plot(ax = ax, column = 'sd_diff', legend = True, cmap = 'Blues', vmin = -10, vmax = 10)
# add snotel coordinates
ax.scatter(x = snotel_coords[0], y = snotel_coords[1], marker = 'x', color = 'black')
# turn off labels
ax.xaxis.label.set_visible(False)
ax.yaxis.label.set_visible(False)
# set title
ax.set_title('Uavsar Snow Depth Inversion vs Field Observations')

## Uncomment this to zoom in on the measured results
# ax.set_xlim(-108.14, -108.23)
# ax.set_ylim(39, 39.05)

plt.show()

Numerical Comparison

We can now extract the snow depth change at each measured point and compare them to the pit values of snow depth change.

# We are going to save our results from rioxarray to a tiff
sd_change_fp = join(data_dir,'gm_phase_dz.tif')
sd_change.rio.to_raster(sd_change_fp)
with rio.open(sd_change_fp) as src:
    coord_list = [(x,y) for x,y in zip(df_both['geometry'].x , df_both['geometry'].y)]
    df_both['uavsar_sd'] = [x[0] for x in src.sample(coord_list)]

f, ax = plt.subplots(figsize = (12,8))
df_both['geometry-str'] = df_both['geometry'].astype(str)
df_dis = df_both.dissolve('geometry-str', aggfunc = 'mean')
field_sd_std = df_both.dissolve('geometry-str', aggfunc = 'std')['sd_diff'].values
ax.errorbar(x = df_dis.uavsar_sd, y = df_dis.sd_diff, yerr = field_sd_std, fmt="o")

lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]
def rmse(predictions, targets):
    return np.sqrt(((predictions - targets) ** 2).mean())
rmse_sd = rmse(df_both['sd_diff'], df_both['uavsar_sd'])
print(f'RMSE between uavsar and field observations is {rmse_sd} cm')

# now plot both limits against each other
ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)
ax.set_aspect('equal')
ax.set_xlim(lims)
ax.set_ylim(lims)
ax.set_xlabel('Uavsar Snow Depth Change')
ax.set_ylabel('Field Measured Snow Depth Change')
plt.show()

RMSE between uavsar and field observations is 7.903946685428741 cm

/usr/share/miniconda3/envs/hackweek/lib/python3.10/site-packages/numpy/core/_methods.py:44: RuntimeWarning: invalid value encountered in reduce
  return umr_minimum(a, axis, None, out, keepdims, initial, where)
/usr/share/miniconda3/envs/hackweek/lib/python3.10/site-packages/numpy/core/_methods.py:40: RuntimeWarning: invalid value encountered in reduce
  return umr_maximum(a, axis, None, out, keepdims, initial, where)

Comparison to Lidar

# Create figures and subplots
fig, axes = plt.subplots(3, 1, figsize = (12,8))

lidar = rxa.open_rasterio(join(data_dir, 'sd_lidar.tif'))

diff = lidar.copy()
diff = diff - sd_change

vmin, vmax = sd_change.quantile([0.1,0.9])
sd_change_masked = sd_change.copy()
sd_change_masked.data[np.isnan(lidar).data] = np.nan
sd_change_masked.plot(ax = axes[0], vmin = vmin, vmax = vmax, cmap = 'Blues', zorder = 1, alpha = 0.7)
lidar.plot(ax = axes[1], vmin = vmin, vmax = vmax, cmap = 'Blues', zorder = 1, alpha = 0.7)
diff.plot(ax = axes[2], vmin = vmin, vmax = vmax, cmap = 'Blues', zorder = 1, alpha = 0.7)

for ax in axes:
    ax.xaxis.set_visible(False)
    ax.yaxis.set_visible(False)

axes[0].set_title('Uavsar Snow Depth Change')
axes[1].set_title('Lidar Snow Depth Change')
axes[2].set_title('Snow Depth Difference')
plt.tight_layout()
plt.show()

/usr/share/miniconda3/envs/hackweek/lib/python3.10/site-packages/matplotlib/colors.py:621: RuntimeWarning: overflow encountered in multiply
  xa *= self.N

plt.figure(figsize = (12,8))
diffs = diff.values.ravel()
diffs = diffs[diffs < 100]
diffs = diffs[diffs > -100]
plt.hist(diffs, bins = 100, density = True, label = 'Uavsar sd change')
# plt.axvline(sd_change_masked.mean().values, label = 'Uavsar Mean Snow Depth Change', color = 'green')
lidar_vals = lidar.astype(np.float64).values[~lidar.isnull().values]
lidar_vals = lidar_vals[lidar_vals < 100]
lidar_vals = lidar_vals[lidar_vals > -100]
mean_lidar = np.nanmean(lidar_vals)
plt.axvline(mean_lidar, color = 'red', linewidth = 5, label = 'mean lidar sd change')
# plt.axvline(mean_lidar, label = 'Lidar Mean Snow Depth Change', color = 'red')
rmse = np.sqrt(((diffs) ** 2).mean())
print(f'Lidar mean depth change: {sd_change_masked.mean().values} cm, uavsar mean depth change: {mean_lidar} cm')
print(f'Mean difference: {np.nanmean(diffs)} cm, rmse = {rmse} cm')
plt.legend(loc = 'lower left')
plt.xlabel('Snow Depth Change (cm)')
plt.show()

Lidar mean depth change: -2.154362201690674 cm, uavsar mean depth change: -2.049166001104784 cm
Mean difference: 0.1797051727771759 cm, rmse = 8.220202445983887 cm

Interactive Comparison

Examine some around Grand Mesa. Where is the uavsar doing a good job of identifying snow depth changes, where is it doing worse? Why do you think these areas are different?

Some terrain features to explore:

• steep rolls overs

• different aspects (follow the edge of the mesa to see at E, W, S)

• roads

• lakes

• areas with snow drifting vs snow scouring

lidar = rxa.open_rasterio(join(data_dir, 'sd_lidar.tif'))
lidar.name = 'lidar'

sd_change_masked = sd_change.copy()
sd_change_masked.data[np.isnan(lidar).data] = np.nan
sd_change_masked.name = 'UAVSAR SD change'

cor = rxa.open_rasterio(join(data_dir, f'cor.tif'))
cor.name = 'coherence'

tiles = hv.element.tiles.EsriUSATopo().opts()

n = 5 # decimate lidar by a factor of 5
uavsar = hv.Image(hv.Dataset(sd_change_masked.rio.reproject('EPSG:3857')[0,::1,::1], kdims=['x','y'])).opts(cmap = 'Blues', colorbar=False, xaxis = None, yaxis = None, title = 'Uavsar SD Change', clim = (-10, 10))
lidar_img = hv.Image(hv.Dataset(lidar.rio.reproject('EPSG:3857')[0,::n,::n], kdims=['x','y'])).opts(cmap = 'Blues', colorbar=True, xaxis = None, yaxis = None, title= 'Lidar SD Change', clim = (-10, 10))
cor = hv.Image(hv.Dataset(cor.rio.reproject('EPSG:3857')[0,::n,::n], kdims=['x','y'])).opts(cmap = 'Greys', colorbar=True, xaxis = None, yaxis = None, title= 'Coherence', clim = (0, 1))
uv_tile = tiles  * uavsar
lidar_tile = tiles  * lidar_img
cor_tile = tiles * cor
topo = tiles* hv.Image(hv.Dataset(lidar.rio.reproject('EPSG:3857')[0,::n,::n], kdims=['x','y'])).opts(alpha = 0, title= 'Topo', colorbar = False, xaxis = None, yaxis = None)

hv.Layout([uv_tile, lidar_tile, cor_tile, topo]).opts(width = 1000, height = 900)

previous

UAVSAR Data Products

next

SnowEx Mission and Dataset Descriptions