MATLAB Examples

measures_interp documentation

The measures_interp interpolates MEaSUREs ice surface velocity data at any location(s), along a path, or onto a new grid.

Contents

Before you use this function:

This function requires the MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2 netcdf dataset which can be downloaded here: https://nsidc.org/data/NSIDC-0484.

Syntax

speedi = measures_interp('speed',lati_or_xi,loni_or_yi)
[vxi,vyi] = measures_interp('velocity',lati_or_xi,loni_or_yi)
err = measures_interp('error',lati_or_xi,loni_or_yi)
count = measures_interp('count',lati_or_xi,loni_or_yi)
[u,v] = measures_interp('uv',lati_or_xi,loni_or_yi)
[along,across] = measures_interp('track',lati_or_xi,loni_or_yi)
[...] = measures_interp(...,'method',interpMethod)
[...] = measures_interp(...,'fill')

Description

speedi = measures_interp('speed',lati_or_xi,loni_or_yi) returns local surface velocity at the location(s) given by lati,loni or xi,yi. Input coordinates are determined as geographic or polar stereographic automatically by the islatlon function.

[vxi,vyi] = measures_interp('velocity',lati_or_xi,loni_or_yi) returns the polar stereographic (true latitude 71S) x and y components of velocity.

err = measures_interp('error',lati_or_xi,loni_or_yi) returns a scalar value of uncertainty estimates presented with the dataset.

count = measures_interp('count',lati_or_xi,loni_or_yi) returns the count of scenes used per pixel.

[u,v] = measures_interp('uv',lati_or_xi,loni_or_yi) returns the geographic zonal (positive eastward) and meridional (positive northward) components of velocity.

[along,across] = measures_interp('track',lati_or_xi,loni_or_yi) returns velocity components along or across a specified track. This can be used to estimate flow through (the across component) a flux gate or grounding line. For the across-track component, positive values are to the right.

[...] = measures_interp(...,'method',interpMethod) specifies any interpolation supported by interp2. Default interpMethod is 'linear'.

Example 1: South Pole station drift

For a single location, interpolation is easy:

measures_interp('speed',-90,0)
ans =
          7.32

which of course is exactly the same as

measures_interp('speed',-90,100)
ans =
          7.32

and is very similar to

measures_interp('speed',-90,0,'method','cubic')
ans =
          7.25

But differences in interpolation methods are tiny compared to the uncertainty of InSAR or any satellite-based surface velocity measurements. The self-reported uncertainty at the South Pole is about 11 m/yr, see:

measures_interp('error',-90,100)
ans =
         11.01

Example 2: Ice speed on a custom grid

You may be working with some other data set that has its own lat/lon grid. Here we consider a 500 km wide grid at 300 m resolution centered on the Siple Coast:

% Create a grid:
[x,y] = psgrid('siple coast',500,0.3,'xy');

% Get speed data at each grid point:
speed = measures_interp('speed',x,y);

% Plot:
pcolor(x,y,log10(speed))
shading interp
axis equal
hold on

It might also be nice to overlay some velocity vectors. Simply use Matlab's built-in quiver function. But before we call quiver, let's downsample the dataset because 1667x1667 little arrows would just end up looking like a big black square, as each arrow would have to be less than a pixel wide. There are a few ways to downsample. You can simply create a much more coarse grid with psgrid, however there would be one tiny problem with that, which is aliasing--probably not much of an issue, but a proper way to deal with aliasing is to perform a low-pass filter before sampling. If you have Matlab's Image Processing Toolbox, the imresize function automatically performs antialiasing that way, so we'll use imresize and to scale the 300 m grid to 2% of its resolution, or 15 km resolution. Like this:

% Get velocity components:
[vx,vy] = measures_interp('velocity',x,y);

% Plot quiver arrows on a 15 km grid:
quiver(imresize(x,0.02),imresize(y,0.02),imresize(vx,0.02),imresize(vy,0.02),'k');

