请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 238453|回复: 221

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

  [复制链接]

新浪微博达人勋

发表于 2016-9-7 18:28:48 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 Lighting 于 2020-5-17 16:40 编辑



写在最前

最近重新整理了一下白化的程序,把支持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.     addRequired(p, 'shapefile', validShape)
  38.     addRequired(p, 'masktype', @islogical);
  39.     addParameter(p, 'longitudes', defaultLon, validLLFcn);
  40.     addParameter(p, 'latitudes', defaultLat, validLLFcn);
  41.     addParameter(p, 'multiprovinces', defaultProv, @iscell);
  42.     addParameter(p, 'reversemask', defaultReve, @islogical);
  43.     addParameter(p, 'facecolor', defaultFace, validFaceFcn)
  44.     addParameter(p, 'edgecolor', defaultLine, validEdgeFcn)
  45.     addParameter(p, 'linewidth', defaultWidth, @isnumeric)

  46.     parse(p, shapefile, masktype, varargin{:});

  47.     if p.Results.masktype
  48.         longitudes = p.Results.longitudes;
  49.         latitudes = p.Results.latitudes;
  50.         if isempty(longitudes) || isempty(latitudes)
  51.             error('需指定经纬度范围!')
  52.         end
  53.         loninc = longitudes(1):longitudes(end);
  54.         latinc = latitudes(1):latitudes(end);
  55.         lone = ones(1,length(latinc)-2)*longitudes(end);
  56.         lons = ones(1,length(latinc)-1)*longitudes(1);
  57.         lats = ones(1,length(loninc)-1)*latitudes(1);
  58.         late = ones(1,length(loninc)-2)*latitudes(end);
  59.         LON = [loninc lone loninc(end:-1:1) lons NaN];
  60.         LAT = [lats latinc late latinc(end:-1:1) NaN];
  61.         [LON,LAT] = poly2cw(LON,LAT);
  62.       
  63.         if ischar(p.Results.shapefile)
  64.             s = shaperead(p.Results.shapefile);
  65.             if isfield(s,'NAME')
  66.                 allprovs = {s.NAME}';
  67.                 allprovnum = length(allprovs);
  68.                 maskindex = zeros(allprovnum,1);
  69.             elseif isfield(s,'NAME99')
  70.                 allprovs = {s.NAME99}';
  71.                 allprovnum = length(allprovs);
  72.                 maskindex = zeros(allprovnum,1);
  73.             end
  74.             long = [s.X];  lati = [s.Y];
  75.             % 获取指定的多个省份的边界点
  76.             if ~isempty(p.Results.multiprovinces)
  77.                 provinces = p.Results.multiprovinces;
  78.                 pronum = length(p.Results.multiprovinces);
  79.                 for ii = 1:pronum
  80.                     sindex = strfind(allprovs, provinces{ii});
  81.                     if any(cell2mat(sindex))
  82.                         for jj = 1: allprovnum
  83.                             if isempty(sindex{jj}) ; sindex{jj} = 0; end;
  84.                         end
  85.                         maskindex = maskindex | cell2mat(sindex);
  86.                     else
  87.                         error('所选省份 %s 不在指定的shp文件中或省份名错误!',provinces{ii})
  88.                     end
  89.                 end
  90.                 if p.Results.reversemask
  91.                     maskindex = ~maskindex;  % 指定省份以外的省份在shapefile中的索引
  92.                 end
  93.                 LON = [LON s(maskindex).X];
  94.                 LAT = [LAT s(maskindex).Y];
  95.             end
  96.         elseif iscell(p.Results.shapefile)
  97.             s = shaperead(p.Results.shapefile{1});
  98.             long = s.X;  lati = s.Y;
  99.             for jj = 2:length(p.Results.shapefile)
  100.                 s = shaperead(p.Results.shapefile{jj});
  101.                 [long, lati] = polybool('union',long,lati, s.X, s.Y);
  102.             end
  103.         end
  104.       
  105.         [xc, yc] = polybool('xor', LON,LAT , long, lati);
  106.         [f, v] = poly2fv(xc, yc);
  107.         %此处edgecolor必须为none,否则facecolor属性值将失效
  108.         patch('Faces', f, 'Vertices', v, 'FaceColor', p.Results.facecolor, 'edgecolor', 'none');
  109.         mapshow(long ,lati, 'color', p.Results.edgecolor, 'linewidth', p.Results.linewidth);
  110.         set(gca, 'xlim', p.Results.longitudes,'ylim', p.Results.latitudes)
  111.     else
  112.         if ischar(p.Results.shapefile)
  113.             mapshow(p.Results.shapefile,'displaytype','polygon','facecolor',p.Results.facecolor,'edgecolor',p.Results.edgecolor,'linewidth',p.Results.linewidth)
  114.         else
  115.             error('仅支持单文件白化!')
  116.         end
  117.     end
  118.     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. loninc = longitudes(1):longitudes(end);
  18. latinc = latitudes(1):latitudes(end);
  19. lone = ones(1,length(latinc)-2)*longitudes(end);
  20. lons = ones(1,length(latinc)-1)*longitudes(1);
  21. lats = ones(1,length(loninc)-1)*latitudes(1);
  22. late = ones(1,length(loninc)-2)*latitudes(end);
  23. LON.index = [loninc lone loninc(end:-1:1) lons NaN];
  24. LAT.index = [lats latinc late latinc(end:-1:1) NaN];
  25. LON.index = [LON.index ncst(:,1)'];
  26. LAT.index = [LAT.index ncst(:,2)'];

  27. coast.long = LON.index;
  28. coast.lat = LAT.index;

  29. 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. R=makerefmat('RasterSize',size(z),'Lonlim', lon,'Latlim',lat);
  10. contourfm(z,R,'linestyle','none')
  11. geoshow(flipud(coast.lat),flipud(coast.long),'DisplayType','polygon','facecolor','w')
复制代码


                               
登录/注册后可看大图

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

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

评分

参与人数 12金钱 +172 贡献 +29 体力 +50 收起 理由
li201657 + 2 很给力!
卷毛蓝鼻子 + 6 666
bqrrdzj + 5 很给力!
杜超杰 + 20 + 1 赞一个!
天坑一地坑 + 20 + 2 花点时间去泡妞
斥鷃 + 20 + 2 居然才看到~罪过罪过!
nuist2015 + 20 + 2
chenweibecome + 20 + 2 非常实用,点个赞!
mofangbao + 15 + 5
wlzhongouc + 20 + 2 很好,学习了
二爷名声在外 + 10 + 5 很给力!
kongfeng0824 + 14 + 8 + 50 很给力!

查看全部评分

本帖被以下淘专辑推荐:

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-12-6 10:17:29 | 显示全部楼层
本帖最后由 fengzhimu 于 2016-12-8 16:49 编辑
Lighting 发表于 2016-11-9 11:43
需要使用正文提供的maskMap函数进行白化。2012版本较低,请升级到2014以上,最好升级到2016,因为2014版 ...


楼主棒棒哒!我用2014a的白化外边的等值线还在,2016a的没法安装,因为电脑是32位的。。。改成32位的2015a的,试了下,也是可以完美白化哒!哎。。。如果用2014a的不画等值线,看上去没什么问题,但如果要画等值线,还是不行
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2016-9-7 18:50:31 | 显示全部楼层
学习!               
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 32430
发表于 2016-9-7 18:56:42 | 显示全部楼层
棒棒哒!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2016-9-8 13:29:44 | 显示全部楼层
好贴 学习了 好东西
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 0
发表于 2016-9-8 15:14:00 | 显示全部楼层
高手啊
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2016-9-8 15:38:18 | 显示全部楼层
学习一下,楼主好厉害
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-9-8 20:41:03 | 显示全部楼层
还收徒弟?
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2016-9-8 21:11:54 | 显示全部楼层

大神不要闹
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-9-8 22:42:02 | 显示全部楼层
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2016-9-9 10:43:02 | 显示全部楼层

我MATLAB 还一直没加载中国地图呢,改天向大神请教下
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表