MATLAB Examples

Plotting post-glacial rebound (PGR)

Post-glacial rebound (also known as glacial isostatic adjustment or GIA) describes the elastic relaxation of land masses that has been occurring since the weight of continental ice sheets began to wane some 20,000 years ago. This example describes how to get PGR data and plot it in Matlab. The Mapping Toolbox makes plotting intuitive, but examples are provided for users without the Mapping Toolbox.

Contents

Data Download

PGR data from the GRACE Tellus project are described here. This example downloads current PGR trends using a 200-km smoothing radius, 60-degree filter. If you prefer a different smoothing radius, or even geoid or mass data, edit the downloaded file name accordingly. In the steps below, PGR data are downloaded as a text file, which is then read into Matlab and re-saved as a .mat file in gridded format.

urlwrite('ftp://podaac-ftp.jpl.nasa.gov/allData/tellus/L3/pgr/GIA_n60_uplift_200km.txt',...
    'GIA_n60_uplift_200km.txt') ;

fid = fopen('GIA_n60_uplift_200km.txt');
GIA = textscan(fid,'%f %f %f','headerlines',18);
fclose(fid);

lat = flipud(reshape(GIA{2},[180 360]));
lon = flipud(reshape(GIA{1},[180 360]));
gia = flipud(reshape(GIA{3},[180 360]));

clear GIA
save gia_n60_uplift_200km

Plot PGR without the Mapping Toolbox

After downloading and re-saving the data as described above, you should have a file called gia_n60_uplift_200km.mat. Plotting the data without the Mapping Toolbox is straightforward. First we'll load the .mat file, then we'll plot it. The only tricky part is remembering that although we tend to say "lat, lon", we will actually need to plot "lon, lat" because longitude describes location on the horizontal axis.

load gia_n60_uplift_200km

figure
pcolor(lon,lat,gia)
shading interp
xlabel('longitude')
ylabel('latitude')
cb = colorbar;
ylabel(cb,'PGR mm/yr')

Plot PGR with the Mapping Toolbox

If you have Matlab's Mapping Toolbox, plotting PGR is quite similar to the steps above.

load gia_n60_uplift_200km

figure
worldmap('world') % initializes map
pcolorm(lat,lon,gia)
xlabel('longitude')
ylabel('latitude')
cb = colorbar('location','southoutside');
xlabel(cb,'PGR mm/yr')

Customizing the map

The map we just plotted is nice, but it leaves a little something to be desired. For example, a seam in the data is visible at the prime meridian, a colormap divergent about zero would better depict where land masses are moving up versus down, a colormap with more than 64 colors would bring the graphics into the 21st century, and the ghostly outlines of continents are cool, but we may want to define them more clearly.

First we'll start with the seam problem. The map looks a bit like my Uncle Roy's belly testing the tensile strength of his shirt buttons after Thanksgiving dinner. Let's give the Earth a slightly larger shirt by overlapping some data. For those of you who are just joining us, we'll start at the step where we load the data which was previously saved in the Data Download section.

load gia_n60_uplift_200km

lon(:,361)=lon(:,1);
lat(:,361)=lat(:,1);
gia(:,361)=gia(:,1);

figure
worldmap('world') % initializes map
pcolorm(lat,lon,gia)
xlabel('longitude')
ylabel('latitude')
cb = colorbar('location','southoutside');
xlabel(cb,'PGR mm/yr')

The seam is gone! For a simple colormap divergent about zero, we can use the b2r function available on the File Exchange site here:

colormap(b2r(-4,14))

In addition to making the colormap divergent about 0 mm/yr PGR rate, the b2r function also created a smooth color gradient with many more than 64 colors. Alternatively, if you have my rgbmap function, you can create a colormap similar to the one we just created; however, the rgbmap function does not automatically center the colormap about zero so you'd have to do it manually like this

caxis([-14 14])
rgbmap('blue','white','red')

Now you may wish to show coastlines. We can do that using the built-in coast data set. Note that loading the coast data set will overwrite lat and create a new long variable.

load coast
plotm(lat,long,'k')

Interpolate PGR using pgr_interp

The interpolation function does not require the Mapping Toolbox. Simply feed pgr_interp a single location, track line, or grid of lat/lon points and it'll return interpolated PGR rates at those points. Note that this interpolation assumes that the underlying PGR data are gridded equally in space, however in reality they are gridded equally in lat/lon coordinates. Given the broad smoothing radius of this data, the difference should be negligible.

Now to answer the question that been on all of our minds, what is the current rate of post-glacial rebound in Ukkusiksalik National Park, Canada?

UkkusiksalikLat = 65.34;
UkkusiksalikLon = -87.31;

pgr_interp(UkkusiksalikLat,UkkusiksalikLon)
ans =

    5.6568

That Laurentide Ice Sheet was nothing to shake a stick at! Ukkusiksalik National Park is still rising over 5 mm per year there.

Author Info

This example file and corresponding pgr_interp function were created by Chad A. Greene of the Institute for Geophysics at the University of Texas at Austin, July 2014. I claim no credit for the PGR data retrieved by this script.