MATLAB Examples

Jacobs2013Getz

Jacobs et al. did a nice job of combining several data sets in their 2013 JGR Oceans paper titled "Getz Ice Shelf melting response to changes in ocean forcing." Figure 4 of their paper may look a bit dry at first glance, but it so nicely illustrates how small changes in the depth of a thermocline can have the capacity to affect a tremendous area of an ice shelf's base, depending on where the thermocline lies in relation to the vertical distribution of ice shelf basal area. In this example we recreate some of the figures from the Jacobs et al. paper.

Contents

Requirements

This example requires

Map the region

Make a map of the region similar to the Jacobs et al. Figure 1. Start by initializing an 800 kilometer wide map centered on Getz Ice Shelf. Plot bed elevations from the Bedmap2 dataset, and if you'd like you can approximate the Jacobs et al. colormap with rgbmap.

figure(1)

% Zoom an 800 km wide map to Getz Ice Shelf and place inset in upper right corner:
mapzoom('getz ice shelf',800,'inset','ne')

% Plot bed elevations:
bedmap2('bed','colorbar','south')

% Set color axis limits:
caxis([-1000 0])

% Mimic Jacobs et al. colormap:
rgbmap('darkish blue','aqua','yellow','red')

% Plot grounded ice and ice shelves as patch objects:
bedmap2('patchshelves','facecolor','w')
bedmap2 patchgl

% Place a grid of lats and lons:
antmap('graticlue','lats',-79:-69,'lons',-135:5:115,'color',.4*[1 1 1])

% Label locations of interest:
places = {'Martin Peninsula','Wright Island',...
    'Duncan Peninsula','Carney Island','Siple Island',...
    'Dean Island','Grant Island'};
scarlabel(places,'fontsize',10)

Plot CTD locations

Plot all the oceanographic stations in the Southern Ocean Database like this:

sodb('stations','ko','markerfacecolor','blue')

You may notice that many of these loations do not exactly match Figure 1 of the Jacobs et al. paper. That's because Jacobs et al. showed stations from a few select years, while the SODB has stations from many years, but most SODB data in this region are from the 1990's.

We can make something similar to the inset map in Jacobs et al.'s Figure 1 like this. Below latnbp,lonnbp is a crude estimate of the Nathaniel B. Palmer's cruise track which I defined by tracing approximate locations with inputm. The fakedata.mat file is included in the Glaciology Fake Book.

figure(11)

% Create 3500 km wide map of West Antarctica:
mapzoom('west antarctica',3500)

% Mark the continental slope with a single contour at 1500 m depth:
bedmap2('bedc','zvals',-1500)

% Show ice shelves and grounded ice as patch objects:
bedmap2 patchshelves
bedmap2 patchgl

% Plot an approximate Palmer cruise track:
load fakedata.mat latnbp lonnbp
plotm(latnbp,lonnbp,'r','linewidth',3)

Interpolate along a cruise track

Jacobs et al.'s Figure 2 shows measured temperature and salinity profiles along a Palmer cruise track. We don't have measured data along that track, but we can make a plot similar to Figure 2 with sodb_profile, which interpolates SODB gridded data along any arbitrary transect. Below we'll specify contour levels to match Jacobs et al.'s Figure 2 and make all contour lines black with the 'k' option.

figure(2)
set(gcf,'pos',[100 100 750 800])

% Plot potential temperature profile:
subplot(211)
sodb_profile('ptm',latnbp,lonnbp,'contour',[-1.8 -1.5 -1 -0.5 0 0.5 1],'k','label')
ylim([-1200 0])

% Plot salinity profile:
subplot(212)
sodb_profile('sal',latnbp,lonnbp,'contour',33.8:.2:34.8,'k','label')
ylim([-1200 0])

You'll see that our interpolated data are discontinuous because the gridded SODB dataset is undefined close to the coast and near major bathymetric features, but the contour plots above contain some of the major features we want to see. Check out the -1.8° loop at about 2800 km along track, and the sloping 34.8 salinity line from 3800 km to the right-hand edge of the domain.

Get select CTD station data

Jacobs et al. distinguish between stations along the ice front versus the outer continental shelf, and they also distinguish between eastern and western stations. I used inputm to freehand-select some bounding boxes for each of the four regions and I got these locations for the bounding boxes:

