MATLAB利用inpolygon筛选特定省市数据

毕业两年的小白

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

                               
登录/注册后可看大图

目标：绘制黑吉辽三省9月的降雨总量

1. clear  % 清除变量
   
2. clc    % 清屏
   
3. 
   
4. %% 1、读取数据
   
5. info = ncinfo('precip_hourly_Sep.nc'); % 查看变量
   
6. lon = ncread('precip_hourly_Sep.nc','longitude'); % 读取经度
   
7. lat = ncread('precip_hourly_Sep.nc','latitude');  % 读取纬度
   
8. time = ncread('precip_hourly_Sep.nc','time');  % 读取时间
   
9. pr = ncread('precip_hourly_Sep.nc','tp');  % 读取9月的数据
   
10. % shp = shaperead('D:\Work\My_MATLABTools\basic_map\bou2_4p.shp');
    
11. % 注意这里要对原始shp手动处理，删掉不想要的省市。我这里只保留了黑吉辽三省，并另存成了shpPart。
    
12. load('shpPart.mat')
    
13. 
    
14. %% 2、数据筛选
    
15. % 创建cLon cLat两个二维数组，用于存储经纬度，注意，都是二维的。
    
16. [cLat,cLon]=meshgrid(lat,lon);
    
17. 
    
18. % 获取shpPart内的点
    
19. [hide]% in = inpolygon(xq,yq,xv,yv)，xq yq是坐标点矩阵，二者大小一致；xv yv是边界矩阵，也要大小一致
    
20. % 返回的in是一个布尔值矩阵，与xq yq大小一致，即每个点是否在xv yv组成的图形内
    
21. % 这只是最基础的用法。更多奇技淫巧参考 https://ww2.mathworks.cn/help/matlab/ref/inpolygon.html
    
22. boolIn=false(length(lon),length(lat));
    
23. for i=1:length(shpPart)
    
24.      in=inpolygon(cLon,cLat,shpPart(i).X, shpPart(i).Y);
    
25.      % shpPart里其实有多组图形范围文件，用循环进行布尔值累加，获得所有东三省范围内的点
    
26.      % 注意这里的boolIn经过加减运算自动转换成了double型，但其值仍然只有0/1
    
27.      boolIn=boolIn+in;
    
28. end[/hide]
    
29. 
    
30. % 求降水的话，用下面办法很简单，直接利用布尔值false即为0的特点，将区外数值变成0即可
    
31. % 但最好的办法还是变成nan，否则0值也会带来麻烦，尤其是温度等跨零数据
    
32. sumPr=zeros(length(lat),length(lon));
    
33. for i=1:length(time)
    
34.      pr1(:,:,i)=pr(:,:,i).*boolIn*1000; %放大1000倍，由“米”转换为“毫米”
    
35.      sumPr=sumPr+(pr1(:,:,i))';
    
36. end
    
37.       
    
38. %将shpPart外的值设为nan
    
39. sumPr2=zeros(length(lat),length(lon));
    
40. for i=1:length(time)
    
41.      for j=1:size(pr,1)
    
42.         for k=1:size(pr,2)
    
43.            if boolIn(j,k)==0
    
44.                pr2(j,k,i)=nan;
    
45.            else
    
46.                pr2(j,k,i)=pr(j,k,i)*1000;%放大1000倍，由“米”转换为“毫米”
    
47.            end
    
48.         end
    
49.      end
    
50.      sumPr2=sumPr2+(pr2(:,:,i))';
    
51. end
    
52. 
    
53. 
    

复制代码

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

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

                               
登录/注册后可看大图

最后说说官方文档中的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 都转换为一维数组的方法肯定是可行的，只要保证一一对应即可

                               
登录/注册后可看大图

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

天坑一地坑

S=shaperead(shp);
bou2_4lx=[S(:).X];
bou2_4ly= [S(:).Y];

in = inpolygon(lon,lat,bou2_4lx,bou2_4ly);
ini = double(in);
ini(ini==0) = nan;

lon1 = lon.*ini;
lat1 = lat.*ini;
rain1 = rain.*ini;
楼主参考下，少用循环

angetina

很棒，码住，学习学习。

毕业两年的小白

angetina 发表于 2021-2-19 12:06
很棒，码住，学习学习。

谢谢，共同学习

新八不是废材男

边界外的是不是掩抹处理掉会更好看点

毕业两年的小白

新八不是废材男 发表于 2021-2-24 20:47
边界外的是不是掩抹处理掉会更好看点

想要出图效果好的话确实需要掩膜处理。不过我这边只需要数据，后续还要求和求平均等运算。所以凑合看看吧

baishizhou

请问降水数据是国家气象中心注册后下载的数据还是类似欧洲ERA5之类在网上下载的数据？

772615260

只保留边界内数据的最好函数我还是觉得是regionpatch，你提到的这个还有其他的函数我用过，主要是存在边界线上受数据分辨率不能完美裁切的问题。我说的那个函数不会存在这个问题

淘气包马小跳

请问如何得到特定省市的shp文件？

毕业两年的小白

baishizhou 发表于 2021-2-26 12:50
请问降水数据是国家气象中心注册后下载的数据还是类似欧洲ERA5之类在网上下载的数据？

era5上的，国家局的数据好像是站点数据吧？记不太清了

毕业两年的小白

淘气包马小跳 发表于 2021-3-18 09:09
请问如何得到特定省市的shp文件？

我的是笨办法：找一个全国的，手动筛选一下，把不要的删掉就好了。