matlab中地图边界与掩膜(去掉边界外区域)的实现（基于shape文件）

斥鷃

多谢@我和小云 后面的补充以及版主们的编辑，大家可以看看一楼的推荐帖，的确比我的方法简单很多，在此推荐！另裂墙推荐——有关无锯齿白化问题可以参考@Lighting 的MATLAB对地图进行白化
论坛有很多人问到过这个问题，解决这个问题也有助于matlab在地学中的应用，前天LZ花了几个小时看help map的文字（英语要学好啊！！！），粗略地研究了一下，总算找到了几个有用的函数，现在把函数和做法分享给大家，以下是一张效果图（网格大小的关系，还是有些地方空了多了~）：

首先是论坛里的一些资源，其实绘制地理类的图形，我一般都推荐用grads、ArcGis（不是很熟悉就不提NCL、surfer和MeteoInfo等等啦，论坛里关于这几个软件的帖子也很多，不一一举例了）~有关这两个软件底图制作和绘制在论坛里面有很详细的教程，大家可以参考（grads中maskout方法，ArcGis底图制作，ArcGis站点绘图参考），而有关matlab底图方法，我估计m_map工具箱中会有一些解决方法，但是也不是很熟悉，我就用matlab自带的map工具箱中的函数进行处理了（不清楚各版本有没有装，改天看到工具箱下载的时候我共享一下）。

一、预备知识：
1.地理参考（R）：地理参考（具体内容可在help大部分map工具箱中的函数得到）是指定位矩阵用的矩阵，通常为3*2的矩阵，通过公式
                          [lon lat] = [row col 1] * R（3*2）    (具体实例见第三部分中程序框)
对矩阵的行列坐标与经纬度进行转换（也有1*3的形式，有兴趣的朋友可以自己了解一下），是map工具箱处理地理数据的基础；
2.矢量数据：一般矢量数据分为多边形、线和点数据，我们这里采用的shape文件的多边形数据是利用线数据定义边界的，在matlab里shape文件以结构数组的形式保存，其中经纬度信息保存在X与Y属性中；
3.2-D矩阵数据：由行标(row)和列标(col)来表示位置的数据，需要通过坐标系来对数据进行对应进行绘图，常用的绘图参考方式有meshgrid函数生成的参考系和map中的地理参考(R)；
4.matlab中NAN（空值）的运算规则：NAN值除赋值外的任何数学运算返回值为NAN，如NAN+1=NAN，NAN*0=NAN；
5.在map工具箱中，注意应用部分函数习惯是先引用纬度值，再引用经度值，需要注意分别此类函数。

二、边界线的制作：
1.读取文件：map=shaperead(filename)：读取shape文件，这里建议读取多边形文件，通用些。shape文件中国部分由于受到边境线问题制约，需要用我们国家自己的数据，国外部分的获取可以由清风提供的网址~另外要再做气候区或者若干地区的底图，用ArcGis分析工具等都可以做，可以见ArcGis相关教程，这里顺便贡献一个ArcGis的中文帮助网址。
2.绘制你所要绘制的图形：这个不多说，要画什么画什么，一般是绘制等值线图，用contour或contourf画，有人也许会问需不需要用contourm画，其实随便，两者的区别是前者用meshgrid生成的坐标系来绘图，后者通过地理参考(R)或者经纬度数据绘图，对matlab来说成图都一样。
3.加入边界线：
hold on
plot(map.X,map.Y,option)
hold off
就是普通的加线了，没什么特别的，option是可选的参数，美化一下绘图而已。

三、掩膜的制作：
1.读取文件，同上1；
2.生成掩膜文件：R=makerefmat('Rastersize',size(Z),'Latlim',[lat_first lat_last],'Lonlim',[lon_first lon_last])：其中Z是你要绘图的矩阵，Latlim和lonlim两个参数确定你要画全图的经纬度范围。这个函数是在后面转换数据时候用的，它的函数形式不止这一个（见help），但是我选这个用，具体原因大家可以尝试对比一下，因为直接使用定义经纬度起始点和经纬度点密度的方法生成的矩阵会出现少边界点的情况，而如上方式定义地理参考会产生一个缩减量使size(R)=size(Z)；
注意：鉴于有朋友提醒，2008版前的matlab中makerefmat可能没有'RasterSize'参数的设置，我试了一下，以下提供上面所用方式与R=makerefmat(x0,y0,dx,dy)通用方式的转换：
转换要点：
1.步长dx、dy由绘图范围和绘图格点总数的比值确定；
2.起始经纬度为绘图范围的左下角经纬度值加上半步长值。
例如：
R=makerefmat('RasterSize',size(Z),'Lonlim',[97 107],'Latlim',[21 30])
等价于
R=makerefmat(97+(107-97)/size(Z,1)/2,21+(30-21)/size(Z,2)/2,(107-97)/size(Z,1),(30-21)/size(Z,2))

1. >> R=makerefmat('RasterSize',size(Z),'Lonlim',[97 107],'Latlim',[21 30])
   
2. 
   
3. R =
   
4. 
   
5.          0    0.0989
   
6.     0.0990         0
   
7.    96.9505   20.9505
   
8. 
   
9. >> R=makerefmat(97,21,0.1,0.1)
   
10. 
    
11. R =
    
12. 
    
13.          0    0.1000
    
14.     0.1000         0
    
15.    96.9000   20.9000

复制代码

有关投影矩阵的内容在《经纬度和矩阵行列标的转换（基于MAP工具箱）及一个副高西伸脊点程序》里还有介绍，不过关注度有点低{:soso_e154:}，在这里打个广告吧~

