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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 260|回复: 3

[程序设计] MATLAB利用inpolygon筛选特定省市数据

[复制链接]

新浪微博达人勋

发表于 6 天前 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 毕业两年的小白 于 2021-2-19 16:02 编辑

下载到的.nc数据,有时候需要筛选出指定省市的数据进行作图。网上搜了一圈,MATLAB主要用inpolygon这个函数,但是没有详细教程。
这里分享一下自己的学习成果。


                               
登录/注册后可看大图

目标:绘制黑吉辽三省9月的降雨总量
  1. clear  % 清除变量
  2. clc    % 清屏

  3. %% 1、读取数据
  4. info = ncinfo('precip_hourly_Sep.nc'); % 查看变量
  5. lon = ncread('precip_hourly_Sep.nc','longitude'); % 读取经度
  6. lat = ncread('precip_hourly_Sep.nc','latitude');  % 读取纬度
  7. time = ncread('precip_hourly_Sep.nc','time');  % 读取时间
  8. pr = ncread('precip_hourly_Sep.nc','tp');  % 读取9月的数据
  9. % shp = shaperead('D:\Work\My_MATLABTools\basic_map\bou2_4p.shp');
  10. % 注意这里要对原始shp手动处理,删掉不想要的省市。我这里只保留了黑吉辽三省,并另存成了shpPart。
  11. load('shpPart.mat')

  12. %% 2、数据筛选
  13. % 创建cLon cLat两个二维数组,用于存储经纬度,注意,都是二维的。
  14. [cLat,cLon]=meshgrid(lat,lon);

  15. % 获取shpPart内的点
  16. [hide]% in = inpolygon(xq,yq,xv,yv),xq yq是坐标点矩阵,二者大小一致;xv yv是边界矩阵,也要大小一致
  17. % 返回的in是一个布尔值矩阵,与xq yq大小一致,即每个点是否在xv yv组成的图形内
  18. % 这只是最基础的用法。更多奇技淫巧参考 https://ww2.mathworks.cn/help/matlab/ref/inpolygon.html
  19. boolIn=false(length(lon),length(lat));
  20. for i=1:length(shpPart)
  21.      in=inpolygon(cLon,cLat,shpPart(i).X, shpPart(i).Y);
  22.      % shpPart里其实有多组图形范围文件,用循环进行布尔值累加,获得所有东三省范围内的点
  23.      % 注意这里的boolIn经过加减运算自动转换成了double型,但其值仍然只有0/1
  24.      boolIn=boolIn+in;
  25. end[/hide]

  26. % 求降水的话,用下面办法很简单,直接利用布尔值false即为0的特点,将区外数值变成0即可
  27. % 但最好的办法还是变成nan,否则0值也会带来麻烦,尤其是温度等跨零数据
  28. sumPr=zeros(length(lat),length(lon));
  29. for i=1:length(time)
  30.      pr1(:,:,i)=pr(:,:,i).*boolIn*1000; %放大1000倍,由“米”转换为“毫米”
  31.      sumPr=sumPr+(pr1(:,:,i))';
  32. end
  33.       
  34. %将shpPart外的值设为nan
  35. sumPr2=zeros(length(lat),length(lon));
  36. for i=1:length(time)
  37.      for j=1:size(pr,1)
  38.         for k=1:size(pr,2)
  39.            if boolIn(j,k)==0
  40.                pr2(j,k,i)=nan;
  41.            else
  42.                pr2(j,k,i)=pr(j,k,i)*1000;%放大1000倍,由“米”转换为“毫米”
  43.            end
  44.         end
  45.      end
  46.      sumPr2=sumPr2+(pr2(:,:,i))';
  47. end


复制代码

效果图,没有美化,只是验证一下数据

效果图,没有美化,只是验证一下数据

效果图,没有细调,只是看看数据区域对不对

                               
登录/注册后可看大图


最后说说官方文档中的xq(in)算法
xq(in) yq(in) 直接获取了xq yq中的区内点坐标,但是,其弊端在于:
如果xq yq是二维数组,xq(in) yq(in)返回的不是二维数组,而是将xq yq 按 列 重 组 成的一维数组
对于本例,pr数据和lat lon数据是割裂的三个变量,用in方法只能判断lat lon是否在区内,无法判断pr
所以,为了理解简单,上面采用了逐点判断的方法
当然,用in方法把pr lat lon 都转换为一维数组的方法肯定是可行的,只要保证一一对应即可


                               
登录/注册后可看大图



原创不易,点个赞再走呗
转载请注明出处。



评分

参与人数 1金钱 +10 收起 理由
kongfeng0824 + 10 很给力!

查看全部评分

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

新浪微博达人勋

发表于 6 天前 | 显示全部楼层

回帖奖励 +1 金钱

很棒,码住,学习学习。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 6 天前 | 显示全部楼层
angetina 发表于 2021-2-19 12:06
很棒,码住,学习学习。

谢谢,共同学习
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 昨天 20:47 | 显示全部楼层
边界外的是不是掩抹处理掉会更好看点
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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