map.html

A mapping package for Matlab:这是一款功能十分强大的地理绘图工具包

<pre>subplot(211);<br>Slongs=[-100 0;-75 25;-5 45; 25 145;45 100;145 295;100 290];<br>Slats= [&nbsp; 8 80;-80&nbsp; 8; 8 80;-80&nbsp;&nbsp; 8; 8&nbsp; 80;-80&nbsp;&nbsp; 0;&nbsp; 0&nbsp; 80];<br>for l=1:7,<br>&nbsp;m_proj('sinusoidal','long',Slongs(l,:),'lat',Slats(l,:));<br>&nbsp;m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 'ytick',[-80:20:80],'yticklabels',[],'linest','-','color',[.9 .9 .9]);<br>&nbsp;m_coast('patch','g');<br>end;<br>xlabel('Interrupted Sinusoidal Projection of World Oceans');<br>% In order to see all the maps we must undo the axis limits set by m_grid calls:<br>set(gca,'xlimmode','auto','ylimmode','auto');<br><br>subplot(212);<br>Slongs=[-100 43;-75 20; 20 145;43 100;145 295;100 295];<br>Slats= [&nbsp; 0&nbsp; 90;-90&nbsp; 0;-90&nbsp;&nbsp; 0; 0&nbsp; 90;-90&nbsp;&nbsp; 0;&nbsp; 0&nbsp; 90];<br>for l=1:6,<br>&nbsp;m_proj('mollweide','long',Slongs(l,:),'lat',Slats(l,:));<br>&nbsp;m_grid('fontsize',6,'xticklabels',[],'xtick',[-180:30:360],...<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 'ytick',[-80:20:80],'yticklabels',[],'linest','-','color','k');<br>&nbsp;m_coast('patch',[.6 .6 .6]);<br>end;<br>xlabel('Interrupted Mollweide Projection of World Oceans');<br>set(gca,'xlimmode','auto','ylimmode','auto');</pre>
<center><img src="./private/exsinus.gif"> </center>
&nbsp;
<h4><a name="5._Oblique_Mercator_Projection_with"></a> 5. Oblique
Mercator
Projection with quiver and contour data</h4>
<pre>%% Nice looking data<br>[lon,lat]=meshgrid([-136:2:-114],[36:2:54]);<br>u=sin(lat/6);<br>v=sin(lon/6);<br><br>m_proj('oblique','lat',[56 30],'lon',[-132 -120],'aspect',.8);<br><br>subplot(121);<br>m_coast('patch',[.9 .9 .9],'edgecolor','none');<br>m_grid('tickdir','out','yaxislocation','right',...<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 'xaxislocation','top','xlabeldir','end','ticklen',.02);<br>hold on;<br>m_quiver(lon,lat,u,v);<br>xlabel('Simulated surface winds');<br><br>subplot(122);<br>m_coast('patch',[.9 .9 .9],'edgecolor','none');<br>m_grid('tickdir','out','yticklabels',[],...<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 'xticklabels',[],'linestyle','none','ticklen',.02);<br>hold on;<br>[cs,h]=m_contour(lon,lat,sqrt(u.*u+v.*v));<br>clabel(cs,h,'fontsize',8);<br>xlabel('Simulated something else');</pre>
<center><img src="./private/exquiv.gif"> </center>
&nbsp;
<h4><a name="6._Miller_Projection_with_Great_Circle"></a> <a name="e6"></a>6.
Miller Projection with Great Circle</h4>
<pre>% Plot a circular orbit<br>lon=[-180:180];<br>lat=atan(tan(60*pi/180)*cos((lon-30)*pi/180))*180/pi;<br><br>m_proj('miller','lat',82);<br>m_coast('color',[0 .6 0]);<br>m_line(lon,lat,'linewi',3,'color','r');<br>m_grid('linestyle','none','box','fancy','tickdir','out');</pre>
<center><img src="./private/exmiller.gif"> </center>
&nbsp;
<h4><a name="7._Lambert_Conformal_Projection_with_Med"></a> 7. Lambert
Conformal
Projection with high-resolution bathymetry of Western Mediterranean</h4>
<pre>m_proj('lambert','lon',[-10 20],'lat',[33 48]);<br>m_tbase('contourf');<br>m_grid('linestyle','none','tickdir','out','linewidth',3);</pre>
<center><img src="./private/extbase.gif"> </center>
&nbsp;
<h4><a name="8._Demonstration_of_fancy_vectors"></a> <a name="e8"></a>8.
Demonstration of fancy vectors</h4>
<pre>m_vec&nbsp;&nbsp; % See code in m_vec.m for details</pre>
<center><img src="./private/exvec.gif"> </center>
&nbsp;
<h4><a name="9._Zoom_in_on_Prince_Edward_Island_to_co"></a> <a
 name="e9"></a>9.
