[支持m_map白化] MATLAB对地图进行白化

Lighting

写在最前：

最近重新整理了一下白化的程序，把支持m_map工具箱的功能完善了一下。内容都发在 https://i-lightning.cn/2020/05/matlab_mask/ 这里了，代码也都开源了。此外还加入了一个简单的函数m_mapshow，支持添加任意shapefile文件到地图中。


==================================================
       感谢那些对程序的完善提出宝贵意见的朋友。
        @kc121   对程序的进一步修改提出了宝贵的意见。在此表示由衷的感谢！
==================================================

       最近看到python板块提供了不少对地图进行完美白化的方法。
       比如@平流层的萝卜     python完美白化
              @非对称   利用scipy和sympy扩大三角网格形成的轮廓
       也看到了MATLAB版对白化功能的讨论，也给出了一些方法。
       比如 @斥鷃     提供的        matlab中地图边界与掩膜的实现
               @raul928    提供的        也提供一个去地图外边界的方法
       上述两种方法都能很好的对地图进行白化。
       本来想将@非对称 的方法用matlab实现的，大部分代码已经写好，但有一点不是太明白就暂时放置在一旁了。然后在无意之间发现了mapshow函数利用 facecolor属性可以设置区域内白化，经过一番探索，便有了下面的白化方法。

       1）  如果需要白化某一区域内的话，只需要使用mapshow函数即可，以shapefile文件为例                     

1. z = peaks(1000);
   
2. lon = [60 150];
   
3. lat = [10 60];
   
4. [LON,LAT] = meshgrid(linspace(lon(1),lon(2),1000), linspace(lat(1), lat(2),1000));
   
5. c = contourf(LON,LAT,z,'linestyle', 'none');
   
6. mapshow('E:\MATLAB\shp\bou2_4p.shp','displaytype','polygon', 'facecolor','w')

复制代码

               
         mapshow 的 displaytype 属性默认为 'line'，这里要设置为'polygon'。

                               
登录/注册后可看大图

      2） 白化区域外部分

   

1. function maskMap(shapefile, masktype, varargin)
   
2.     %  对所绘制图形进行白化
   
3.     %  输入参数：
   
4.     %      shapefile  :  shapefile文件。 字符串型或元胞型
   
5.     %              为元胞数组时可通过指定多个省份的shp文件进行白化。
   
6.     %              比如要白化江苏省，江西省，黑龙江省。使用元胞数组白化多个省份，需要手动添加中国边界地图
   
7.     %              shapefile = {'path\江苏省.shp', 'path\江西省.shp', 'path\黑龙江省.shp'}
   
8.     %                  其中path为shp文件所在路径
   
9.     %       masktype  : 逻辑值。
   
10.     %           true  :  对区域外进行白化
    
11.     %                可选参数：
    
12.     %                    longitudes ： 图形的经度范围，二元素向量。 数值型。
    
13.     %                    latitudes  :  图形的纬度范围，二元素向量。 数值型。
    
14.     %                    multiprovinces : 元胞数组。
    
15.     %                    用于指定要白化的多个区域。使用此参数需要shapefile为字符型且要包含省份名称
    
16.     %                    此处建议大家使用 bou2_4p.shp 进行白化
    
17.     %                              比如要白化江苏省，江西省，黑龙江省，则
    
18.     %                         multipro = {'江苏省','江西省','黑龙江省'}
    
19.     %                    reversemask  :  指定要白化的是 multiprovinces 元胞中指定的区域还是
    
20.     %                                  元胞中指定的区域以外的部分。 逻辑值。
    
21.     %                                true  :  表示白化元胞中指定的区域。 默认值。
    
22.     %                                false :  表示白化元胞中指定区域以外的部分。
    
23.     %           false :  对区域内进行白化
    
24.     %      共享可选参数值对：
    
25.     %           facecolor :  用于指定白化区域的颜色。 默认为白色。 字符型或数值型。
    
26.     %                        取值为 'none'时，不起作用。数值型时为RGB三元组。
    
27.     %           edgecolor :  指定的边界线的颜色。 默认为黑色。 字符型或数值型。同上。
    
28.     %           linewidth :  边界线宽度。 默认为1. 数值型。
    
29.     %%
    
30.     p = inputParser;
    
31.     validShape = @(x) ischar(x) || iscell(x);
    
32.     validLLFcn = @(x) isvector(x) && length(x) == 2;
    
33.     validFaceFcn = @(x) (ischar(x) && ~strcmp(x, 'none')) || (isvector(x) && length(x) == 3);
    
34.     validEdgeFcn = @(x) ischar(x) || (isvector(x) && length(x) == 3);
    
35.     defaultLon = [];  defaultLat = []; defaultProv = {}; defaultReve = true;
    
36.     defaultFace = 'w';  defaultLine = 'k';  defaultWidth = 1;
    
