Greetings,
This weeks post will cover the basics of generating maps in Matlab. Julie’s recent post showed how to do some of this in Python, but, Matlab is also widely used by the community. You can get a lot done with Matlab, but in this post we’ll just cover a few of the basics.
We’ll start off by plotting a map of the continental United States, with the states. We used three this with three commands: usamap, shaperead, and geoshow. usamap creates an empty map axes having the Lambert Projection covering the area of the US, or any state or collection of states. shaperead reads shapefiles (duh) and returns a Matlab geographic data structure, composed of both geographic data and attributes. This Matlab data structure then interfaces really well with various Matlab functions (duh). Finally, geoshow plots geographic data, in our case on the map axes we defined. Here’s some code putting it all together.
hold on
figure1 = figure;
ax = usamap('conus');
set(ax, 'Visible', 'off')
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
states = shaperead('usastatehi',...
'UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [0.5 0.5 0.5])
tightmap
hold off
Note that ‘usastatehi’ is a shapefile containing the US states (duh) that’s distributed with Matlab. The above code generates this figure:
Now, suppose we wanted to plot some data, say a precipitation forecast, on our CONUS map. Let’s assume our forecast is being made at many points (lat,long). To interpolate between the points for plotting we’ll use Matlab’s griddata function. Once we’ve done this, we use the Matlab’s contourm command. This works exactly like the normal contour function, but the ‘m’ indicates it plots map data.
xi = min(x):0.5:max(x);
yi = min(y):0.5:max(y);
[XI, YI] = meshgrid(xi,yi);
ZI = griddata(x,y,V,XI,YI);
hold on
figure2 = figure;
ax = usamap('conus');
set(ax, 'Visible', 'off')
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
states = shaperead('usastatehi',...
'UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [0.5 0.5 0.5])
contourm(YI,-1*XI,ZI)
tightmap
hold off
Here x, y, and V are vectors of long, lat, and foretasted precipitation respectively. This code generates the following figure:
Wow! Louisiana is really getting hammered! Let’s take a closer look. We can do this by changing the entry to usamap to indicate we want to consider only Louisiana. Note, usamap accepts US postal code abbreviations.
ax = usamap('LA');
Making that change results in this figure:
Neat! We can also look at two states and add annotations. Suppose, for no reason in particular, you’re interested in the location of Tufts University relative to Cornell. We can make a map to look at this with the textm and scatterm functions. As before, the ‘m’ indicates the functions plot on a map axes.
hold on
figure4 = figure;
ax = usamap({'MA','NY'});
set(ax, 'Visible', 'off')
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
states = shaperead('usastatehi',...
'UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [0.5 0.5 0.5])
scatterm(42.4075,-71.1190,100,'k','filled')
textm(42.4075+0.2,-71.1190+0.2,'Tufts','FontSize',30)
scatterm(42.4491,-76.4842,100,'k','filled')
textm(42.4491+0.2,-76.4842+0.2,'Cornell','FontSize',30)
tightmap
hold off
This code generates the following figure.
Cool! Now back to forecasts. NOAA distributes short term Quantitative Precipitation Forecasts (QPFs) for different durations every six hours. You can download these forecasts in the form of shapefiles from a NOAA server. Here’s an example of a 24-hour rainfall forecast made at 8:22 AM UTC on April 29.
Wow, that’s a lot of rain! Can we plot our own version of this map using Matlab! You bet! Again we’ll use usamap, shaperead, and geoshow. The for loop, (0,1) scaling, and log transform are simply to make the color map more visually appealing for the post. There’s probably a cleaner way to do this, but this got the job done!
figure5 = figure;
ax = usamap('conus');
S=shaperead('94q2912','UseGeoCoords',true);
set(ax, 'Visible', 'off')
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
states = shaperead('usastatehi',...
'UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [0.5 0.5 0.5])
p = colormap(jet);
N = max(size(S));
d = zeros(N,1);
for i = 1:N
d(i) = log(S(i).QPF);
end
y=floor(((d-min(d))/range(d))*63)+1;
col = p(y,:);
for i = 1:N
geoshow(S(i),'FaceColor',col(i,:),'FaceAlpha',0.5)%,'SymbolSpec', faceColors)
end
This code generates the following figure:
If you are not plotting in the US, Matlab also has a worldmap command. This works exactly the same as usamap, but now for the world (duh). Matlab is distibuted with a shapefile ‘landareas.shp’ which contains all of the land areas in the world (duh). Generating a global map is then trivial:
figure6 = figure;
worldmap('World')
land = shaperead('landareas.shp', 'UseGeoCoords', true);
geoshow(land, 'FaceColor', [0.15 0.5 0.15])
Which generates this figure.
Matlab also comes with a number of other included that might be of interest. For instance, shapefiles detailing the locations of major world cities, lakes, and rivers. We can plot those with the following code:
figure7 = figure;
worldmap('World')
land = shaperead('landareas.shp', 'UseGeoCoords', true);
geoshow(land, 'FaceColor', [0.15 0.5 0.15])
lakes = shaperead('worldlakes', 'UseGeoCoords', true);
geoshow(lakes, 'FaceColor', 'blue')
rivers = shaperead('worldrivers', 'UseGeoCoords', true);
geoshow(rivers, 'Color', 'blue')
cities = shaperead('worldcities', 'UseGeoCoords', true);
geoshow(cities, 'Marker', '.', 'Color', 'red')
Which generates the figure:
But suppose we’re interested in one country or a group of countries. worldmap works in the same usamap does. Also, you can plot continents, for instance Europe.
worldmap('Europe')
Those are the basics, but there are many other capabilities, including 3-D projections. I can cover this in a later post if there is interest.
That’s it for now!
Curvatures 
Water Programming: A Collaborative Research Blog 