.ipynb .pdf

repository open issue suggest edit

Binder

Contents

Why would I use Uavsar for snow?

L-band radar imagery 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 energy 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 test data from GitHub

import os
os.chdir('/tmp/')
!git clone https://github.com/SnowEx/uavsar-tutorial-data.git

Cloning into 'uavsar-tutorial-data'...

remote: Enumerating objects: 23, done.
remote: Counting objects:   6% (1/15)
remote: Counting objects:  13% (2/15)
remote: Counting objects:  20% (3/15)
remote: Counting objects:  26% (4/15)
remote: Counting objects:  33% (5/15)
remote: Counting objects:  40% (6/15)
remote: Counting objects:  46% (7/15)
remote: Counting objects:  53% (8/15)
remote: Counting objects:  60% (9/15)
remote: Counting objects:  66% (10/15)
remote: Counting objects:  73% (11/15)
remote: Counting objects:  80% (12/15)
remote: Counting objects:  86% (13/15)
remote: Counting objects:  93% (14/15)
remote: Counting objects: 100% (15/15)
remote: Counting objects: 100% (15/15), done.
remote: Compressing objects:   8% (1/12)
remote: Compressing objects:  16% (2/12)
remote: Compressing objects:  25% (3/12)

remote: Compressing objects:  33% (4/12)

remote: Compressing objects:  41% (5/12)

remote: Compressing objects:  50% (6/12)

remote: Compressing objects:  58% (7/12)
remote: Compressing objects:  66% (8/12)
remote: Compressing objects:  75% (9/12)
remote: Compressing objects:  83% (10/12)
remote: Compressing objects:  91% (11/12)
remote: Compressing objects: 100% (12/12)

remote: Compressing objects: 100% (12/12), done.
Receiving objects:   4% (1/23)
Receiving objects:   8% (2/23)

Receiving objects:  13% (3/23)

Receiving objects:  17% (4/23), 23.80 MiB | 47.59 MiB/s

Receiving objects:  17% (4/23), 49.10 MiB | 49.09 MiB/s

Receiving objects:  21% (5/23), 49.10 MiB | 49.09 MiB/s

Receiving objects:  26% (6/23), 74.71 MiB | 49.80 MiB/s

Receiving objects:  26% (6/23), 99.75 MiB | 49.87 MiB/s

Receiving objects:  30% (7/23), 99.75 MiB | 49.87 MiB/s

Receiving objects:  34% (8/23), 125.89 MiB | 50.35 MiB/s
Receiving objects:  39% (9/23), 125.89 MiB | 50.35 MiB/s
Receiving objects:  43% (10/23), 125.89 MiB | 50.35 MiB/s
Receiving objects:  47% (11/23), 125.89 MiB | 50.35 MiB/s
Receiving objects:  52% (12/23), 125.89 MiB | 50.35 MiB/s
Receiving objects:  56% (13/23), 125.89 MiB | 50.35 MiB/s

Receiving objects:  60% (14/23), 125.89 MiB | 50.35 MiB/s
Receiving objects:  65% (15/23), 125.89 MiB | 50.35 MiB/s

Receiving objects:  65% (15/23), 147.05 MiB | 49.02 MiB/s

Receiving objects:  69% (16/23), 147.05 MiB | 49.02 MiB/s
Receiving objects:  73% (17/23), 147.05 MiB | 49.02 MiB/s

Receiving objects:  78% (18/23), 147.05 MiB | 49.02 MiB/s

Receiving objects:  82% (19/23), 167.08 MiB | 47.74 MiB/s
Receiving objects:  86% (20/23), 167.08 MiB | 47.74 MiB/s
Receiving objects:  91% (21/23), 167.08 MiB | 47.74 MiB/s

Receiving objects:  91% (21/23), 190.45 MiB | 47.61 MiB/s
Receiving objects:  95% (22/23), 190.45 MiB | 47.61 MiB/s

remote: Total 23 (delta 1), reused 13 (delta 0), pack-reused 8
Receiving objects: 100% (23/23), 190.45 MiB | 47.61 MiB/s
Receiving objects: 100% (23/23), 211.71 MiB | 47.98 MiB/s, done.
Resolving deltas:   0% (0/1)
Resolving deltas: 100% (1/1)
Resolving deltas: 100% (1/1), done.