Example 3: Flux gate calculation of Thwaites Glacier mass loss

To calculate mass flux across a grounding line, we first define a grounding line. Here we use the grounding line from the Bedmap2 Toolbox and clip the data to include only the region of Thwaites Glacier:

load bedmap2gl

gllat = flipud(gllat{1});
gllon = flipud(gllon{1});
gllat = gllat(gllon>-108.5&gllon<-104.5);
gllon = gllon(gllon>-108.5&gllon<-104.5);

For context, we can plot this grounding line on a lima image with velocity vectors overlaid.

% Initialize a base map:
figure
lima('thwaites glacier',150,'xy')
hold on

% Overlay black velocity vectors:
measuresps('vel','k')

% Plot the segment of grounding line we're using:
plotps(gllat,gllon,'g','linewidth',2)

% Add text labels:
textps(gllat(1),gllon(1),'start','color','g',...
    'fontweight','bold','backgroundcolor','w')
textps(gllat(end),gllon(end),'end','color','r',...
    'fontweight','bold','backgroundcolor','w')

If we think of the grounding line as a path along which you might walk, the only component of ice flow that contributes to continental mass loss is the cross-path component. To estimate the volume of ice crossing the grounding line, interpolate the cross-track component of ice velocity along the entire path of the grounding line, then multiply velocity by thickness for each unit length along the grounding line.

First we define x as the distance you'd walk along the grounding line and dx is the approximate distance between points along the grounding line. Here we use the pathdistps function to calculate the distance traveled along the grounding line.

% Distance along grounding line:
d = pathdistps(gllat,gllon);

% Distance between each grounding line data point:
dx = gradient(d);

crossTrackVelocity = measures_interp('cross',gllat,gllon);

thickness = bedmap2_interp(gllat,gllon,'thickness');
thickness(isnan(thickness))=0; % zero thickness where undefined

flowAcrossGL = crossTrackVelocity.*thickness;

figure
plot(d/1000,flowAcrossGL)
xlabel('distance along grounding line (km)')
ylabel('flow across gl (m^3/yr per meter along grounding line)')
box off; axis tight;

Above, we see that sometimes ice flow is negative--that happens where a sinuous grounding line lets ice go over the ocean, then reground, then cross the grounding line again. We can easily distill all this rich information down to a single value of mass loss if we ignore firn density and say that everything flowing across the grounding line is pure ice. massBalanceGT is taken as the negative to indicate mass loss and multiplied by 1e-12 to convert from kg to GT.

totalVolFlow = sum(flowAcrossGL.*dx);
iceDensity = 917; % kg/m3

massBalanceGT = totalVolFlow*iceDensity*1e-12
massBalanceGT =
        113.97

This value is in close agreement with 113.5 GT/yr found by Rignot et al., 2013.

Citing these datasets

VELOCITY DATA: Rignot, E., J. Mouginot, and B. Scheuchl. 2017. MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2. [Indicate subset used]. Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. doi: http://dx.doi.org/10.5067/D7GK8F5J8M8R.

A LITERARY REFERENCE FOR THE VELOCITY DATA: Rignot, E., J. Mouginot, and B. Scheuchl. 2011.Ice Flow of the Antarctic Ice Sheet, Science, Vol. 333(6048): 1427-1430. doi:10.1126/science.1208336.

ANTARCTIC MAPPING TOOLS: Greene, C.A., Gwyther, D.E. and Blankenship, D.D., 2016. Antarctic Mapping Tools for Matlab. Computers & Geosciences. http://dx.doi.org/10.1016/j.cageo.2016.08.003

File history

July 2014: First version written.

August 2014: Updated as a plugin for Antarctic Mapping Tools.

October 2016: Fully rewritten--Now reads data from the new measures_data function, which reads the .nc file directly.

May 2017: Updated for data version 2.

Author Info

This function was written by Chad A. Greene of the University of Texas Institute for Geophysics (UTIG), July 2014. Rewritten October 2016 for efficiency and usability.