.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. 6 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'):
    !git clone https://github.com/SnowEx/uavsar-tutorial-data.git

# list files downloaded
data_dir = '/tmp/uavsar-tutorial-data/'
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/'

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

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'}

for i, type in enumerate(vis_dic.keys()):
    ax = axes[i]
    img = rxa.open_rasterio(join(data_dir, f'{type}.tif'))
    vmin, vmax = img.quantile([0.1,0.9])
    img.plot(ax = ax, vmin = vmin, vmax = vmax, cmap = vis_dic[type], zorder = 1, alpha = 0.7)
    ax.set_xlim(-108.28,-108)
    ax.set_ylim(38.98, 39.08)
    cx.add_basemap(ax, crs=img.rio.crs, alpha = 0.8, source = cx.providers.USGS.USTopo) #cx.providers.USGS.USTopo)
    ax.xaxis.label.set_visible(False)
    ax.yaxis.label.set_visible(False)

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
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

# 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()

Snow Depth Change Inversion from Phase

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()

Setting the Zero Phase Change point

Unwrapped phase has one unknown - the zero phase point. This means we have an unknown absolute scene wide shift we can control. We will use the snotel depth change between the two flights to set this unknown and get our absolute phase change.

# read in snotel data. I got this from the USDA site:
# https://wcc.sc.egov.usda.gov/nwcc/site?sitenum=622&state=co
df = pd.read_csv(join(data_dir, 'mesa_snotel.csv'), skiprows = 3, index_col=['Date'], parse_dates=['Date'])
# Calculate the snow depth change at the snotel between the two uavsar flights
snotel_sd_delta = (df[df.index == '2020-02-01']['SNWD.I-1 (in) ']*0.0254).values[0] - (df[df.index == '2020-02-12']['SNWD.I-1 (in) ']*0.0254).values[0]
# Open the incidence angle and find the incidence angle at the snotel
with rio.open(join(data_dir, 'inc.tif')) as src:
    for val in src.sample([snotel_coords]): 
        snotel_inc = val[0]
# This will function from Uavsar_pytools calculates a phase change from a known
# depth change and density of new snow (175 kg/m3)
snotel_sd_phase_from_sd_change = phase_from_depth(snotel_sd_delta, snotel_inc, density = 175)
# These are the coordinates of the snotel for plotting. Again these come from the USDA
snotel_coords = (-108.05, 39.05)
# Open the unwrapped phase and get the pahse
with rio.open(join(data_dir, 'unw.tif')) as src:
    for val in src.sample([snotel_coords]): 
        snotel_phase = val[0]
unw = rxa.open_rasterio(join(data_dir, 'unw.tif'))
# Print out what we have found
print(f'Snotel snow depth change: {snotel_sd_delta} cm. Phase should be {snotel_sd_phase_from_sd_change} and is currently {snotel_phase}')

# This subtracts the phase from our known snow depth from the unwrapped phase
# to shift our unwrapped phase to the right level
unw = unw - (snotel_phase - snotel_sd_phase_from_sd_change)

No permittivity data provided -- calculating permittivity from snow density using method guneriussen2001.
Snotel snow depth change: 0.0 cm. Phase should be 0.0 and is currently -0.9738530516624451

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 zero-point corrected 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'))
# set our wavelength for L-band Uavsar and the value of pi
wavelength = 0.238403545 # meters
pi = np.pi

###############  CODE TO CALCULATE SNOW DEPTH CHANGE  ####################
# take a moment to code this up yourself using the formula above.
# note that the resulting variable should be called sd_change, the unwrapped phase
# is called unw and the incidence angle is inc (from directly above), the final
# piece is the dielectic constant of snow which you can use either of what we 
# calculated right above here es_measured or es_estimate








###############  CODE TO CALCULATE SNOW DEPTH CHANGE  ####################
# 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))

sd_change.plot(ax = ax, cmap = 'Blues', vmin = -10, vmax = 10)

df_both.plot(ax = ax, color = 'black', markersize = 90)
df_both.plot(ax = ax, column = 'sd_diff', legend = True, cmap = 'Blues', vmin = -10, vmax = 10)
ax.scatter(x = snotel_coords[0], y = snotel_coords[1], marker = 'x', color = 'black')
ax.xaxis.label.set_visible(False)
ax.yaxis.label.set_visible(False)
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
]

# 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)

/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)

(-15.922229385375976, 19.55)

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?

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'

tiles = hv.element.tiles.EsriUSATopo().opts()
n = 5 # decimate lidar by a factor of 5
uavsar = hv.Image(hv.Dataset(sd_change_masked[0,::1,::1], kdims=['x','y'])).opts(cmap = 'Blues', colorbar=False, xaxis = None, yaxis = None, title = 'Uavsar SD Change', clim = (-10, 10))
lidar = hv.Image(hv.Dataset(lidar[0,::n,::n], kdims=['x','y'])).opts(cmap = 'Blues', colorbar=True, xaxis = None, yaxis = None, title= 'Lidar SD Change', clim = (-10, 10))
uv_tile = tiles  * uavsar
unw_tile = tiles  * lidar

hv.Layout([uavsar, lidar]).opts(width = 1000, height = 900)

previous

UAVSAR Data Products

next

SnowEx Mission and Dataset Descriptions