Zoom in on Prince Edward Island to show different coastline resolutions</h4>
<pre>% Example showing the default coastline and all of the different resolutions&nbsp;<br>% of GSHHS coastlines as we zoom in on a section of Prince Edward Island.<br><br>clf<br>axes('position',[.35 .6 .37 .37]);<br>m_proj('albers equal-area','lat',[40 60],'long',[-90 -50],'rect','on');<br>m_coast('patch',[0 1 0]);<br>m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right');<br>m_text(-69,41,'Standard coastline','color','r','fontweight','bold');<br><br>axes('position',[.09 .5 .37 .37]);<br>m_proj('albers equal-area','lat',[40 54],'long',[-80 -55],'rect','on');<br>m_gshhs_c('patch',[.2 .8 .2]);<br>m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top');<br>m_text(-80,52.5,'GSHHS\_C (crude)','color','m','fontweight','bold','fontsize',14);<br><br>axes('position',[.13 .2 .37 .37]);<br>m_proj('albers equal-area','lat',[43 48],'long',[-67 -59],'rect','on');<br>m_gshhs_l('patch',[.4 .6 .4]);<br>m_grid('linest','none','linewidth',2,'tickdir','out');<br>m_text(-66.5,43.5,'GSHHS\_L (low)','color','m','fontweight','bold','fontsize',14);<br><br>axes('position',[.35 .05 .37 .37]);<br>m_proj('albers equal-area','lat',[45.8 47.2],'long',[-64.5 -62],'rect','on');<br>m_gshhs_i('patch',[.5 .6 .5]);<br>m_grid('linest','none','linewidth',2,'tickdir','out','yaxisloc','right');<br>m_text(-64.4,45.9,'GSHHS\_I (intermediate)','color','m','fontweight','bold','fontsize',14);<br><br>axes('position',[.55 .23 .37 .37]);<br>m_proj('albers equal-area','lat',[46.375 46.6],'long',[-64.2 -63.7],'rect','on');<br>m_gshhs_h('patch',[.6 .6 .6]);<br>m_grid('linest','none','linewidth',2,'tickdir','out','xaxisloc','top','yaxisloc','right');<br>m_text(-64.18,46.58,'GSHHS\_H (high)','color','m','fontweight','bold','fontsize',14);</pre>
<center><img src="./private/exgshhs.gif"> </center>
&nbsp;
<h4><a name="10._Tracklines_and_UTM_projection"></a> <a name="e10"></a>10.
Tracklines and UTM projection</h4>
<pre>m_proj('UTM','long',[-72 -68],'lat',[40 44]);<br>m_gshhs_i('color','k');<br>m_grid('box','fancy','tickdir','in');<br><br>% fake up a trackline<br>lons=[-71:.1:-67];<br>lats=60*cos((lons+115)*pi/180);<br>dates=datenum(1997,10,23,15,1:41,zeros(1,41));<br><br>m_track(lons,lats,dates,'ticks',0,'times',4,'dates',8,...<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 'clip','off','color','r','orient','upright');</pre>
<center><img src="./private/track1.gif"> </center>
&nbsp;
<h4><a name="11._Range_rings"></a> <a name="e11"></a>11. Range rings</h4>
<pre>&nbsp;&nbsp;&nbsp; m_proj('hammer','clong',170);<br>&nbsp;&nbsp;&nbsp; m_grid('xtick',[],'ytick',[],'linestyle','-');<br>&nbsp;&nbsp;&nbsp; m_coast('patch','g');<br>&nbsp;&nbsp;&nbsp; m_line(100.5,13.5,'marker','square','color','r');<br>&nbsp;&nbsp;&nbsp; m_range_ring(100.5,13.5,[1000:1000:15000],'color','b','linewi',2);<br>&nbsp;&nbsp;&nbsp; xlabel('1000km range rings from Bangkok');</pre>
<center><img src="./private/rring.gif"><br>
<br>
</center>
<h4><a name="12._Speckle"></a> <a name="e12"></a>12. Speckled boundary</h4>
<pre>
    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&amp;W presentation (best in postscript format)','fontsize',14);
    m_text(-128,48,5,{'Pacific','Ocean'},'fontsize',18);