% ice front east:
ifelat = [ -74.03  -73.72  -73.70  -73.82  -74.53  -74.43];
ifelon = [-114.24 -115.53 -117.70 -119.96 -119.81 -115.69];

% ice front west:
ifwlat = [ -74.68  -73.73  -73.89  -74.23  -75.06];
ifwlon = [-126.39 -126.64 -132.99 -135.37 -133.06];

% outer continental shelf east:
ocelat = [ -72.95  -72.19  -71.66  -72.29  -72.98  -73.23];
ocelon = [-112.99 -112.16 -117.49 -125.26 -124.47 -113.75];

% outer continental shelf west:
ocwlat = [ -73.36  -72.84  -72.72  -74.07  -74.44];
ocwlon = [-127.39 -126.83 -128.02 -137.70 -136.47];

Plot the bounding boxes like this. I'm using blue for eastern stations and my rgb function to match the green of the Jacobs et al. paper for western locations:

figure(1)
plotm(ifelat,ifelon,'b-o','linewidth',2)
plotm(ifwlat,ifwlon,'-o','color',rgb('green'),'linewidth',2)
plotm(ocelat,ocelon,'b-o','linewidth',2)
plotm(ocwlat,ocwlon,'-o','color',rgb('green'),'linewidth',2)

Get data from stations within each of the four polygons defined above using sodb_inpolygon. You can either define the lats and lons in script as we're doing here, or you can call sodb_inpolygon without any input arguments to manually select polygons via mouse clicks. Give this a few seconds to process:

ife = sodb_inpolygon(ifelat,ifelon);
ifw = sodb_inpolygon(ifwlat,ifwlon);
oce = sodb_inpolygon(ocelat,ocelon);
ocw = sodb_inpolygon(ocwlat,ocwlon);

Jacobs et al. invesitgated stations from the years 2000 and 2007. Unfortunately, the SODB does not contain these datasets. But we have some data collected in other years--primarily 1994. Check out the data collection times from the stations in the ice front east polygon:

datestr(ife.t)
ans =

01-Mar-1994 06:52:48
01-Mar-1994 11:28:00
01-Mar-1994 20:04:00
01-Mar-1994 05:29:36
01-Mar-1994 10:28:00
01-Mar-1994 14:29:36
01-Mar-1994 18:16:00

Station data profiles

So we cannot exactly recreate the Jacobs et al. results because the SODB does not have all the same data, but we can still make a plot similar to Jacobs et al.'s Figure 4. We'll start with Figure 4a, which shows temperature and salinity profiles from the ice front polygons. One way to do this is to loop through each station. We'll start with salinity:

figure(41)
hold on

% For each ice front east station, plot salinity:
for k = 1:length(ife.sta)
    plot(ife.sal{k},ife.prs{k},'.','color',rgb('gray'),'markersize',7)
end

% Now do the same for ice front west:
for k = 1:length(ifw.sta)
    plot(ifw.sal{k},ifw.prs{k},'.','color',rgb('gray'),'markersize',7)
end

% Format the axes:
axis([33.5 34.8 0 1200])
set(gca,'ydir','reverse')
xlabel('salinity')

There are a few ways to plot a second x axis. Here we'll get the position gcap of the current axes, then create a new set of axes in the exact same position and set the background color of the new axes to 'none' so it won't cover up our salinity plot:

% Get position of salinity axes:
gcap = get(gca,'pos');

% Create new set of axes atop salinity axes:
axes('pos',gcap,'color','none')
hold on

% Ice front east potential temperature:
for k = 1:length(ife.sta)
    plot(ife.ptm{k},ife.prs{k},'b.','markersize',7)
end

% Ice front west potential temperature:
for k = 1:length(ifw.sta)
    plot(ifw.ptm{k},ifw.prs{k},'.','color',rgb('green'),'markersize',7)
end

set(gca,'XAxisLocation','top','ydir','reverse')
axis([-2 1.5 0 1200])
ylabel 'Pressure (dbar)'
xlabel 'Temperature ({\circ}C)'

text(0.4,100,'1994E','color','b','fontweight','bold')
text(0.4,200,'1994W','color',rgb('green'),'fontweight','bold')

To recreate Figure 4b we'll do pretty much the same thing as for 4a, but using outer continental shelf stations oce and ocw:

figure(42)
hold on

