Python截取不规则区域内格点

Masterpiece

Python截取不规则区域内格点

做EOF的时候，一般都要区域内的格点进行分析，以避免区域外的信息分担EOF权重。
取区域内的格点的做法是我的舍友涌哥教我的方法，大家可以借鉴一下思路和做法吧！

首先要准备好区域范围的shp文件，文件提供了边界的经纬度，在程序中用shapefile库进行读取:

1. sf=shapefile.Reader("D:\\dt\\SC\\Sichuan_province")
   
2. Shapes = sf.shapes()
   
3. boundary = Shapes[0].points

复制代码

第二步就是关键了，这步是用matplotlib根据得到的经纬度点，建立一个path：

1. from matplotlib.path import Path
   
2. p= Path(boundary)

复制代码

这时候得到的p就是一个path边界的类变量，其中包含类函数contains_points来判断所给格点是否在区域内。例如:

1. x=100
   
2. y=30
   
3. p.contains_points([[x,y]])

复制代码

如果这个（100,30）的点在path区域内，within得到值True，反之是False。

并且contains_points可以将复数个位置点以(n x 2)的数组的形式输入，返回一个对应的布尔数组，又例如：

1. grid=np.array( [[[100,30],[100,50]],[[100,30],[100,30]]] )
   
2. within=p.contains_points(grid.reshape(4,2))

复制代码

这时类函数返回within得到一个布尔数组：

1. [ True False  True  True]

复制代码

这时候就可以利用这个布尔数组对numpy数组进行取数了。
比如我要画散点图，那必须知道散点对应的经纬度位置。先利用nc文件的lon、lat信息建立二维的经纬度格点：

1. lons, lats = np.meshgrid(lon, lat)

复制代码

这时候lons和lats就是一个(n,m)的数组，对应着每个格点的经度或纬度，以lons为例：

1. [[  0.    2.5   5.  ... 352.5 355.  357.5]
   
2. [  0.    2.5   5.  ... 352.5 355.  357.5]
   
3. [  0.    2.5   5.  ... 352.5 355.  357.5]
   
4. ...
   
5. [  0.    2.5   5.  ... 352.5 355.  357.5]
   
6. [  0.    2.5   5.  ... 352.5 355.  357.5]
   
7. [  0.    2.5   5.  ... 352.5 355.  357.5]]

复制代码

然后再利用布尔数组mask，将其中为True的对应格点的经度和纬度返回来，并且做一次坐标投影变换：

1. grib=np.stack((lons,lats),axis=-1)
   
2. mask=p.contains_points(grib.reshape(lat_num*lon_num,2))
   
3. mask=mask.reshape(lat_num,lon_num)
   
4. x,y=m(lons[mask],lats[mask])

复制代码

这时候在边界内的格点的图上投影坐标x,y就得到了，于是画散点图！

1. m.scatter(x,y,s=50, marker='o', color='g', zorder=6)

复制代码

这样就可以在图上标出在区域内的格点位置了。
所以可以使用类似思路将数据提取出来。

开篇图片所使用的代码，以四川省为例。
所用到的文件我放在附件里边。

1. import matplotlib.pyplot as plt
   
2. from mpl_toolkits.basemap import Basemap
   
3. import numpy as np
   
4. import netCDF4 as nc
   
5. from matplotlib.path import Path
   
6. import shapefile
   
7. 
   
8. af=nc.Dataset('test.nc')
   
9. lat=af.variables['latitude'][:]
   
10. lon=af.variables['longitude'][:]
    
11. p72=af.variables['p72.162'][0,:,:]/100
    
12. lat_num=len(lat)
    
13. lon_num=len(lon)
    
14. 
    
15. sf=shapefile.Reader("D:\\dt\\SC\\Sichuan_province")
    
16. Shapes = sf.shapes()
    
17. boundary = Shapes[0].points
    
18. p= Path(boundary)
    
19. 
    
20. lons, lats = np.meshgrid(lon, lat)
    
21. 
    
22. grib=np.stack((lons,lats),axis=-1) #变为(n,m,2)的经纬度的点
    
23. mask=p.contains_points(grib.reshape(lat_num*lon_num,2)) #变为(n*m,2)的经纬度的点,否则出错
    
24. # 如果点在路径之内，对应的mask赋值为False，否则为Ture
    