</pre>
<center><img src="./private/extspeckle.gif"><br>
<br>
</center>
<h4><a name="13._Blueocean"></a> <a name="e13"></a>13. Blue Ocean</h4>
<pre>    m_proj('miller','lat',[-75 75]);<br>    <br>    set(gca,'color',[.9 .99 1]);     % Trick is to set this *before* the patch call.<br>    <br>    m_coast('patch',[.7 1 .7],'edgecolor','none');<br>    m_grid('box','fancy','linestyle','none');<br>       <br>    cities={'Cairo','Washington','Buenos Aires'};<br>    lons=[ 30+2/60  -77-2/60   -58-22/60];<br>    lats=[ 31+21/60  38+53/60  -34-45/60];<br>    <br>    for k=1:3,<br>      [range,ln,lt]=m_lldist([-123-6/60 lons(k)],[49+13/60  lats(k)],40);<br>      m_line(ln,lt,'color','r','linewi',2);<br>      m_text(ln(end),lt(end),sprintf('%s - %d km',cities{k},round(range)));<br>    end;<br>      <br>    title('Great Circle Routes','fontsize',14,'fontweight','bold');<br></pre>
<center><img src="./private/extblueocean.gif"><br>
<br>
</center>
<h2><a name="satellite_examples"></a>Examples of satellite data
manipulation</h2>
<h4><a name="1._Global_SST"></a>1. Global SST (or any variable on a
global
Lat/Long grid)</h4>
<blockquote><tt>% NOAA/NASA Pathfinder AVHRR SST product</tt><br>
  <tt>% http://podaac.jpl.nasa.gov/sst/</tt><br>
  <br>
  <tt>[P,map]=imread('../m_mapWK/199911h54ma-gdm.hdf');</tt><br>
  <br>
  <tt>% Documentation for the 54km dataset gives</tt><br>
  <tt>% this formula for temperature</tt><br>
  <tt>P=0.15*double(P)-3; % deg C</tt><br>
  <br>
  <tt>%...and defines this Lat/Long grid for the data</tt><br>
  <tt>Plat=90-.25-[0:359]*.5;Plon=-180+.25+[0:719]*.5;</tt><br>
  <br>
  <tt>% Since the grid is rectangluar in lat/long (i.e. not</tt><br>
  <tt>% really a projection at all, althouhg it is included in</tt><br>
  <tt>% m_map under the name 'equidistant cyldindrical'), we </tt><br>
  <tt>% don't want to use the 'image' technique. Instead... </tt><br>
  <tt>% Create a grid, offsetting by half a grid point to account</tt><br>
  <tt>% for the flat pcolor</tt><br>
  <tt>[Plg,Plt]=meshgrid(Plon-0.25,Plat+0.25);</tt><br>
  <br>
  <tt>m_proj('hammer-aitoff','clongitude',-150);</tt><br>
  <br>
  <tt>% Rather than rearranging the data so its limits match the</tt><br>
  <tt>% plot I just draw it twice (you can see the join at 180W</tt><br>
  <tt>% because of the quirks of flat pcolor) (Note that </tt><br>
  <tt>% all the global projections have 360 deg ambiguities)</tt><br>
  <tt>m_pcolor(Plg,Plt,P);shading flat;colormap(map);</tt><br>
  <tt>hold on;</tt><br>
  <tt>m_pcolor(Plg-360,Plt,P);shading flat;colormap(map);</tt><br>
  <br>
  <tt>m_coast('patch',[.6 1 .6]);</tt><br>
  <tt>m_grid('xaxis','middle');</tt><br>
  <br>
  <tt>% add a standard colorbar.</tt><br>
  <tt>h=colorbar('h');</tt><br>
  <tt>set(get(h,'title'),'string','AVHRR SST Nov 1999');</tt><br>
</blockquote>
<pre> </pre>
<div align="center"><img src="private/ex_sst.gif" alt="SST pic"
 height="376" width="553">
