A Mapping package for Matlab (2024)

You have collected your data, loaded it into Matlab , analyzed everything todeath, and now you want to make a simple map showing how it relates tothe world.

But you can't.

Instead you have to figure out how to save all your data, and thenread it into another program (like, for example GMT ), and then spend all thatextratime figuring out why it doesn't give you what you expected itwould...oryou can invest in Matlab's own mapping toolbox (with a similarly steeplearningcurve)... or not!

Announcing M_Map v1.4e! (updated Oct/09)

M_Map is a set of mapping tools written for Matlabv5 and later. These include:

1. Routines to project data in 19 different spherical projections(and determine inverse mappings)
2. A grid generation routine to make nice axes with limits eitherin lat/long terms or in planar X/Y terms.
3. A coastline database (with 1/4 degree resolution)
4. A global elevation database (1 degree resolution)
5. Hooks into freely available high-resolution coastline andbathymetry databases

New in release 1.4e are

1. GSHHS Now comes with the WDB rivers and borders database, and links to them have been added via m_gshhs.m
2. Ellipsoidal Albers and Lambert conic projections included

New in release 1.4d are

1. Modifications to work around bugs in matlab7 contourf
2. Robinson projection
3. A few compatibility issues with current matlab versions
4. Upgraded hooks into some databases.

New in release 1.4b are

1. m_hatch for hatched and speckled patches
2. old-fashioned speckled coastlines (good for B&W pics - see Example 13).
   
3. m_lldist now also returns points on great circle geodesics.
4. m_fdist, m_idist, and m_geodesic for geodesics on an ellipsoidalearth.
   

New in release 1.4 are

1. m_pcolor
2. m_coord (to allow for geographic and geomagnetic coordinatesystems)
3. A very few minor bug fixes.
4. Some hints about and examples of adding satellite image data toyourmaps.
   

How to get M_Map

You can download the M_Map toolbox either as a gzipped tar-file ,or as zip archive(Click on these links to download). If you are unpacking the zipfile MAKE SURE YOUALSO UNPACK SUBDIRECTORIES! Both are around 650k in size. Once you havethisarchive, read the Getting started sectionof the User's guide to correctlyinstall this toolbox, and sections 9and 10.1 to install TerrainBaseand GSHHS respectively.

A number of examples have been given tohighlight the various capabilities of M_Map (thumbnails are shownabove).

You can also get m_namebox (a set of utilities forautomatically adding names to your map), through its home page athttp://www.nersc.no/~even/.

User's guide

1. Getting started
2. Specifying projections
1. Azimuthal projections
2. Cylindrical andPseudo-cylindrical projections
3. Conic projections
4. Miscellaneous globalprojections
5. Yeah, but which projectionshould I use?
6. Map scales
7. Map coordinate systems-geographic and geomagnetic
   
3. Coastlines and Bathymetry
1. Coastline options
2. Topography/Bathymetryoptions
4. Customizing the axes
1. Grid lines and labels
2. Titles and x/ylabels
3. Legend Boxes
5. Adding your own data
1. Drawing lines, text,arrows, patches, hatches, speckles, and contours
2. Drawing images and p_color
3. Drawing tracklines
4. Drawing range rings andgeodesics
5. Convertinglongitude/latitude to projection coordinates
6. Converting projectioncoordinates to longitude/latitude
7. Computing distances betweenpoints
6. More complex plots
7. Removing data from a plot
8. Adding your own coastlines
1. DCW political boundaries
9. Adding your owntopography/bathymetry
1. Sandwell and SmithBathymetry
10. Using TerrainBase 5-minuteor ETOPO2 2-minute global bathymetry/topography
11. Using the GSHHShigh-resolution coastline database
1. Installing GSHHS
2. Using GSHHS effectively
12. M_Map toolbox contents anddescription
13. Known Problems and Bugs
14. Changes since lastrelease

For information, help, suggestions, or bug reports, contact Rich Pawlowicz ,( rich@eos.ubc.ca )

Acknowledgements

A number of people have helped out with suggestions, code fixes,etc. I am especially grateful for the work done by E. Firing, D.Byrne, M. Mann, J. Pringle, and J. E. Nilsen who have allcontributed code.

Examples

1. M_Map Logo

m_proj('ortho','lat',48','long',-123');
m_coast('patch','r');
m_grid('linest','-','xticklabels',[],'yticklabels',[]);
patch(.55*[-1 1 1 -1],.25*[-1 -1 1 1]-.55,'w');
text(0,-.55,'M\_Map','fontsize',25,'color','b',...
 'vertical','middle','horizontal','center');