Updating files:  66% (10/15)

Updating files:  73% (11/15)

Updating files:  80% (12/15)
Updating files:  86% (13/15)

Updating files:  93% (14/15)

Updating files: 100% (15/15)
Updating files: 100% (15/15), done.

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
%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. It is stored in the gm_mesa subdirectory of this tutorial. 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)
#cx.add_basemap(ax, crs=img.rio.crs, alpha = 0.8, source = cx.providers.Stamen.TerrainLabels)
plt.title('Overview Map')
plt.show()

Using Database to collection 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')
# 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]'})
# 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'}

for i, type in enumerate(vis_dic.keys()):
    ax = axes[i]
    img = rxa.open_rasterio(join(data_dir, f'{type}.tif'))
    if type == 'inc':
        img = np.rad2deg(img)
    vmin, vmax = img.quantile([0.1,0.9])
    im = 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)
    ax.xaxis.label.set_visible(False)
    ax.yaxis.label.set_visible(False)
axes[0].set_title('Incidence Angle')
axes[1].set_title('Unwrapped Phase Change')
plt.show()

Setting the Zero Phase Point

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

df = pd.read_csv(join(data_dir, 'mesa_snotel.csv'), skiprows = 3, index_col=['Date'], parse_dates=['Date'])
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]
with rio.open(join(data_dir, 'inc.tif')) as src:
    for val in src.sample([snotel_coords]): 
        snotel_inc = val[0]
snotel_sd_phase_from_sd_change = phase_from_depth(snotel_sd_delta, snotel_inc, density = 175)
snotel_coords = (-108.05, 39.05)
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(f'Snotel snow depth change: {snotel_sd_delta} cm. Phase should be {snotel_sd_phase_from_sd_change} and is currently {snotel_phase}')
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)[https://ieeexplore.ieee.org/document/957273]:

\[ 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)
# Filter to snow depths
qry = qry.filter(LayerData.type == 'permittivity')
# Then filter on it first date. We are gonna get one day either side of our flight date
qry = qry.filter(LayerData.date >= date(2020, 1, 31))
qry = qry.filter(LayerData.date <= date(2020, 2, 2))
qry = qry.limit(100)
df = query_to_geopandas(qry, engine)
es_values = []
for id in np.unique(df.site_id):
    sub = df[df.site_id == id]
    es = float(sub.sort_values(by = 'depth', ascending = False).iloc[0]['value'])
    es_values.append(es)
mean_new_snow_es = np.nanmean(es_values)

# 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))
# Filter to snow density
qry_p = qry.filter(LayerData.type == 'density')
df = query_to_geopandas(qry_p, engine)
p_values = []
for id in np.unique(df.site_id):
    sub = df[df.site_id == id]
    p = float(sub.sort_values(by = 'depth', ascending = False).iloc[0]['value'])
    p_values.append(p)
mean_new_density = np.nanmean(p_values)
density_new_permittivity = 1 + 0.0016*mean_new_density + 1.8e-09*mean_new_density**3

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 = float(sub.sort_values(by = 'depth', ascending = False).iloc[0]['value'])
    es_values.append(es)
mean_new_permittivity = np.nanmean(es_values)

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

New snow measured permittivity: 1.2133225806451613. Permittivity from density: 1.284756617495022

unw = rxa.open_rasterio(join(data_dir, f'unw.tif'))
inc = rxa.open_rasterio(join(data_dir, f'inc.tif'))
# sd_change = depth_from_phase(unw, inc, permittivity = mean_new_snow_es)
sd_change = depth_from_phase(unw, inc, density = mean_new_density)
sd_change = sd_change*100

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

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')
plt.show()

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')
ax.set_xlim(-108.14, -108.23)
ax.set_ylim(39, 39.05)
plt.show()

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)

(-15.517686891555787, 16.522941550353977)

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

diffs = diff.values.ravel()
diffs = diffs[diffs < 100]
diffs = diffs[diffs > -100]
plt.hist(diffs, bins = 100, density = True)
# 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, 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()
plt.show()

Lidar mean depth change: -2.1021041870117188 cm, uavsar mean depth change: -2.049166001104784 cm
Mean difference: 0.1270764023065567 cm, rmse = 8.079691886901855 cm

previous

Uavsar Data Products

next

SnowEx Mission and Dataset Descriptions