<br>
<br>
<div align="left"><br>
<h4><a name="2._SSMI_Ice_cover"></a>2. &nbsp;SSM/I Ice cover (data
provided
on a fixed grid)</h4>
<blockquote><tt>% SSM/I Ice concentration grids</tt><br>
  <tt>% (National Snow and Ice Data Center)</tt><br>
  <br>
  <tt>P=hdfread('/mnt/cdrom/nasateam/northern/1991/feb/910222.tot','8-bit
Raster Image #2');</tt><br>
  <tt>P(P==168)=101; % Land, no coverage.</tt><br>
  <tt>P(P==157)=102; % Bad data</tt><br>
  <br>
  <tt>% SSM/I ice products are in a polar stereographic projection.</tt><br>
  <tt>% and the corner points of the grid are given. Here we just</tt><br>
  <tt>% use those given corner points and 'assume' everything will</tt><br>
  <tt>% work. It's not bad, although their projection actually uses</tt><br>
  <tt>% an ellipsoidal earth (m_map uses a spherical earth).</tt><br>
  <br>
  <tt>m_proj('stereographic','latitude',90,'radius',55,'rotangle',45);</tt><br>
  <br>
  <tt>% Convert bottom and left corner points to screen coords. This</tt><br>
  <tt>% is of course a kludge.</tt><br>
  <tt>[MAPX,dm]=m_ll2xy([279.26 350.03],[33.92 34.35],'clip','off');</tt><br>
  <tt>[dm,MAPY]=m_ll2xy([168.35 279.26],[30.98 33.92],'clip','off');</tt><br>
  <br>
  <tt>clf</tt><br>
  <tt>% Plot data as an image</tt><br>
  <tt>image(MAPX,MAPY,P);set(gca,'ydir','normal');</tt><br>
  <tt>colormap([jet(100);0 0 0;1 1 1]);</tt><br>
  <br>
  <tt>m_coast('patch',[.6 .6 .6]);</tt><br>
  <tt>m_grid('linewi',2,'tickdir','out');</tt><br>
  <tt>title('SSM/I Ice cover Feb 22
1991','fontsize',14,'fontweight','bold');</tt><br>
  <br>
  <tt>h=colorbar('v');</tt><br>
  <tt>set(get(h,'ylabel'),'string','Total Ice Concentration (%)');</tt><br>
</blockquote>
<br>
<div align="center"><img src="private/ex_ssmi.gif" alt="SSMi Ice"
 height="413" width="522">
<br>
</div>
</div>
</div>
<p> </p>
<br>
<br>
<br>
<h4><a name="3._Aerial_photos"></a>3. &nbsp;Aerial photos on an UTM grid</h4>
<blockquote><tt>% This image comes from the TerraServer</tt><br>
  <tt>% (http://terraserver.microsoft.com/)</tt><br>
  <tt>% and has been georeferenced to UTM coords. The UTM projection</tt><br>
  <tt>% uses UTM coordinates on the screen (as long as the ellipse</tt><br>
  <tt>% parameter is set to something other than the default).</tt><br>
  <tt>[P,map]=imread('../m_mapWK/oncehome.jpeg');</tt><br>
  <br>
  <tt>% Set the projection limits to the lat/long of image</tt><br>
  <tt>% corners.</tt><br>
  <tt>m_proj('UTM','long',[-71-6/60-30/3600 -71-4/60-43/3600],...</tt><br>
  <tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
'lat',[42+21/60+13/3600&nbsp; 42+22/60+7/3600],'ellipse','wgs84');</tt><br>
  <br>
  <tt>clf;</tt><br>
  <tt>image([326400 328800],[4692800
4691200],P);set(gca,'ydir','normal');</tt><br>
  <tt>m_grid('tickdir','out','linewi',2,'fontsize',14);</tt><br>
  <tt>title('A home for certain nerds','fontsize',16);</tt><br>
</blockquote>
<div align="center"><img src="private/ex_terra.gif" alt="nerdhome"
 height="401" width="542">
<br>
</div>
<br>
<br>
<h4><a name="4._A_subset_of_a_global_dataset"></a>4. &nbsp;A subset of
a
global dataset (HDF format)</h4>
<blockquote><font size="-1"><tt>% Ocean colour data from
http://seawifs.gsfc.nasa.gov/SEAWIFS.html</tt></font><br>
  <font size="-1"><tt>%</tt></font><br>
  <font size="-1"><tt>% Take a 4km weakly average dataset and plot a
map
for the Strait of</tt></font><br>
  <font size="-1"><tt>% Georgia and outer coast. Note that most of this
code
is used</tt></font><br>
  <font size="-1"><tt>% for reading in and subsetting the data.</tt></font><br>
  <br>
  <font size="-1"><tt>LATLIMS=[47 51];</tt></font><br>
  <font size="-1"><tt>LONLIMS=[-130 -121];</tt></font><br>
  <br>
  <font size="-1"><tt>% Note - This is probably not the most efficient
way
to read and</tt></font><br>
  <font size="-1"><tt>%&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
handle
HDF data, but I don't usually do this...</tt></font><br>
  <font size="-1"><tt>%</tt></font><br>
  <font size="-1"><tt>% First, get the attribute data</tt></font><br>
  <font size="-1"><tt>PI=hdfinfo('A20040972004104.L3m_8D_CHLO_4KM');</tt></font><br>
  <font size="-1"><tt>% And write it into a structure</tt></font><br>
  <font size="-1"><tt>pin=[];</tt></font><br>
  <font size="-1"><tt>for k=1:59,</tt></font><br>
  <font size="-1"><tt>&nbsp; nm=PI.Attributes(k).Name;nm(nm==' ')='_';</tt></font><br>
  <font size="-1"><tt>&nbsp; if isstr(PI.Attributes(k).Value),</tt></font><br>
  <font size="-1"><tt>&nbsp;&nbsp;&nbsp;
pin=setfield(pin,nm,PI.Attributes(k).Value);</tt></font><br>
  <font size="-1"><tt>&nbsp; else</tt></font><br>
  <font size="-1"><tt>&nbsp;&nbsp;&nbsp;
pin=setfield(pin,nm,double(PI.Attributes(k).Value));</tt></font><br>
  <font size="-1"><tt>&nbsp; end</tt></font><br>
  <font size="-1"><tt>end;&nbsp; </tt></font><br>
  <br>
  <font size="-1"><tt>% lon/lat of grid corners</tt></font><br>
  <font size="-1"><tt>lon=[pin.Westernmost_Longitude:pin.Longitude_Step:pin.Easternmost_Longitude];</tt></font><br>
  <font size="-1"><tt>lat=[pin.Northernmost_Latitude:-pin.Latitude_Step:pin.Southernmost_Latitude];</tt></font><br>
  <br>
  <font size="-1"><tt>% Get the indices needed for the area of interest</tt></font><br>
  <font size="-1"><tt>[mn,ilt]=min(abs(lat-max(LATLIMS)));</tt></font><br>
  <font size="-1"><tt>[mn,ilg]=min(abs(lon-min(LONLIMS)));</tt></font><br>
  <font size="-1"><tt>ltlm=fix(diff(LATLIMS)/pin.Latitude_Step);</tt></font><br>
  <font size="-1"><tt>lglm=fix(diff(LONLIMS)/pin.Longitude_Step);</tt></font><br>
  <br>
  <font size="-1"><tt>% load the subset of data needed for the map
limits
given</tt></font><br>
  <font size="-1"><tt>P=hdfread('A20040972004104.L3m_8D_CHLO_4KM','l3m_data','Index',{[ilt
ilg],[],[ltlm lglm]});</tt></font><br>
  <br>
  <font size="-1"><tt>% Convert data into log(Chla) using the equations
given.
Blank no-data.</tt></font><br>
  <font size="-1"><tt>P=double(P);</tt></font><br>
  <font size="-1"><tt>P(P==255)=NaN;</tt></font><br>
  <font size="-1"><tt>P=(pin.Slope*P+pin.Intercept);&nbsp;&nbsp; %
log_10
of chla</tt></font><br>
  <br>
  <font size="-1"><tt>LT=lat(ilt+[0:ltlm-1]);LG=lon(ilg+[0:lglm-1]);</tt></font><br>
  <font size="-1"><tt>[Plg,Plt]=meshgrid(LG,LT);</tt></font><br>
  <br>
  <font size="-1"><tt>% Draw the map...</tt></font><br>
  <br>
  <font size="-1"><tt>clf;</tt></font><br>
  <font size="-1"><tt>m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);</tt></font><br>
  <font size="-1"><tt>m_pcolor(Plg,Plt,P);shading flat;</tt></font><br>
  <font size="-1"><tt>m_gshhs_i('color','k');;</tt></font><br>
  <font size="-1"><tt>m_grid('linewi',2,'tickdir','out');;</tt></font><br>
  <font size="-1"><tt>h=colorbar;</tt></font><br>
  <font size="-1"><tt>set(get(h,'ylabel'),'String','Chla (\mug/l)');</tt></font><br>
  <font size="-1"><tt>set(h,'ytick',log10([.5 1 2 3 5 10 20
30]),'yticklabel',[.5
1 2 3 5 10 20 30],'tickdir','out');</tt></font><br>
  <font size="-1"><tt>title(['MODIS Chla&nbsp;&nbsp; '
datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_Start_Day)
' -&gt; ' ...</tt></font><br>
  <font size="-1"><tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
datestr(datenum(pin.Period_Start_Year,1,0)+pin.Period_End_Day)],...</tt></font><br>
  <font size="-1"><tt>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
'fontsize',14,'fontweight','bold');</tt></font><br>
</blockquote>
<br>
<div align="center"><img src="private/ex_modis.gif" alt="modis"
 height="396" width="540">
<br>
</div>
<hr><i>Last changed 30/Dec/2005. Questions and comments to <a
 href="mailto:[e-mail]">[e-mail]</a></i> <br>
<br>
</body>
</html>