37. 
    
38.     addRequired(p, 'shapefile', validShape)
    
39.     addRequired(p, 'masktype', @islogical);
    
40.     addParameter(p, 'longitudes', defaultLon, validLLFcn);
    
41.     addParameter(p, 'latitudes', defaultLat, validLLFcn);
    
42.     addParameter(p, 'multiprovinces', defaultProv, @iscell);
    
43.     addParameter(p, 'reversemask', defaultReve, @islogical);
    
44.     addParameter(p, 'facecolor', defaultFace, validFaceFcn)
    
45.     addParameter(p, 'edgecolor', defaultLine, validEdgeFcn)
    
46.     addParameter(p, 'linewidth', defaultWidth, @isnumeric)
    
47. 
    
48.     parse(p, shapefile, masktype, varargin{:});
    
49. 
    
50.     if p.Results.masktype
    
51.         longitudes = p.Results.longitudes;
    
52.         latitudes = p.Results.latitudes;
    
53.         if isempty(longitudes) || isempty(latitudes)
    
54.             error('需指定经纬度范围！')
    
55.         end
    
56.         loninc = longitudes(1):longitudes(end);
    
57.         latinc = latitudes(1):latitudes(end);
    
58.         lone = ones(1,length(latinc)-2)*longitudes(end);
    
59.         lons = ones(1,length(latinc)-1)*longitudes(1);
    
60.         lats = ones(1,length(loninc)-1)*latitudes(1);
    
61.         late = ones(1,length(loninc)-2)*latitudes(end);
    
62.         LON = [loninc lone loninc(end:-1:1) lons NaN];
    
63.         LAT = [lats latinc late latinc(end:-1:1) NaN];
    
64.         [LON,LAT] = poly2cw(LON,LAT);
    
65.       
    
66.         if ischar(p.Results.shapefile)
    
67.             s = shaperead(p.Results.shapefile);
    
68.             if isfield(s,'NAME')
    
69.                 allprovs = {s.NAME}';
    
70.                 allprovnum = length(allprovs);
    
71.                 maskindex = zeros(allprovnum,1);
    
72.             elseif isfield(s,'NAME99')
    
73.                 allprovs = {s.NAME99}';
    
74.                 allprovnum = length(allprovs);
    
75.                 maskindex = zeros(allprovnum,1);
    
76.             end
    
77.             long = [s.X];  lati = [s.Y];
    
78.             % 获取指定的多个省份的边界点
    
79.             if ~isempty(p.Results.multiprovinces)
    
80.                 provinces = p.Results.multiprovinces;
    
81.                 pronum = length(p.Results.multiprovinces);
    
82.                 for ii = 1:pronum
    
83.                     sindex = strfind(allprovs, provinces{ii});
    
84.                     if any(cell2mat(sindex))
    
85.                         for jj = 1: allprovnum
    
86.                             if isempty(sindex{jj}) ; sindex{jj} = 0; end;
    
87.                         end
    
88.                         maskindex = maskindex | cell2mat(sindex);
    
89.                     else
    
90.                         error('所选省份 %s 不在指定的shp文件中或省份名错误！',provinces{ii})
    
91.                     end
    
92.                 end
    
93.                 if p.Results.reversemask
    
94.                     maskindex = ~maskindex;  % 指定省份以外的省份在shapefile中的索引
    
95.                 end
    
96.                 LON = [LON s(maskindex).X];
    
97.                 LAT = [LAT s(maskindex).Y];
    
98.             end
    
99.         elseif iscell(p.Results.shapefile)
    
100.             s = shaperead(p.Results.shapefile{1});
     
101.             long = s.X;  lati = s.Y;
     
102.             for jj = 2:length(p.Results.shapefile)
     
103.                 s = shaperead(p.Results.shapefile{jj});
     
104.                 [long, lati] = polybool('union',long,lati, s.X, s.Y);
     
105.             end
     
106.         end
     
107.       
     
108.         [xc, yc] = polybool('xor', LON,LAT , long, lati);
     
109.         [f, v] = poly2fv(xc, yc);
     
110.         %此处edgecolor必须为none，否则facecolor属性值将失效
     
111.         patch('Faces', f, 'Vertices', v, 'FaceColor', p.Results.facecolor, 'edgecolor', 'none');
     
112.         mapshow(long ,lati, 'color', p.Results.edgecolor, 'linewidth', p.Results.linewidth);
     
113.         set(gca, 'xlim', p.Results.longitudes,'ylim', p.Results.latitudes)
     
114.     else
     
115.         if ischar(p.Results.shapefile)
     
116.             mapshow(p.Results.shapefile,'displaytype','polygon','facecolor',p.Results.facecolor,'edgecolor',p.Results.edgecolor,'linewidth',p.Results.linewidth)
     
117.         else
     