set(gcf,'units','inches','position',[2 2 3 3]);
set(gcf,'paperposition',[3 3 3 3]);

2.Lambert Conformal Conic projection of North American Topography

m_proj('lambert','long',[-160 -40],'lat',[30 80]);
m_coast('patch',[1 .85 .7]);
m_elev('contourf',[500:500:6000]);
m_grid('box','fancy','tickdir','in');
colormap(flipud(copper));

3.Stereographicprojection of North Polar regions

% Note that coastline is drawn OVER the grid because of the order in which
% the two routines are called
m_proj('stereographic','lat',90,'long',30,'radius',25);
m_elev('contour',[-3500:1000:-500],'edgecolor','b');
m_grid('xtick',12,'tickdir','out','ytick',[70 80],'linest','-');
m_coast('patch',[.7 .7 .7],'edgecolor','r');

4.Two Interrupted Projections of the World's Oceans

subplot(211);
Slongs=[-100 0;-75 25;-5 45; 25 145;45 100;145 295;100 290];
Slats= [ 8 80;-80 8; 8 80;-80 8; 8 80;-80 0; 0 80];
for l=1:7,
m_proj('sinusoidal','long',Slongs(l,:),'lat',Slats(l,:));
m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...
 'ytick',[-80:20:80],'yticklabels',[],'linest','-','color',[.9 .9 .9]);
m_coast('patch','g');
end;
xlabel('Interrupted Sinusoidal Projection of World Oceans');
% In order to see all the maps we must undo the axis limits set by m_grid calls:
set(gca,'xlimmode','auto','ylimmode','auto');
subplot(212);
Slongs=[-100 43;-75 20; 20 145;43 100;145 295;100 295];
Slats= [ 0 90;-90 0;-90 0; 0 90;-90 0; 0 90];
for l=1:6,
m_proj('mollweide','long',Slongs(l,:),'lat',Slats(l,:));
m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...
 'ytick',[-80:20:80],'yticklabels',[],'linest','-','color','k');
m_coast('patch',[.6 .6 .6]);
end;
xlabel('Interrupted Mollweide Projection of World Oceans');
set(gca,'xlimmode','auto','ylimmode','auto');
See Also
Write a Map Function
- MATLAB & Simulink
- MathWorks América LatinaObject that maps unique keys to valuesCreate axesm-based map for world regionWrite a Map Function
- MATLAB & Simulink
- MathWorks France

5. ObliqueMercatorProjection with quiver and contour data

%% Nice looking data
[lon,lat]=meshgrid([-136:2:-114],[36:2:54]);
u=sin(lat/6);
v=sin(lon/6);
m_proj('oblique','lat',[56 30],'lon',[-132 -120],'aspect',.8);
subplot(121);
m_coast('patch',[.9 .9 .9],'edgecolor','none');
m_grid('tickdir','out','yaxislocation','right',...
 'xaxislocation','top','xlabeldir','end','ticklen',.02);
hold on;
m_quiver(lon,lat,u,v);
xlabel('Simulated surface winds');
subplot(122);
m_coast('patch',[.9 .9 .9],'edgecolor','none');
m_grid('tickdir','out','yticklabels',[],...
 'xticklabels',[],'linestyle','none','ticklen',.02);
hold on;
[cs,h]=m_contour(lon,lat,sqrt(u.*u+v.*v));
clabel(cs,h,'fontsize',8);
xlabel('Simulated something else');

6.Miller Projection with Great Circle

% Plot a circular orbit
lon=[-180:180];
lat=atan(tan(60*pi/180)*cos((lon-30)*pi/180))*180/pi;
m_proj('miller','lat',82);
m_coast('color',[0 .6 0]);
m_line(lon,lat,'linewi',3,'color','r');
m_grid('linestyle','none','box','fancy','tickdir','out');

7. LambertConformalProjection with high-resolution bathymetry of Western Mediterranean

m_proj('lambert','lon',[-10 20],'lat',[33 48]);
m_tbase('contourf');
m_grid('linestyle','none','tickdir','out','linewidth',3);

8.Demonstration of fancy vectors

m_vec % See code in m_vec.m for details

9.Zoom in on Prince Edward Island to show different coastline resolutions