% Salinities:
for k = 1:length(ife.sta)
    plot(oce.sal{k},oce.prs{k},'.','color',rgb('gray'),'markersize',7)
end
for k = 1:length(ifw.sta)
    plot(ocw.sal{k},ocw.prs{k},'.','color',rgb('gray'),'markersize',7)
end
axis off
axis([33.5 34.8 0 600])
set(gca,'ydir','reverse')
gcap = get(gca,'pos');

% Temperatures:
axes('pos',gcap,'color','none')
hold on
for k = 1:length(oce.sta)
    plot(oce.ptm{k},oce.prs{k},'b.','markersize',7)
end
for k = 1:length(ocw.sta)
    plot(ocw.ptm{k},ocw.prs{k},'.','color',rgb('green'),'markersize',7)
end
set(gca,'ydir','reverse')
axis([-2 1.5 0 600])
ylabel 'Pressure (dbar)'
xlabel 'Temperature ({\circ}C)'

T-S diagram

T-S diagrams are easy to make in Matlab. We can loop through each station as above, or, since we just want a simple scatter plot without preserving information about each station, we can concatenate data from multiple stations into single arrays like this:

oceptm = cell2mat(oce.ptm(:));
ocesal = cell2mat(oce.sal(:));
ocwptm = cell2mat(ocw.ptm(:));
ocwsal = cell2mat(ocw.sal(:));

Plotting concatenated data is a little bit more elegant can be used to mimic Figure 4c like this:

figure(43)
plot(ocesal,oceptm,'b.','markersize',7)
hold on
plot(ocwsal,ocwptm,'.','color',rgb('green'),'markersize',7)
axis([33.5 34.8 -2 1.5])
xlabel 'Salinity'
set(gca,'YAxisLocation','right')
ylabel 'Temperature ({\circ}C)'

Ice shelf basal elevation histogram

Figure 4d plots a histogram of basal elevations of Getz Ice Shelf. To do this, start by defining an area of interest. We'll need an outline of Getz ice shelf, and the Bedmap2 Toolbox function outlineashelf can give us just that:

[~,getzlat,getzlon] = outlineashelf('getz');

Load Bedmap2 DEM data

With the bedmap2_data function we can import the entire Bedmap2 dataset, but loading an entire continent's worth of data would be unnecessarily computationally expensive. If you enter some array of latitudes and longitudes into the bedmap2_data function, it'll only load enough Bedmap2 data to fully encompass the geolocations you entered. That's what we'll do.

There is no Bedmap2 basal elevation dataset, so we will have to build one by subtracting a Bedmap2 thickness from the Bedmap2 surface. Let's do that now:

% Load gridded Bedmap2 data encompassing Getz:
[lat,lon,sfz] = bedmap2_data('surface',getzlat,getzlon);
[~,~,thck] = bedmap2_data('thickness',getzlat,getzlon);

% Basal elevations = surface - thickness:
base = sfz-thck;

Note that Matlab sometimes gets confused if a variable is named sfc, so above we avoid this potential problem by naming our surface elevations sfz.

Mask Bedmap2 DEM data from Getz outline

The outline of Getz Ice Shelf is given by getzlat,getzlon. To find only the basal elevation DEM datapoints within this outline, use inpolygon:

getzind = inpolygon(lat,lon,getzlat,getzlon);

Switching to polar stereographic coordinates

Now getzind contains indices of all the 31,860 gridded datapoints within the outline given by getzlat,getzlon. However! We have given Matlab the impression that one degree of latitude is equal in distance to one degree of longitude. By that I mean, if you place the lat/lon coordinates of the Getz Ice Shelf outline on an equally-spaced lat/lon set of axes, the shape of the polygon will be distorted.

So for best accuracy, we switch from lat/lon coordinates to polar stereographic coordinates with ll2ps:

[getzx,getzy] = ll2ps(getzlat,getzlon);

Just for fun, let's compare the shapes of the Getz polygon in equally-spaced lat/lon coordinates with equally-spaced polar sterographic coordinates:

figure
subplot(211)
patch(getzlon,getzlat,'b')
xlabel('degrees longitude')
ylabel('degrees latitude')
axis tight equal
box off

subplot(212)
patch(getzx,getzy,'b')
xlabel('meters easting')
ylabel('meters northing')
axis tight equal
box off

