Python解决任意区域的格点包含和面积平均问题

平流层的萝卜

气象中应该会经常遇到算一个区域内的平均值问题。有时候区域是矩形，则grads的aave可以直接解决；有时候区域是不规则的，比如说河南省地区。如果用格点资料算河南省的平均降水量，就涉及到不规则区域的面积平均问题。

本帖记录一下使用Python解决任意polygon（任意区域）的简单面积平均问题。

PS:所谓简单面积平均，即不涉及面积内的经纬度和权重等复杂运算，仅仅是面积内的格点的值的平均。

具体分为以下几步：
1）确定要算平均的区域，将该区域的point所围成的polygon多边形转换为matplotlib中的path。可以从shp文件中获取感兴趣的区域，比如西藏自治区（注意藏南和阿克赛钦部分），北京大兴县，江浙沪三地总辖区等。这里用shp文件提取感兴趣区时，使用的代码片段出自我以前的python maskout帖子：http://bbs.06climate.com/forum.php?mod=viewthread&tid=42437
2）判断有哪些点落在这个path里。此处利用Path里的contains_point方法，得到一个在或不在区域内的bool型array。
3）利用上一步的array，将区域外的降水量（或者其他变量）的值转化为np.nan。结合np.nanmean方法可以排除np.nan值，从而直接计算降水量（或者其他变量）的简单面积均值。

由于手头的降水数据较大，所以给一个计算海拔高度的面积平均的例子，代码和示例文件如下：

1. #coding=utf-8
   
2. import numpy as np
   
3. import shapefile
   
4. import netCDF4 as nc
   
5. from matplotlib.path import Path
   
6. from mpl_toolkits.basemap import Basemap
   
7. import matplotlib.pyplot as plt
   
8. 
   
9. ncfile=nc.Dataset(r'surface_z_0.25_46_36_115_128.nc')
   
10. lons=ncfile.variables['longitude'][:]
    
11. lats=ncfile.variables['latitude'][:]
    
12. z=ncfile.variables['z'][0,:,:]
    
13. xx,yy=np.meshgrid(lons,lats)
    
14. x,y=xx.flatten(),yy.flatten()
    
15. points=np.vstack((x,y)).T
    
16. 
    
17. # print yy
    
18. # exit()
    
19. 
    
20. sf=shapefile.Reader(r'bou2_4p',' ')
    
21. vertices = []
    
22. codes = []
    
23. for shape_rec in sf.shapeRecords():
    
24.     # print record
    
25.     if shape_rec.record[4] in [210000,110000,120000]:   ###辽宁、北京和天津的行政编码，百度可查
    
26.         pts = shape_rec.shape.points
    
27.         prt = list(shape_rec.shape.parts) + [len(pts)]
    
28.         for i in range(len(prt) - 1):
    
29.             for j in range(prt[i], prt[i+1]):
    
30.                 vertices.append((pts[j][0], pts[j][1]))
    
31.             codes += [Path.MOVETO]
    
32.             codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
    
33.             codes += [Path.CLOSEPOLY]
    
34.         clip = Path(vertices, codes)
    
35. 
    
36. polygon_of_Liaoning_beijing_tianjin= clip
    
37. # exit()
    
38. 
    
39. bool_grid_1d = polygon_of_Liaoning_beijing_tianjin.contains_points(points)
    
40. bool_grid_original_shape=bool_grid_1d.reshape(xx.shape)
    
41. 
    
42. 
    
43. ###########################################################################
    
44. ###### Highlighting grids within the area using scatterplot################
    
45. ###########################################################################
    
46. m = Basemap(llcrnrlat=lats[-1], llcrnrlon=lons[0], urcrnrlat=lats[0], urcrnrlon=lons[-1])
    
47. m.readshapefile(r'E:\shpmap\bou2_4p', ' ', color='grey', linewidth=2)
    
48. m.drawparallels(np.linspace(lats[0],lats[-1],10), labels=[1, 0, 0, 0], linewidth=0.5, fontsize=20)
    
49. m.drawmeridians(np.linspace(lons[0],lons[-1],10), labels=[0, 0, 0, 1], linewidth=0.5, fontsize=20)
    
50. x[~bool_grid_1d]=np.nan
    
51. y[~bool_grid_1d]=np.nan
    
52. m.scatter(x,y,s=10,color='red')
    
53. 
    
54. ###########################################################################
    
55. ###### Calculating area mean by np,nanmean#################################
    
56. ###########################################################################
    
57. z[~bool_grid_original_shape]=np.nan
    
58. print np.nanmean(z)/9.8
    
59. 
    
60. plt.show()
    

复制代码

点的位置的示意图如下：


最后得到的计算结果如下：

可见，北京市、天津市和辽宁省三地的平均地形海拔高度为236m。

脚本、海拔高度数据和地图边界文件打包上传了，可下载使用： area_mean.rar (626.66 KB, 下载次数: 66)

2020-10-11 16:12 上传

点击文件名下载附件

likailing2008

楼主厉害，写的帖子都是我想要的，赞赞赞

zane

有趣有趣 谢谢分享

Alliris

谢谢分享，又有新收获

qzf688018

path函数的使用，matlab里的inpolygon函数应该也能做到。

平流层的萝卜

qzf688018 发表于 2020-10-12 09:40
path函数的使用，matlab里的inpolygon函数应该也能做到。

matlab是肯定可以的，毕竟matlab的科学计算的功能比python强大很多

dongyi

棒棒哒，也是想要的

feiji158

大大，请教以下问题，1.此方法是否可以实现某个物理量等值线间的地理面积的计算？2.在用再分析资料时，往往某地海拔高度下仍有数值，以前的方法基本都是贴个狗皮膏（尤其是grads），很不美观，根据您的mask方法，似乎应该可以实现海拔高度下进行nan处理，这需要如何实现？望解答，谢谢！

平流层的萝卜

feiji158 发表于 2021-12-8 12:55
大大，请教以下问题，1.此方法是否可以实现某个物理量等值线间的地理面积的计算？2.在用再分析资料时，往往 ...

你好，第一个问题还真不了解。第二个问题，可以构建一个经纬度二维的mask_array的数组，判断一下，当某个经纬度的高度低于当地海拔时，设置mask_array该位置的元素为nan，否则为1，最后用mask_array与原始的数组相乘，nan乘出来就是nan，1乘出来就是原始值。

feiji158

平流层的萝卜 发表于 2021-12-09 09:42
你好，第一个问题还真不了解。第二个问题，可以构建一个经纬度二维的mask_array的数组，判断一下，当某个经纬度的高度低于当地海拔时，设置mask_array该位置的元素为nan，否则为1，最后用mask_array与原始的数组相乘，nan乘出来就是nan，1乘出来就是原始值。

厉害

                               
登录/注册后可看大图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

，谢谢大佬解答.

来自苹果APP客户端