Individual glacier surface velocity analysis

Individual glacier surface velocity analysis#

This notebook will build upon the data access and inspection steps in the earlier notebooks and demonstrate basic data analysis and visualization of surface velocity data at the scale of an individual glacier using xarray.

Learning goals#

Concepts#

Visualizing statistical distributions
Vectorized calculation of summary statistics over large dimensions
Working with velocity component vectors
Calculating magnitude of displacement from velocity component vectors

Techniques#

Using xr.reduce() to apply scipy statistical functions designed to ingest numpy arrays to Xarray objects
Re-organize Xarray objects using xr.sortby()
Temporal resampling using xr.resample()
Vectorized computation and reductions using xr.groupby() and xr.map()
Visualization using FacetGrid objects

Other useful resources#

These are resources that contain additional examples and discussion of the content in this notebook and more.

Xarray user guide
How do I… this is very helpful!
Xarray High-level computational patterns discussion of concepts and associated code examples
Parallel computing with dask Xarray tutorial demonstrating wrapping of dask arrays

Software + Setup#

%xmode minimal

Exception reporting mode: Minimal

import os
import json
import urllib.request
import numpy as np
import xarray as xr
import rioxarray as rxr
import geopandas as gpd
import pandas as pd
import scipy

import matplotlib.pyplot as plt

from shapely.geometry import Polygon
from shapely.geometry import Point


import flox

%config InlineBackend.figure_format='retina'

Visualize and examine distributions of different variables#

Let’s first take a look at the velocity components and magnitude of velocity. Plotting a histogram of vx, vy, and v we see that they both have Gaussian distributions. In contrast, magnitude of velocity is positively skewed. The distribution of magnitude of velocity follows a Rician distribution. Because magnitude of velocity is the norm of the two normally-distributed variables (vx and vy) it follows a Rician distribution. To avoid propagating noise in the dataset, we will work primarily with velocity component variables and calculate the magnitude of velocity after reductions have been performed.

To construct these plots, we use a combination of xarray plotting functionality and matplotlib object-oriented plotting

fig,axs=plt.subplots(ncols=3, figsize=(20,5))
hist_y = sample_glacier_raster.vy.plot.hist(ax=axs[0], bins=100)
cumulative_y = np.cumsum(hist_y[0])
axs[0].plot(hist_y[1][1:], cumulative_y, color='orange', linestyle='-', alpha=0.5)

hist_x = sample_glacier_raster.vx.plot.hist(ax=axs[1], bins=100)
cumulative_x = np.cumsum(hist_x[0])
axs[1].plot(hist_x[1][1:], cumulative_x, color='orange', linestyle='-', alpha=0.5)

hist_v = sample_glacier_raster.v.plot.hist(ax=axs[2], bins=100)
cumulative_v = np.cumsum(hist_v[0])
axs[2].plot(hist_v[1][1:], cumulative_v, color='orange', linestyle='-', alpha=0.5)


axs[0].set_title('VY')
axs[1].set_title('VX')
axs[2].set_title('V (magnitude)')

#axs[0].set_xlim(0,50)
#axs[1].set_xlim(0,400)
axs[0].set_ylabel('# obs')
axs[1].set_ylabel('# obs')
axs[2].set_ylabel('# obs')

axs[0].set_xlabel('magnitude of velocity (m/y)')
axs[1].set_xlabel('magnitude of velocity (m/y)')
axs[2].set_xlabel('magnitude of velocity (m/y)')
fig.suptitle('Histogram (blue) and cumulative distribution function (orange) of velocity components and magnitude');

_images/bcd2860459766ca554dabbcbb63b9f23f4f78a28980478af805e36b50a9b3b47.png

We can also look at the skew of each variable. TO do this, we use xr.reduce() and scipy.stats.skew()

print('Skew of vy: ',sample_glacier_raster.vy.reduce(func=scipy.stats.skew, nan_policy='omit', dim=['x','y','mid_date']).data)
print('Skew of vx: ',sample_glacier_raster.vx.reduce(func=scipy.stats.skew, nan_policy='omit', dim=['x','y','mid_date']).data)
print('Skew of v: ',sample_glacier_raster.v.reduce(func=scipy.stats.skew, nan_policy='omit', dim=['x','y','mid_date']).data)

Skew of vy:  -0.4540084857562447
Skew of vx:  0.19356243875503296
Skew of v:  5.005839671572479

