MATLAB Examples

How to make a map of sea level rise

This example file describes how anyone, even you! can create a map of sea level rise using data from the University of Colorado's Sea Level Research Group.

Contents

Download data

The first thing you'll need to do is download the text data from this website and put in your current folder or some other folder where Matlab can find it. Note that the text data is a large file. Instead of clicking on the link you can simply right-click and save the data directly.

Import data

Now we'll use textscan to import the data into Matlab. The text data are provided in an eight-column format without a clear georeferencing explanation. After a little tinkering, this seems to be the solution:

fid = fopen('sl.txt');
SL = textscan(fid,'%f %f %f %f %f %f %f %f','headerlines',7);
fclose(fid);

sl = [SL{1} SL{2} SL{3} SL{4} SL{5} SL{6} SL{7} SL{8}]';
sl(sl==999)=NaN;
sl = (reshape(sl,[1440 716]))';

lat = repmat(linspace(89.5,-89.5,716)',1,1440);
lon = repmat(linspace(0,360,1440),716,1);

Plot the data without the Mapping Toolbox

Now we will plot the data. To mimic the color map shown in the CU Sea Level Research Group's site, we can use the rgbmap. If you do not want to use the rgbmap function, you can simply delete the rgbmap line.

figure
pcolor(lon,lat,sl)
shading interp
cb = colorbar('location','southoutside');
caxis([-15 15])
rgbmap('dark blue','light blue','light green','yellow','magenta','white');
xlabel(cb,'sea level rise (mm/yr)')

Plot with the Mapping Toolbox

Matlab's Mapping Toolbox is decent for using georeferenced coordinates. I hear an open-source m_map package is nice too. If you have the Mapping Toolbox, you can do something similar to the above like this:

figure
ax = worldmap('world');                              % initializes a map
land = shaperead('landareas', 'UseGeoCoords', true); % reads land shapes
geoshow(ax, land, 'FaceColor', 'k')                  % shows land shapes

pcolorm(lat,lon,sl)                                  % plots sea level data
cb = colorbar('location','southoutside');
caxis([-15 15])
rgbmap('dark blue','light blue','light green','yellow','magenta','white');
xlabel(cb,'sea level rise (mm/yr)')

Interpolate to arbitrary lat/lon location(s)

Suppose you want to know the current rate of sea level rise at (11 N, 138 E)?

lati = 11;
loni = 138;

slri = interp2(lat',lon',sl',lati,loni)
slri =

         12.11

Wow, that's ~12 mm per year in the western Pacific!

Author Info

This script was written by Chad A. Greene of the Institute for Geophysics at the University of Texas at Austin. June 30, 2014.