Of course they're a bit rotated relative to each other, and that doesn't matter much--what matters is that the shape of the polygon is distorted when we imply that lats and lons are equally spaced, so for accuracy, it's preferable to use inpolygon with polar stereographic coordinates. And that is what we will do:

[x,y] = ll2ps(lat,lon);
getzindps = inpolygon(x,y,getzx,getzy);

All the basal elevations of Getz Ice Shelf can be collected into an array like this:

getzBaseElevation = base(getzindps);

A histogram of grid cells corresponding to specific depths can be obtained with hist:

depth = -1200:20:0; % meters depth roughly equates to dbar pressure
ncount = hist(getzBaseElevation,depth);

A couple of convenient bits that make this math easy:

  1. meters depth in the ocean roughly corresponds to dbar pressure, and
  2. The Bedmap2 data set is a 1 km by 1 km grid, so a count of raster data points corresponds to square kilometers of area.
% Here's how we can make Figure 4d:

figure(44)
barh(depth,ncount,...
    'facecolor',.6*[1 1 1],...
    'edgecolor',.5*[1 1 1])
axis([0 2400 -1200 0])
xlabel('Area (km^2)')
set(gca,'xaxislocation','top','yaxislocation','right')
ylabel 'depth (m) ~dbar'

Figure 4 in its entirety

Putting all the subplots together looks like this:

figure(4)

% Figure 4a:
subplot(3,2,[1 3])
hold on
for k = 1:length(ife.sta)
    plot(ife.sal{k},ife.prs{k},'.','color',rgb('gray'),'markersize',7)
end
for k = 1:length(ifw.sta)
    plot(ifw.sal{k},ifw.prs{k},'.','color',rgb('gray'),'markersize',7)
end
axis([33.5 34.8 0 1200])
set(gca,'ydir','reverse')
xlabel('salinity')
gcap = get(gca,'pos');
axes('pos',gcap,'color','none')
hold on
for k = 1:length(ife.sta)
    plot(ife.ptm{k},ife.prs{k},'b.','markersize',7)
end
for k = 1:length(ifw.sta)
    plot(ifw.ptm{k},ifw.prs{k},'.','color',rgb('green'),'markersize',7)
end
set(gca,'XAxisLocation','top','ydir','reverse')
axis([-2 1.5 0 1200])
ylabel 'Pressure (dbar)'
xlabel 'Temperature ({\circ}C)'
text(0.4,100,'1994E','color','b','fontweight','bold')
text(0.4,200,'1994W','color',rgb('green'),'fontweight','bold')
box on

% Figure 4b:
subplot(3,2,5)
hold on
for k = 1:length(ife.sta)
    plot(oce.sal{k},oce.prs{k},'.','color',rgb('gray'),'markersize',7)
end
for k = 1:length(ifw.sta)
    plot(ocw.sal{k},ocw.prs{k},'.','color',rgb('gray'),'markersize',7)
end
axis off
axis([33.5 34.8 0 600])
set(gca,'ydir','reverse')
gcap = get(gca,'pos');
axes('pos',gcap,'color','none')
hold on
for k = 1:length(oce.sta)
    plot(oce.ptm{k},oce.prs{k},'b.','markersize',7)
end
for k = 1:length(ocw.sta)
    plot(ocw.ptm{k},ocw.prs{k},'.','color',rgb('green'),'markersize',7)
end
set(gca,'ydir','reverse')
axis([-2 1.5 0 600])
ylabel 'Pressure (dbar)'
xlabel 'Temperature ({\circ}C)'
box on

% Figure 4c:
subplot(3,2,6)
plot(ocesal,oceptm,'b.','markersize',7)
hold on
plot(ocwsal,ocwptm,'.','color',rgb('green'),'markersize',7)
axis([33.5 34.8 -2 1.5])
xlabel 'Salinity'
set(gca,'YAxisLocation','right')
ylabel 'Temperature ({\circ}C)'

% Figure 4d:
subplot(3,2,[2 4])
barh(depth,ncount,...
    'facecolor',.6*[1 1 1],...
    'edgecolor',.5*[1 1 1])
axis([0 2400 -1200 0])
xlabel('Area (km^2)')
set(gca,'xaxislocation','top','yaxislocation','right')
ylabel 'depth (m) ~dbar'
set(gcf,'pos',[10 10 650 800])

Author Info

Chad Greene wrote this example file. Chad claims absolutely no credit for the data, the science, or the creative display shown above.