Calculating magnitude of velocity#

We’ll first define a function for calculating magnitude of velocity in two ways. Because we want to calculate the magnitude of the displacement vector after we have already reduced the data along a dimension, we write a function that creates two magnitude of velocity variables, one where magnitude is calculated from the means of the vx and vy vectors in space and one where the magnitude is calculated from the medians of the vx and vy vectors in time. We just need to be careful which variable we use.

def calc_velocity_magnitude(ds):
    #naming convention is the dim that variable still has 
    # use mean for reductions of components because components normally distributed
    ds = ds.copy()
    
    ds['v_mag_time'] = np.sqrt(ds.vx.mean(dim=['x','y'])**2 + ds.vy.mean(dim=['x','y'])**2)

    ds['v_mag_space'] = np.sqrt(ds.vx.mean(dim=['mid_date'])**2 + (ds.vy.mean(dim=['mid_date'])**2))
    return ds

sample_glacier_raster_mag = calc_velocity_magnitude(sample_glacier_raster)

Visualize velocity variability#

Now that we have calculated magnitude of velocity, we will plot the histogram of calculated magnitude of velocity values in space (left) and time (right). Note that the scales of the x-axes of the two plots are different, as are the number of bins used in each. These results indicate that the variability of velocity in time is much greater than the variability of velocity in space.

fig,axs=plt.subplots(ncols=2, figsize=(20,5))
sample_glacier_raster_mag.v_mag_space.plot.hist(ax=axs[0], bins=100)
sample_glacier_raster_mag.v_mag_time.plot.hist(ax=axs[1], bins=500)
axs[0].set_title('space')
axs[1].set_title('time')

#axs[0].set_xlim(0,60)
#axs[1].set_xlim(0,400)

axs[0].set_ylabel('# obs')
axs[1].set_ylabel('# obs')
axs[0].set_xlabel('magnitude of velocity (m/y)')
axs[1].set_xlabel('magnitude of velocity (m/y)')

fig.suptitle('Histogram of magnitude of velocity') ;

_images/b287675fb2bca267aa7a45fb40a7b2106debd29d8fddfc2e8698ac45f841c63c.png

We can look at skewness of the magnitude distributions to better understand how these variables behave in space and time.

print('Skew of vmag reduced along spatial dimensions: ',sample_glacier_raster_mag.v_mag_time.reduce(func=scipy.stats.skew, nan_policy='omit', dim=['mid_date']).data)
print('Skew of vmag reduced along temporal dimension: ',sample_glacier_raster_mag.v_mag_space.reduce(func=scipy.stats.skew, nan_policy='omit', dim=['x','y',]).data)

Skew of vmag reduced along spatial dimensions:  7.368324621443265
Skew of vmag reduced along temporal dimension:  2.0117210137206567

Both v_mag_time (reduced along spatial dimensions, exists along temporal dimension) and v_mag_space (reduced along temporal dimension, exists along spatial dimensions) are positively skewed. We also see that the skewness v_mag_time is much greater than v_mag_space, suggesting that the distribution of magnitude of velocity is much more positively skewed (characterized by large, positive outliers) in time than in space. Physically speaking, magnitude of velocity variability is characterized by more extreme outliers over time than across the surface of the glacier.

Checking coverage#

Now that we’ve reduced the dataset, we can look at the coverage of the magnitude variables using xarray methods. First, we want to know how many observations (not NaNs) exist along the time dimension of v_mag_time. We can use xr.DataArray.count():

sample_glacier_raster_mag.v_mag_time.count(dim='mid_date').data

array(927)

We can verify that by using isnull(), notnull() and sum():

sample_glacier_raster_mag.v_mag_time.isnull().sum().data

array(3047)

sample_glacier_raster_mag.v_mag_time.notnull().sum().data

array(927)

Great, now look at space:

sample_glacier_raster_mag.v_mag_space.count(dim=['x','y']).data

array(707)

And checking against the sum() methods:

sample_glacier_raster_mag.v_mag_space.isnull().sum().data

array(773)

sample_glacier_raster_mag.v_mag_space.notnull().sum().data

array(707)

Next, we will visualize magnitude of velocity calculated by reducing along the x and y dimensions:

