MATLAB Examples

Basal slipperiness and englacial rate factor of the Antarctic Ice Sheet

This example file describes how to import and work with the basal slipperiness and englacial rate factor data from G. Hilmar Gudmundsson.

The dataset is included as a .mat file in this File Exchange submission, but if you'd like to access the raw data and full description, you can find it here.


Load the data

Load the geolocated basal slipperiness (C) and rate factor (A) data like this:

load AntarcticBasalProps.mat

And now you should have variables in your workspace named lat, lon, A, and C. Let's plot 'em! And hopefully you have a relatively recent version of Matlab, because there are over half a million datapoints here, and until a few years ago Matlab used to struggle with large scatter plots.


axis tight              % removes white space
axis off                % removes x, y axes
cb = colorbar;          % places a colorbar
ylabel(cb,'Rate Factor A')
cb.Ruler.Scale = 'log'; % sets colorbar to log scale
caxis([1e-9 1e-8])      % sets colorbar axis limits

And similarly for the slipperiness:

axis tight              % removes white space
axis off                % removes x, y axes
cb = colorbar;          % places a colorbar
ylabel(cb,'Basal slipperiness C')
cb.Ruler.Scale = 'log'; % sets colorbar to log scale
caxis([1e-6 1e-2])      % sets colorbar axis limits

Subsetting to a region

You may have noticed your computer slowing down to plot this large dataset. If you're working on a regional scale, it may help to subset the data and just look within polar stereographic quadrangle. Here's how.

Consider Amery Ice Shelf. Let's say we want to focus on a region 800 km wide by 500 km tall grid centered on Amery. Start by creating a 1 km resolution grid that fills the area.

[X,Y] = psgrid('amery ice shelf',[800 500],1,'xy');

Now get the indices of all the basal properties inside the Amery grid:

ind = inpsquad(lat,lon,X,Y);

How many basal property values are within the Amery grid we created?

ans =

Less than 20 thousand. That's much more manageable than half a million. There are a few ways we can trim down the dataset, so for this example we'll just rewrite each variable as just the data inside the quadrangle:

lat = lat(ind);
lon = lon(ind);
A = A(ind);
C = C(ind);

Ya see, now the dataset's much more manageable:

axis tight off

Gridding and interpolating

From the figure above, you see that the basal property dataset is not on a rectangular grid. This means interpolating won't work until we grid the data ourselves. There are a few ways to grid this data, and one is with the built-in function scatteredInterpolant. Here's how it works.

Begin by converting the geocoordinates to polar stereographic, then create an interpolant. Use the natural neighbor interpolation method.

[x,y] = ll2ps(lat,lon);
F = scatteredInterpolant(x,y,A,'natural');

That's easy enough. And now we can use F to get values of A at any x,y points just by evaluating F(x,y). Let's do that for a grid:

A_gridded = F(X,Y);

% mask out the open ocean:
A_gridded(isopenocean(X,Y)) = nan;

% Plot the gridded data:
shading interp
axis image off % removes white space, sets equal aspect ratio, removes border

Now we have the dataset nicely gridded. Whenever gridding data, I always use the AMT shadem function with the user interface mode to see potential gridding artifacts. Although it may not seem intuitive or physically meaningful to apply hillshade to something like a rate factor, shading a gridded surface can highlight areas of overfitting (usually identified by local pits and spikes) or over smoothing (smeared out areas). Here's what I mean:


In the image above, I definitely see artifacts of gridding, but it's actually from the Bedmap2 data that informed Gudmundsson's data. You see those distinct lines of lumps that look like a pearl necklace? Those are flight lines that fed Bedmap2. And I see some smeared out areas, indicating regions of missing ice thickness data. We don't see any hexagonal patterns from gridding the basal property data itself.

How I created the AntarcticBasalProps.mat file

The dataset in AntarcticBasalProps.mat was generated by G.H. Gudmundsson and can be found here. It's called AC.dat and here's how I converted it into the .mat file:

%  D = importdata('AC.dat');
%  lat = D(:,1);
%  lon = D(:,2);
%  A = D(:,5);
%  C = D(:,6);
%  readme = 'Gudmundsson''s basal properties from';
%  save('AntarcticBasalProps.mat','lat','lon','A','C','readme');

Note that I did NOT save the x or y data present in AC.dat, because although it is described as being in polar stereographic coordinates, the exact parameters of the coordinate transformation are unknown and they don't match the NSIDC or Antarctic Mapping Tools standard. The difference is small, about a 68 m rms difference compared to the NSIDC and Antarctic Mapping Tools standard. I suspect Gudmundsson used a nonstandard value for Earth's eccentricity when converting to polar stereographic coordinates, so for consistency I stick with geo coordinates because they are absolute, then convert to standard polar stereographic as needed.


The dataset is from Gudmundsson 2013, built on Rignot 2011 velocity and Fretwell ice geometry. If this example file or anything in Antarctic Mapping Tools is useful to you, please cite my AMT paper (Greene et al., 2017).

1. Rignot, E., Mouginot, J. & Scheuchl, B. Ice flow of the Antarctic ice sheet. Science 333, 1427-1430 (2011). DOI:10.1126/science.1208336.

2. Fretwell, P. et al. Bedmap2: Improved ice bed, surface and thickness datasets for Antarctica. Cryosphere 7, 375-393 (2013). DOI:10.5194/tc-7-375-2013.

3. Gudmundsson, G. H. Ice-shelf buttressing and the stability of marine ice sheets. Cryosph. 7, 647-655 (2013). DOI:10.1073/pnas.1512482112.

4. Greene, C. A., Gwyther, D. E., Blankenship, D. D. Antarctic Mapping Tools for Matlab. Computers and Geosciences 104, 151-157 (2017). doi:10.1016/j.cageo.2016.08.003.