25. mask=mask.reshape(lat_num,lon_num)# 恢复到 (n,m)
    
26. 
    
27. print(lons[mask])
    
28. lon_leftup=90-2.5/2;lat_leftup=40+2.5/2
    
29. lon_rightdown=120+2.5/2;lat_rightdown=20-2.5/2
    
30. 
    
31. fig, ax = plt.subplots(figsize=(14,9))
    
32. 
    
33. m = Basemap(projection='cyl',
    
34.             llcrnrlat=lat_rightdown, urcrnrlat=lat_leftup,
    
35.             llcrnrlon=lon_leftup, urcrnrlon=lon_rightdown,
    
36.             resolution='i')
    
37. m.drawcoastlines(linewidth=0.3, color='black')
    
38. m.readshapefile('D:\\dt\\SC\\Sichuan_province', 'sc', color='k', linewidth=0.6)
    
39. 
    
40. x,y=m(lons[mask],lats[mask])
    
41. m.scatter(x,y,s=50, marker='o', color='g', zorder=6)
    
42. 
    
43. x,y=m(lons[~mask],lats[~mask])
    
44. m.scatter(x,y,s=50, marker='x', color='r', zorder=6)
    
45. 
    
46. parallels = np.arange(20,45,5) #lat
    
47. m.drawparallels(parallels,labels=[True,False,False,False],linewidth=0.6,dashes=[1,4])
    
48. meridians = np.arange(90,125,5) #lon
    
49. m.drawmeridians(meridians,labels=[False,False,False,True],linewidth=0.6,dashes=[1,4])
    
50. 
    
51. plt.yticks(parallels,len(parallels)*[''])
    
52. plt.xticks(meridians,len(meridians)*[''])
    
53. 
    
54. plt.savefig('grid_within.png',dpi=300,bbox_inches='tight')
    
55. 
    
56. plt.show()

复制代码

--所用文件--
grid_within.rar (136.94 KB, 下载次数: 41)

2020-4-27 19:31 上传

点击文件名下载附件

Masterpiece

言墨 发表于 2020-5-3 21:30
pre=af['pr'][0,mask]*24*3600  # pr(time,lat,lon)
cn=ax.contourf(lons[mask],lats[mask],pre,transform ...

有可能是lons[mask]这个的原因，因为mask取数之后，lons[mask]是一个一维的

lm8005507771

感谢! 收藏了，以后肯定有机会用上的

likailing2008

好文章，以后一定会用到，现在强调精细化，各种形状的判断和操作会越来越多

墨青衫

mark
一下，以后有用

失落的积雨云

楼主你好，想请问一下，按照你的这个方法，想把数据提取出来的话具体应该怎么操作呢？我原本想着用lons[mask]作为数据的经纬度索引，却发现经纬度的索引应该是从零开始的整数，这个思路行不通。

Masterpiece

失落的积雨云 发表于 2020-5-3 00:41
楼主你好，想请问一下，按照你的这个方法，想把数据提取出来的话具体应该怎么操作呢？我原本想着用lons[mas ...

numpy的数组不能用经纬度直接mask
如果要mask数据出来，直接：...['snow_depth'][ : , mask] 就好

失落的积雨云

Masterpiece 发表于 2020-5-3 07:58
numpy的数组不能用经纬度直接mask
如果要mask数据出来，直接：...['snow_depth'][ : , mask] 就好

问题解决了，非常感谢！

言墨

你好  请问22行的目的是什么呢

Masterpiece

言墨 发表于 2020-5-3 20:59
你好  请问22行的目的是什么呢

这个的目的是生成一个经纬度(n,m,2)的数组，前两维是经向和纬向格点数，第三维是经度和纬度
因为contains_points需要是经度纬度成对成对的输入：
[[点A经度,点A纬度]
[点B经度,点B纬度]
[点C经度,点C纬度]
[点D经度,点D纬度]]

言墨

pre=af['pr'][0,mask]*24*3600  # pr(time,lat,lon)
cn=ax.contourf(lons[mask],lats[mask],pre,transform=proj,levels=levs)

楼主打扰了，咨询一个问题，我想把一个省的数据mask出来然后画shade图，这么写一直报错，一直不知道怎么回事，请您指导下。
报错内容为：IndexError: Index cannot be multidimensional