% Example showing the default coastline and all of the different resolutions
% of GSHHS coastlines as we zoom in on a section of Prince Edward Island.
clf
axes('position',[.35 .6 .37 .37]);
m_proj('albers equal-area','lat',[40 60],'long',[-90 -50],'rect','on');
m_coast('patch',[0 1 0]);
m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right');
m_text(-69,41,'Standard coastline','color','r','fontweight','bold');
axes('position',[.09 .5 .37 .37]);
m_proj('albers equal-area','lat',[40 54],'long',[-80 -55],'rect','on');
m_gshhs_c('patch',[.2 .8 .2]);
m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top');
m_text(-80,52.5,'GSHHS\_C (crude)','color','m','fontweight','bold','fontsize',14);
axes('position',[.13 .2 .37 .37]);
m_proj('albers equal-area','lat',[43 48],'long',[-67 -59],'rect','on');
m_gshhs_l('patch',[.4 .6 .4]);
m_grid('linest','none','linewidth',2,'tickdir','out');
m_text(-66.5,43.5,'GSHHS\_L (low)','color','m','fontweight','bold','fontsize',14);
See Also
A Mapping package for Matlab
axes('position',[.35 .05 .37 .37]);
m_proj('albers equal-area','lat',[45.8 47.2],'long',[-64.5 -62],'rect','on');
m_gshhs_i('patch',[.5 .6 .5]);
m_grid('linest','none','linewidth',2,'tickdir','out','yaxisloc','right');
m_text(-64.4,45.9,'GSHHS\_I (intermediate)','color','m','fontweight','bold','fontsize',14);
axes('position',[.55 .23 .37 .37]);
m_proj('albers equal-area','lat',[46.375 46.6],'long',[-64.2 -63.7],'rect','on');
m_gshhs_h('patch',[.6 .6 .6]);
m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right');
m_text(-64.18,46.58,'GSHHS\_H (high)','color','m','fontweight','bold','fontsize',14);

10.Tracklines and UTM projection

m_proj('UTM','long',[-72 -68],'lat',[40 44]);
m_gshhs_i('color','k');
m_grid('box','fancy','tickdir','in');
% fake up a trackline
lons=[-71:.1:-67];
lats=60*cos((lons+115)*pi/180);
dates=datenum(1997,10,23,15,1:41,zeros(1,41));
m_track(lons,lats,dates,'ticks',0,'times',4,'dates',8,...
 'clip','off','color','r','orient','upright');

11. Range rings

 m_proj('hammer','clong',170);
 m_grid('xtick',[],'ytick',[],'linestyle','-');
 m_coast('patch','g');
 m_line(100.5,13.5,'marker','square','color','r');
 m_range_ring(100.5,13.5,[1000:1000:15000],'color','b','linewi',2);
 xlabel('1000km range rings from Bangkok');

12. Speckled boundary

 bndry_lon=[-128.8 -128.8 -128.3 -128 -126.8 -126.6 -128.8]; bndry_lat=[49 50.33 50.33 50 49.5 49 49]; clf; m_proj('lambert','long',[-130 -121.5],'lat',[47 51.5],'rectbox','on'); m_gshhs_i('color','k'); % Coastline... m_gshhs_i('speckle','color','k'); % with speckle added m_line(bndry_lon,bndry_lat,'linewi',2,'color','k'); % Area outline ... m_hatch(bndry_lon,bndry_lat,'single',30,5,'color','k'); % ...with hatching added. m_grid('linewi',2,'linest','none','tickdir','out','fontsize',12); title('Speckled Boundaries for nice B&W presentation (best in postscript format)','fontsize',14); m_text(-128,48,5,{'Pacific','Ocean'},'fontsize',18);

13. Blue Ocean

 m_proj('miller','lat',[-75 75]);
 set(gca,'color',[.9 .99 1]); % Trick is to set this *before* the patch call.
 m_coast('patch',[.7 1 .7],'edgecolor','none');
 m_grid('box','fancy','linestyle','none');
 cities={'Cairo','Washington','Buenos Aires'};
 lons=[ 30+2/60 -77-2/60 -58-22/60];
 lats=[ 31+21/60 38+53/60 -34-45/60];
 for k=1:3,
 [range,ln,lt]=m_lldist([-123-6/60 lons(k)],[49+13/60 lats(k)],40);
 m_line(ln,lt,'color','r','linewi',2);
 m_text(ln(end),lt(end),sprintf('%s - %d km',cities{k},round(range)));
 end;
 title('Great Circle Routes','fontsize',14,'fontweight','bold');

Examples of satellite datamanipulation

1. Global SST (or any variable on aglobalLat/Long grid)

