MATLAB Examples

flowline documentation

The flowline function predicts ice flowpath(s) based on MEaSUREs InSAR-derived surface velocities.



[lat,lon] = flowline(latstart,lonstart)
[lat,lon] = flowline(xstart,ystart)
[lat,lon] = flowline(...,DomainSizekm)
[lat,lon] = flowline(...,'plotgeo',LineProperty,LineValue,...)
[lat,lon] = flowline(...,'plotxy',LineProperty,LineValue,...)
[lat,lon,v,t,d,h] = flowline(...)
h = flowline(...,'plot...)


[lat,lon] = flowline(latstart,lonstart) calculates flow path(s) from seed locations given by geographic coordinate(s) latstart,lonstart. If multiple starting points are specified, output lat and lon will be cell arrays.

[lat,lon] = flowline(xstart,ystart) calculates flow path(s) from seed locations given by polar stereographic coordinate(s) xstart,ystart. Polar stereographic coordinates (ps71) are automatically determined by the islatlon function.

[lat,lon] = flowline(...,DomainSizekm) specifies a domain size in kilometers. Bigger domains take longer to solve. Smaller domains could potentially clip the edge of data. Just make sure your domain size captures an the entire drainage basin(s) of your area of interest and you'll be good. Default DomainSizekm is 2000, which takes about 4 seconds on my laptop. Processing time is roughly proportional to the domain size squared, so a domain size of 1000 km takes about 1 second on my laptop.

[lat,lon,v,t,d,h] = flowline(...) returns corresponding speed v (m/yr), time t (yr), and distance d (m) vectors for each flow line. Time and distance are measured from the starting location(s). The optional output argument h is a handle of plotted objects. Calculation of v, t, and d can slow things down a bit, especially for large grids of many starting locations.

h = flowline(...,'plot...) returns a handle h of plotted flowlines.


This function requires Antarctic Mapping Tools and the MEaSUREs Toolbox.

Example 1: Plot flow lines from gridded seed locations

Perhaps you want to seed the stream line calculation with multiple starting points. Here we use psgrid to create a 500 km wide grid of starting points centered on Pine Island Glacier at 25 km resolution, then solve for the flowpaths.

[latstart,lonstart] = psgrid('pine island glacier',450,25);

[lat,lon] = flowline(latstart,lonstart);

You'll notice that you now have lat and lon, which are 19x19 cell arrays in your workspace. Each cell in those arrays contains the coordinates of a single flowline. If there's only one coordinate in the cell, it means there's no ice motion in that location.

If you want to plot results when you call flowline, syntax is easy:


If you have Matlab's Mapping Toolbox and you prefer to plot in geo coordinates, you can do that too. Also, any arguments that follow 'plotxy' or 'plotgeo' are interpreted as linestyle specification. So below I'll plot red lines on top of a MODIS MOA image of Pine Island Glacier. Also, our domain is pretty small, so I'll specify a DomainWidthkm of 800 to save a couple seconds of processing time.

modismoa('pine island glacier','inset','ne')

Example 2: Travel time and velocity along flow

Sometimes you might want to know not only the path a parcel of ice might take, but also its speed or travel time as a function of distance traveled. The flowline function returns path distance in meters, speed in meters per year, and time in years.

Consider the Lambert Glacier/Amery Ice Shelf system. Below I'm using ramp to plot a Radarsat Mosaic, then I'm overlaying a semitransparent ice speed layer with measures. We have to call freezeColors to allow the colorful ice speed layer to exist in the same axes as the ramp image. ramp, measures and freezeColors are all available on the Mathworks File Exchange website. I'll mark a starting location at (74°S,65.2°E) with a blue circle:

mapzoom('lambert glacier',1200,'inset','nw')

measures('speed','lambert glacier','mapwidth',1200,...


Calculate and plot the flow path:

[lat,lon,v,t,d] = flowline(-74.6,65.2,'plotgeo','b-','linewidth',2);

Plot speed as a function of distance traveled, and also plot speed as a function of travel time. Travel distances are divided by 1000 to convert from meters to kilometers:

box off; axis tight
xlabel 'distance traveled (km)'
ylabel 'speed (m/yr)'

box off; axis tight
xlabel 'travel time (years)'
ylabel 'speed (m/yr)'

Perhaps you don't work with time scales of thousands of years. If you would like to compare a 1 year long GPS record to the flow path predicted by MEaSUREs, simply perform 1D interpolation. Let's assume you have daily GPS records for 1 year. First convert t to days to make the math easy:

t_days = t*365.24;

Interpolate time, distance, and speed to 1:365 days to get the first year of daily predicted motion. There are some undefined points near the calving front of Amery Ice Shelf, so we'll ignore them in the interp1 calculation by including only finite data:

ti = 1:365;
di = interp1(t_days(isfinite(t)),d(isfinite(t)),ti);
vi = interp1(t_days(isfinite(t)),v(isfinite(t)),ti);
lati = interp1(t_days(isfinite(t)),lat(isfinite(t)),ti);
loni = interp1(t_days(isfinite(t)),lon(isfinite(t)),ti);

The first year of GPS data should look something like this:

box off; axis tight
xlabel 'distance traveled (m)'
ylabel 'speed (m/yr)'

box off; axis tight
xlabel 'travel time (years)'
ylabel 'speed (m/yr)'

That's right, this is a pretty slow-moving area, so we don't expect any a parcel of ice to noticeably speed up or slow down in the course of a year.

If you plot your GPS data in polar stereographic meters, compare it to MEaSUREs-predicted drift tracks like this:

% This is the predicted path at daily resolution:
[xi,yi] = ll2ps(lati,loni);

% And let's say this is your daily GPS data:
xGPS = xi+0.1*cumsum([0,randn(1,364)]);
yGPS = yi+0.1*cumsum([0,randn(1,364)]);

% Compare predicted to measured tracks:
hold on
legend('InSAR-predicted ice flow path',...
    'one year of GPS data','location','southwest')
box off
legend boxoff
xlabel 'polar stereographic eastings (m)'
ylabel 'polar stereographic northings (m)'

Author Info

The flowline function and supporting documentation were written by Chad A. Greene of the University of Texas at Austin's Institute for Geophysics (UTIG), 2015.