118.             error('仅支持单文件白化！')
     
119.         end
     
120.     end
     
121.     end

复制代码

        maskMap函数提供了一个逻辑值选项用于选择区域内白化还是区域外白化。如果是区域内白化，只需要提供shapefile和masktype两个参数即可，如果选择区域外白化，除上述两个参数之外，还需要提供‘longitudes’ 和 'latitudes' 两个参数，调用方式和matlab内置函数调用方式类似。   

        

                               
登录/注册后可看大图

        参考@平流层的萝卜 的python完美白化的帖子，也支持了白化几个省份的部分。

                               
登录/注册后可看大图

         无论是白化区域内还是白化区域外，均可以通过设置 facecolor， edgecolor, linewidth等属性值设置白化区域的颜色，边界线颜色，线宽等。鉴于白化时通常只需要改变这几个属性，目前仅支持这几个属性。
         下图是放大之后的，可以看出已经完美白化了，没有任何锯齿存在！

                               
登录/注册后可看大图

            直接调用上述函数即可。使用精确到县级的shp文件。示例代码如下：

1. z = peaks(1000);
   
2. lon = [60 135];
   
3. lat = [10 60];
   
4. [LON,LAT] = meshgrid(linspace(lon(1),lon(2),1000), linspace(lat(1), lat(2),1000));
   
5. c = contourf(LON,LAT,z,'linestyle', 'none');
   
6. maskMap('E:\MATLAB\shp\BOUNT_poly.shp',true,'lon',lon,'lat',lat,'linewidth',0.5, 'edgecolor','b')
   

复制代码

         

                               
登录/注册后可看大图

          可以看出，使用县级文件进行白化也毫无问题！         
=================================================
          2016.9.8   解决调用 maskMap 函数白化区域内出错的问题。程序已经更改，复制替换源程序即可。
     2016.9.9    解决白化时设置edgecolor 属性值为非'none'值时出现一条直线的问题。
     2016.9.14   解决设置等值线线型后无法白化区域外等值线的问题(2016a版本可行,2014a版本不可行。应该是语法有更新导致)；同时通过更改白化方法解决了设置线宽过大时或颜色不为黑色时，坐标轴会出现粗线或多颜色相间的问题。
     2016.9.23   释放使用RGB三元组设置facecolor 和 edgecolor 属性值。

      把重新编辑后的内容以PDF文档的形式上传，其中包括了对函数的详细解释，以及最新更新（当然最新的程序也已更新在 2）中）。白化算法有更新，但是文档内程序未更新，不过用法是相同的。

       注意：如果使用过程中出现以下错误提示，说明你使用的MATLAB版本较低，请更换到更高版本，建议更新到2015版之后。

1. Undefined function 'addParameter' for input arguments of type 'inputParser'.

复制代码

==== 插播更新利用m_map工具箱生成海岸线数据进行白化 ====
     最近看了一下m_map工具箱中关于生成海岸线的部分，发现可以自定义常用的绘图区，结合上面的白化方法，便想着利用m_map工具箱和 gshhs 高精度地形数据制作用于白化的数据。
      废话不多说，直接上程序。
     

1. lon = [60 150]; lat = [10 60];
   
2. FILNAME = 'private/gshhs_i.b';
   
3. h = figure;
   
4. m_proj('mercator', 'longitudes', lon, 'latitudes', lat)
   
5. m_grid('box','on')
   
6. % mu_coast一般情况下是不可调用的，可把private文件夹下的函数复制到另一个文件夹下
   
7. % 把该文件夹加入到搜索路径中
   
8. mu_coast('i',FILNAME,'tag','m_gshhs_i');  
   
9. m_gshhs_i('save','test')  % 保存 ncst , area ,k 数据到 test.mat 文件中
   
10. set(h , 'visible','off')   % 不显示图形。如果想显示图形到窗口可注释此语句
    

复制代码

     上面这部分用于生成绘图海岸线数据。
      下一步就是调用函数生成用于白化的海岸线数据。

1. function coast = creatCoast(longitudes,latitudes, ncst)
   
2. % 生成用于掩膜的经纬度数组
   
3. %   输入参数：
   
4. %       longitudes  : 经度起止坐标. 二元素向量.
   
5. %       latitudes   : 纬度起止坐标. 二元素向量.
   
6. %          ncst     : 由m_gshhs_* 函数生成的高精度经纬度坐标信息
   
7. %                     N X 2 数组。 其中第一列为经度信息
   
8. %                                      第二列为纬度信息
   
9. %   输出参数：
   
10. %        coast 为结构体变量. 包含了用于白化的经纬度信息
    
11. %            coast.long  存储的是经度信息
    
12. %            coast.lat   存储的是纬度信息
    
13. %%
    
14. validateattributes(longitudes,{'double'}, {'vector','numel',2})
    