% NOAA/NASA Pathfinder AVHRR SST product
% http://podaac.jpl.nasa.gov/sst/[P,map]=imread('../m_mapWK/199911h54ma-gdm.hdf');% Documentation for the 54km dataset gives
% this formula for temperature
P=0.15*double(P)-3; % deg C%...and defines this Lat/Long grid for the data
Plat=90-.25-[0:359]*.5;Plon=-180+.25+[0:719]*.5;% Since the grid is rectangluar in lat/long (i.e. not
% really a projection at all, althouhg it is included in
% m_map under the name 'equidistant cyldindrical'), we
% don't want to use the 'image' technique. Instead...
% Create a grid, offsetting by half a grid point to account
% for the flat pcolor
[Plg,Plt]=meshgrid(Plon-0.25,Plat+0.25);m_proj('hammer-aitoff','clongitude',-150);% Rather than rearranging the data so its limits match the
% plot I just draw it twice (you can see the join at 180W
% because of the quirks of flat pcolor) (Note that
% all the global projections have 360 deg ambiguities)
m_pcolor(Plg,Plt,P);shading flat;colormap(map);
hold on;
m_pcolor(Plg-360,Plt,P);shading flat;colormap(map);m_coast('patch',[.6 1 .6]);
m_grid('xaxis','middle');% add a standard colorbar.
h=colorbar('h');
set(get(h,'title'),'string','AVHRR SST Nov 1999');

3. Aerial photos on an UTM grid

% This image comes from the TerraServer
% (http://terraserver.microsoft.com/)
% and has been georeferenced to UTM coords. The UTM projection
% uses UTM coordinates on the screen (as long as the ellipse
% parameter is set to something other than the default).
[P,map]=imread('../m_mapWK/oncehome.jpeg');% Set the projection limits to the lat/long of image
% corners.
m_proj('UTM','long',[-71-6/60-30/3600 -71-4/60-43/3600],...
'lat',[42+21/60+13/3600 42+22/60+7/3600],'ellipse','wgs84');clf;
image([326400 328800],[46928004691200],P);set(gca,'ydir','normal');
m_grid('tickdir','out','linewi',2,'fontsize',14);
title('A home for certain nerds','fontsize',16);

4. A subset ofa*global dataset (HDF format)

% Ocean colour data fromhttp://seawifs.gsfc.nasa.gov/SEAWIFS.html
%
% Take a 4km weakly average dataset and plot amapfor the Strait of
% Georgia and outer coast. Note that most of thiscodeis used
% for reading in and subsetting the data.LATLIMS=[47 51];
LONLIMS=[-130 -121];% Note - This is probably not the most efficientwayto read and
%handleHDF data, but I don't usually do this...
%
% First, get the attribute data
PI=hdfinfo('A20040972004104.L3m_8D_CHLO_4KM');
% And write it into a structure
pin=[];
for k=1:59,
nm=PI.Attributes(k).Name;nm(nm==' ')='_';
if isstr(PI.Attributes(k).Value),
pin=setfield(pin,nm,PI.Attributes(k).Value);
else
pin=setfield(pin,nm,double(PI.Attributes(k).Value));
end
end; % lon/lat of grid corners
lon=[pin.Westernmost_Longitude:pin.Longitude_Step:pin.Easternmost_Longitude];
lat=[pin.Northernmost_Latitude:-pin.Latitude_Step:pin.Southernmost_Latitude];% Get the indices needed for the area of interest
[mn,ilt]=min(abs(lat-max(LATLIMS)));
[mn,ilg]=min(abs(lon-min(LONLIMS)));
ltlm=fix(diff(LATLIMS)/pin.Latitude_Step);
lglm=fix(diff(LONLIMS)/pin.Longitude_Step);% load the subset of data needed for the maplimitsgiven
P=hdfread('A20040972004104.L3m_8D_CHLO_4KM','l3m_data','Index',{[iltilg],[],[ltlm lglm]});% Convert data into log(Chla) using the equationsgiven.Blank no-data.
P=double(P);
P(P==255)=NaN;
P=(pin.Slope*P+pin.Intercept); %log_10of chlaLT=lat(ilt+[0:ltlm-1]);LG=lon(ilg+[0:lglm-1]);
[Plg,Plt]=meshgrid(LG,LT);% Draw the map...clf;
m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);
m_pcolor(Plg,Plt,P);shading flat;
m_gshhs_i('color','k');;
m_grid('linewi',2,'tickdir','out');;
h=colorbar;
set(get(h,'ylabel'),'String','Chla (\mug/l)');
set(h,'ytick',log10([.5 1 2 3 5 10 2030]),'yticklabel',[.51 2 3 5 10 20 30],'tickdir','out');
title(['MODIS Chla 'datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_Start_Day)' -> ' ...
datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_End_Day)],...
'fontsize',14,'fontweight','bold');

Last changed 30/Dec/2005. Questions and comments to rich@eos.ubc.ca