map.html

m_map是加拿大学者编写的一个在matlab上绘地图的软件

  <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>