15. validateattributes(latitudes,{'double'}, {'vector', 'numel',2})
    
16. validateattributes(ncst,{'double'}, {'ncols',2})
    
17. 
    
18. loninc = longitudes(1):longitudes(end);
    
19. latinc = latitudes(1):latitudes(end);
    
20. lone = ones(1,length(latinc)-2)*longitudes(end);
    
21. lons = ones(1,length(latinc)-1)*longitudes(1);
    
22. lats = ones(1,length(loninc)-1)*latitudes(1);
    
23. late = ones(1,length(loninc)-2)*latitudes(end);
    
24. LON.index = [loninc lone loninc(end:-1:1) lons NaN];
    
25. LAT.index = [lats latinc late latinc(end:-1:1) NaN];
    
26. LON.index = [LON.index ncst(:,1)'];
    
27. LAT.index = [LAT.index ncst(:,2)'];
    
28. 
    
29. coast.long = LON.index;
    
30. coast.lat = LAT.index;
    
31. 
    
32. end

复制代码

         注意： 从头到尾所使用的经纬度的起止范围要是相同的，不要频繁更改经纬度起止范围。

1. z = peaks(1000);
   
2. lon = [60 150];
   
3. lat = [10 60];
   
4. [LON,LAT] = meshgrid(linspace(lon(1),lon(2),1000), linspace(lat(1), lat(2),1000));
   
5. figure
   
6. contourf(LON,LAT,z,'linestyle' , 'none')
   
7. load test.mat
   
8. coast = creatCoast(lon,lat,ncst);
   
9. % displaytype 要设置为 'polygon'
   
10. mapshow(coast.long,coast.lat, 'Displaytype', 'polygon', 'facecolor', 'w', 'edgecolor', 'k', 'linewidth', 1);

复制代码

       执行后会发现有一点问题！竟然凭空的多出一块出来！应该是自定义海岸线时，多了一部分，但是究竟是怎么回事暂时不清楚。

                               
登录/注册后可看大图

不用担心，解决这个问题很简单！只要执行以下语句限制坐标轴范围就行了！

1. set(gca, 'xlim', lon, 'ylim', lat)

复制代码

                               
登录/注册后可看大图

          除直接限制坐标轴范围外，还可以进行数据筛选，把超出范围部分找出，然后去除即可。
        以上部分仅供参考。

=========  提供一种 MATLAB Central 中提到的一种白化方法 ==============
        首先附上链接： ocean mask
        
        文中提到的白化方法是利用MATLAB自带的coast.mat文件进行白化，主要是反转全球海岸线经纬度向量的顶点方向进行白化。尝试着反转中国的shp数据和gshhs数据进行白化，但是均以失败告终。如果大家有兴趣可以尝试一下。成功了的话欢迎科普

        测试代码：

1. coast = load('coast');
   
2. z = peaks(1000);
   
3. lon = [60 150];
   
4. lat = [10 60];
   
5. figure('Color','w')
   
6. axesm('mercator','MapLatLimit',lat,'MapLonLimit', lon)
   
7. axis off; framem on; gridm on; mlabel on; plabel on;
   
8. setm(gca,'MLabelLocation',60)
   
9. 
   
10. R=makerefmat('RasterSize',size(z),'Lonlim', lon,'Latlim',lat);
    
11. contourfm(z,R,'linestyle','none')
    
12. geoshow(flipud(coast.lat),flipud(coast.long),'DisplayType','polygon','facecolor','w')

复制代码

                               
登录/注册后可看大图

      使用内置的海岸线数据无法准确显示中国的边界线。

maskMap.pdf (780.95 KB, 下载次数: 477)

2016-9-11 16:19 上传

点击文件名下载附件 chinamap.rar (32.89 MB, 下载次数: 1391)

2016-9-18 13:19 上传

点击文件名下载附件
边界图

fengzhimu

Lighting 发表于 2016-11-9 11:43
需要使用正文提供的maskMap函数进行白化。2012版本较低，请升级到2014以上，最好升级到2016，因为2014版 ...

楼主棒棒哒！我用2014a的白化外边的等值线还在，2016a的没法安装，因为电脑是32位的。。。改成32位的2015a的，试了下，也是可以完美白化哒！哎。。。如果用2014a的不画等值线，看上去没什么问题，但如果要画等值线，还是不行

kongfeng0824

学习！               

二爷名声在外

棒棒哒！

aas314

好贴 学习了 好东西

jiajian0512

高手啊

雾失流年

学习一下，楼主好厉害

天坑一地坑

还收徒弟？

Lighting

天坑一地坑 发表于 2016-9-8 20:41
还收徒弟？

大神不要闹

tangjian

天坑一地坑

Lighting 发表于 2016-9-8 21:11
大神不要闹

我MATLAB 还一直没加载中国地图呢，改天向大神请教下