fig, ax= plt.subplots()
sample_glacier_raster_mag['v_mag_time'].plot.scatter(marker='o', alpha=0.9)
ax.set_title('Median magnitude of velocity over x,y dimensions, 2017-2020')
ax.set_ylabel('Magnitude of velocity (m/y)')
ax.set_xlabel('Time');

_images/37aeeb0db949d819cc07f3f91b041670de0fecfae1038ba6f3108c671015581c.png

Calculating velocity anomalies#

To do this, we will use xarray .groupby() and .map()

following example from xarray tutorial

We first define a function that subtracts the long-term mean from a single observation.

def remove_time_mean(x):
    return x-x.mean(dim='mid_date')

We then group the dataset by month and apply the function to calculate the anomaly on each group. Because we are using the anomalies of the velocity component vectors we will calculate the mean:

glacier_anom = sample_glacier_raster.groupby('mid_date.month').map(remove_time_mean)

We can visualize the anomalies of the x and y-components of velocity averaged over the glacier surface over time (below). Note that if we tried to compute the magnitude of velocity anomalies, they would all be positive because magnitude is a scalar quantity. Thus, this would show us the magnitude of the anomaly’s deviation from the mean but not the direction. To discern whether a magnitude of velocity anomaly is positive or negative, we would need to examine the signs of the velocity component anomalies. For simplicity, this demonstration will focus on the component anomalies.

fig, ax = plt.subplots(ncols=2, figsize=(18,5))
glacier_anom.mean(dim=['x','y']).vx.plot(ax=ax[0], marker='o', linewidth=0, markersize=5, alpha=0.5)
glacier_anom.mean(dim=['x','y']).vy.plot(ax=ax[1], marker='o', linewidth=0, markersize=5, alpha=0.5)

ax[0].set_title('x-component')
ax[1].set_title('y-component')
ax[0].set_ylabel('m/y')
ax[1].set_ylabel('m/y')
ax[0].set_xlabel('Time')
ax[1].set_xlabel('Time')

fig.suptitle('Spatially-averaged velocity component anomalies');

_images/5ec4939b08901abbca0288443f31ccc9121a25f5fb2bbcfab22c2da402eacc4d.png

Extracting and visualizing data at a single point#

We can use xarray’s .sel() to extract velocity data at a single point or within a subset along given dimensions. In this example, we use .sel() to compare the magnitude of velocity of ice flow at a point in the glacier’s accumulation zone to the mean magnitude of velocit of the entire glacier.

#choose coordinates for a point
y = 3.3125e6
x = 786500

point_v_resample = np.sqrt(glacier_resample_1mo.sel(x=x,y=y,method='nearest').vx**2 + glacier_resample_1mo.sel(x=x,y=y,method='nearest').vy**2)
point_v = np.sqrt(sample_glacier_raster.sel(x=x,y=y,method='nearest').vx**2 + sample_glacier_raster.sel(x=x,y=y,method='nearest').vy**2)

# ax1:
fig,axs=plt.subplots(ncols=3, figsize=(20,5))
#plot raw (not resampled data)
calc_velocity_magnitude(sample_glacier_raster).v_mag_time.plot(ax=axs[0])
point_v.plot(ax=axs[0])
#plot temporal mean and point location
b = calc_velocity_magnitude(glacier_resample_1mo).v_mag_space.plot(ax=axs[1])
axs[1].axvline(x=x, c= 'red')
axs[1].axhline(y=y, c='red')

# spatial mean average velocity over time
calc_velocity_magnitude(glacier_resample_1mo).v_mag_time.plot(ax=axs[2])
# point velocity over time
point_v_resample.plot(ax=axs[2])

fig.suptitle('Comparing velocity at point on glacier (orange) v. spatial average across glacier (blue)')
axs[0].set_title('Raw observations')
axs[1].set_title('Temporal average velocity with point location highlighted')
axs[2].set_title('Monthly resampled observations')

axs[0].set_ylabel('m/yr')
axs[2].set_ylabel('m/yr')
axs[0].set_xlabel('Time')
axs[2].set_xlabel('Time')
b.colorbar.set_label('m/yr');

_images/b7e39e42e0d3e4ed2ca9396ef3620ffed80b89f0c6978ccafda95730fc817089.png

Conclusion#

This notebook demonstrated trimming a full ITS_LIVE granule down to the spatial extent of a single glacier outline. We then worked through some statistical investigation of the data and preliminary data visualization steps. The next notebook will look at grouped analysis of multiple glaciers.