MASK=vec2mtx(map.Y,map.X,Z,R,'filled')：产生了一个与绘图矩阵Z同阶的，map多边形区域外赋值为2，多边形边界赋值为1，多边形内赋值为0的矩阵，同样，vec2mtx函数也有其他输入形式，不过会有与上面地理参考生成时同样的问题，具体见help；
最后，将MASK中值>1的区域设置为nan，边界设置为0，掩膜制作完成；
3.绘图方法：绘图时绘制矩阵为(Z+MASK)，我这里使用的是加法方法，在2中掩膜设置方法改变一下就可以设置乘法方法的掩膜，原理同样是基于NAN与任何数的数学运算结果为NAN。

附上绘制上面那幅图的源程序，其中EOF_used是降水数据EOF的第一向量，gy_locat是站点的站点号和经纬度信息（测试数据见下，数据已直接提取了插值后的Z值了，跳过了代码的第二段）：


[lon lat]=meshgrid([97:0.1:107],[21:0.1:30]);
% Z=griddata(gy_locat(:,2),gy_locat(:,3),EOF_used(:,1),lon,lat,'v4');
yunnan=shaperead('yunnan.shp');

R=makerefmat('RasterSize',size(Z),'Lonlim',[97 107],'Latlim',[21 30]);
MASK=vec2mtx(yunnan.Y,yunnan.X,Z,R,'filled');
MASK(find(MASK>1))=nan;
MASK(find(MASK==1))=0;

contourf(lon,lat,Z+MASK,30);
shading flat
colorbar

hold on
plot(yunnan.X,yunnan.Y,'-k','linewidth',3)
hold off

不过在绘图过程中shading flat有时候会出现报错信息，这个就不是很清楚原因了，请各路大神指教~不过转换方法本身应该没什么问题。

写在最后：大家要是有什么好的matlab掩膜或者边界线制作的方法都跟上一贴吧，我也学习学习，交流一下。这个是折腾好久才弄出来的，只是map工具箱的冰山一角，希望大家也多有探索精神，利用身边的资源哈~相信matlab在地学绘图中的应用会走得更远。


应大家的要求，上传了一份处理过的测试数据：

我和小云

首先感谢楼主的分享，两年前的帖子帮了我大忙。其次关于MASK的方法，我在搜索国外的相关网站时貌似发现了更简单的做法，贴出来供大家参考下。用楼主之前的程序加以修改，黑体为关键部分。至于inpolygon函数大家自行help咯。
[lon lat]=meshgrid([97:0.1:107],[21:0.1:30]);
% Z=griddata(gy_locat(:,2),gy_locat(:,3),EOF_used(:,1),lon,lat,'v4');
yunnan=shaperead('yunnan.shp');

isin=inpolygon(lon,lat,yunnan.Lon,yunnan.Lat);
Z(~isin)=NaN;
contourf(lon,lat,Z,30);
shading flat
colorbar

hold on
plot(yunnan.X,yunnan.Y,'-k','linewidth',3)
hold off

zeroes

楼主有无试试regionpatch
代码如下（x和y为边界数据）：

function hc=regionpatch(varargin)
% patch the region outside the define region
% Usage: hc=regionpatch(|h,x,y,|color)
% h is axes handle. gca by default.
% color is the facecolor of out of region, white by default.

mv=0;temp=varargin{1};
if length(temp)==1&&ishandle(temp),  h=temp;mv=mv+1;end
x=varargin{mv+1};y=varargin{mv+2};
if nargin>mv+2, c=varargin{mv+3};end
if ~exist('h','var'), h=gca;end
if ~exist('c','var'), c='w';end

axes(h);hold on;
range=axis;

x=x(:);y=y(:);
[xi,i]=max(x);
x=[x(1:i);range([2 2 1 1 2 2])';x(i:end)];
y=[y(1:i);y(i);range([4 4 3 3])';y(i);y(i:end)];
hc=patch(x,y,c);
set(hc,'linestyle','none','edgecolor','none');

wlzhongouc

很棒啊 学习了！感谢楼主的分享！

斥鷃

wlzhongouc 发表于 2013-1-30 18:01

                               
登录/注册后可看大图

很棒啊 学习了！感谢楼主的分享！

很高兴得到大神夸奖，哈哈~

斥鷃

最得意一篇帖都没人回~求讨论~~~

shadowcathy

这个真的太强大了，一定要顶起，感谢楼主无私分享。

njzqxt

图真心漂亮。

kongfeng0824

matlab读取shpfile就可以使得底图数据瞬间丰富起来呀。

斥鷃

shadowcathy 发表于 2013-1-30 23:52

                               
登录/注册后可看大图

这个真的太强大了，一定要顶起，感谢楼主无私分享。

matlab软件很强大啦~那个下午就help map然后花几个小时一点点看找map工具箱里的函数，因为有ArcGis的基础，认为matlab肯定会有矢量和栅格数据互转功能，所以抱着一定能找到的信念去找，总算让我发现了，哈哈。不过英语真心要学好，太重要了~~~

斥鷃

njzqxt 发表于 2013-1-31 03:02

                               
登录/注册后可看大图

图真心漂亮。

O(∩_∩)O谢谢，matlab的colormap比较丰富，很适合绘图，所以特别不想浪费这个功能，以前都是在grads里画的，用着matlab不绘图心痛啊~

斥鷃

kongfeng0824 发表于 2013-1-31 07:00

                               
登录/注册后可看大图

matlab读取shpfile就可以使得底图数据瞬间丰富起来呀。

不过像你提的很多问题一样，数据读取方面还是一个大问题的，要弄通了就可以解决大部分问题了。我看map工具箱里面也有很多数据读取的函数，不过没有做深入研究，什么时候有空研究下读取吧~现在发现help比很多网上的资料都给力~~~