M_Map使用手册

M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] M_Map: A mapping package for Matlab You have collected your data, loaded it into Matlab , analyzed everything to death, and now you want to make a simple map showing how it relates to the world. But you can't. Instead you have to figure out how to save all your data, and then read it into another program (like, for example GMT ), and then spend all that extra time figuring out why it doesn't give you what you expected it would...or you can invest in Matlab's own mapping toolbox (with a similarly steep learning curve)... or not! Announcing M_Map v1.4f! (updated Dec/2011) M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] M_Map is a set of mapping tools written for Matlab v5 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 either in 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 and bathymetry databases New in release 1.4f are 1. m_shaperead.m to read SHAPEFILES! (so Natural Earth can be used for political boundaries instead of DCW, among other opportunities) 2. gshhs version 2 support 3. partial OCTAVE compatibility 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 ellipsoidal earth. New in release 1.4 are 1. m_pcolor 2. m_coord (to allow for geographic and geomagnetic coordinate systems) 3. A very few minor bug fixes. 4. Some hints about and examples of adding satellite image data to your maps. 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 zip file MAKE SURE YOU ALSO UNPACK SUBDIRECTORIES! Both are M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] around 650k in size. Once you have this archive, read the Getting started section of the User's guide to correctly install this toolbox, and sections 9 and 10.1 to install TerrainBase and GSHHS respectively. A number of examples have been given to highlight the various capabilities of M_Map (thumbnails are shown above). You can also get m_namebox (a set of utilities for automatically adding names to your map), through its home page at http://www.nersc.no/~even/ . User's guide 1. Getting started 2. Specifying projections 1. Azimuthal projections 2. Cylindrical and Pseudo-cylindrical projections 3. Conic projections 4. Miscellaneous global projections 5. Yeah, but which projection should I use? 6. Map scales 7. Map coordinate systems -geographic and geomagnetic 3. Coastlines and Bathymetry 1. Coastline options 2. Topography/Bathymetry options 4. Customizing the axes 1. Grid lines and labels 2. Titles and x/ylabels 3. Legend Boxes 4. Scale bar 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 and geodesics 5. Converting longitude/latitude to projection coordinates 6. Converting projection coordinates to longitude/latitude 7. Computing distances between points 6. More complex plots 7. Removing data from a plot 8. Adding your own coastlines 1. DCW political boundaries 2. ESRI Shapefiles 9. Adding your own topography/bathymetry 1. Sandwell and Smith Bathymetry 10. Using TerrainBase 5-minute or ETOPO2 2-minute global bathymetry/topography 11. Using the GSHHS high-resolution coastline database 1. Installing GSHHS 2. Using GSHHS effectively 12. M_Map toolbox contents and description 13. Known Problems and Bugs 14. OCTAVE Compatibility Issues 15. Changes since last release For information, help, suggestions, or bug reports, contact Rich Pawlowicz , ( rich@eos.ubc.ca ) M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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 all contributed 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)); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 3. Stereographic projection 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]; M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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'); 5. Oblique Mercator Projection 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_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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'); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 7. Lambert Conformal Projection 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 M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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'); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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'); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] Examples of satellite data manipulation 1. Global SST (or any variable on a global Lat/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'); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 2. SSM/I Ice cover (data provided on a fixed grid) % SSM/I Ice concentration grids % (National Snow and Ice Data Center) P=hdfread('/mnt/cdrom/nasateam/northern/1991/feb/910222.tot','8-bit Raster Image #2'); P(P==168)=101; % Land, no coverage. P(P==157)=102; % Bad data % SSM/I ice products are in a polar stereographic projection. % and the corner points of the grid are given. Here we just % use those given corner points and 'assume' everything will % work. It's not bad, although their projection actually uses % an ellipsoidal earth (m_map uses a spherical earth). m_proj('stereographic','latitude',90,'radius',55,'rotangle',45); % Convert bottom and left corner points to screen coords. This % is of course a kludge. [MAPX,dm]=m_ll2xy([279.26 350.03],[33.92 34.35],'clip','off'); [dm,MAPY]=m_ll2xy([168.35 279.26],[30.98 33.92],'clip','off'); clf % Plot data as an image image(MAPX,MAPY,P);set(gca,'ydir','normal'); colormap([jet(100);0 0 0;1 1 1]); m_coast('patch',[.6 .6 .6]); m_grid('linewi',2,'tickdir','out'); title('SSM/I Ice cover Feb 22 1991','fontsize',14,'fontweight','bold'); h=colorbar('v'); set(get(h,'ylabel'),'string','Total Ice Concentration (%)'); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 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],[4692800 4691200],P);set(gca,'ydir','normal'); m_grid('tickdir','out','linewi',2,'fontsize',14); title('A home for certain nerds','fontsize',16); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] 4. A subset of a global dataset (HDF format) % Ocean colour data from http://seawifs.gsfc.nasa.gov/SEAWIFS.html % % Take a 4km weakly average dataset and plot a map for the Strait of % Georgia and outer coast. Note that most of this code is used % for reading in and subsetting the data. LATLIMS=[47 51]; LONLIMS=[-130 -121]; % Note - This is probably not the most efficient way to read and % handle HDF 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 map limits given P=hdfread('A20040972004104.L3m_8D_CHLO_4KM','l3m_data','Index',{[ilt ilg],[],[ltlm lglm]}); % Convert data into log(Chla) using the equations given. Blank no-data. P=double(P); P(P==255)=NaN; P=(pin.Slope*P+pin.Intercept); % log_10 of chla LT=lat(ilt+[0:ltlm-1]);LG=lon(ilg+[0:lglm-1]); [Plg,Plt]=meshgrid(LG,LT); M_Map: A Mapping package for Matlab http://www.eos.ubc.ca/~rich/map.html[2012/3/2 16:38:33] % 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 20 30]),'yticklabel',[.5 1 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 www.eos.ubc.ca M_Map: A Mapping package for Matlab