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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 27191|回复: 15

[经验总结] Python解决任意区域的格点包含和面积平均问题

[复制链接]

新浪微博达人勋

发表于 2020-10-11 16:19:09 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 平流层的萝卜 于 2020-10-15 10:40 编辑

气象中应该会经常遇到算一个区域内的平均值问题。有时候区域是矩形,则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. ncfile=nc.Dataset(r'surface_z_0.25_46_36_115_128.nc')
  9. lons=ncfile.variables['longitude'][:]
  10. lats=ncfile.variables['latitude'][:]
  11. z=ncfile.variables['z'][0,:,:]
  12. xx,yy=np.meshgrid(lons,lats)
  13. x,y=xx.flatten(),yy.flatten()
  14. points=np.vstack((x,y)).T

  15. # print yy
  16. # exit()

  17. sf=shapefile.Reader(r'bou2_4p',' ')
  18. vertices = []
  19. codes = []
  20. for shape_rec in sf.shapeRecords():
  21.     # print record
  22.     if shape_rec.record[4] in [210000,110000,120000]:   ###辽宁、北京和天津的行政编码,百度可查
  23.         pts = shape_rec.shape.points
  24.         prt = list(shape_rec.shape.parts) + [len(pts)]
  25.         for i in range(len(prt) - 1):
  26.             for j in range(prt[i], prt[i+1]):
  27.                 vertices.append((pts[j][0], pts[j][1]))
  28.             codes += [Path.MOVETO]
  29.             codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
  30.             codes += [Path.CLOSEPOLY]
  31.         clip = Path(vertices, codes)

  32. polygon_of_Liaoning_beijing_tianjin= clip
  33. # exit()

  34. bool_grid_1d = polygon_of_Liaoning_beijing_tianjin.contains_points(points)
  35. bool_grid_original_shape=bool_grid_1d.reshape(xx.shape)


  36. ###########################################################################
  37. ###### Highlighting grids within the area using scatterplot################
  38. ###########################################################################
  39. m = Basemap(llcrnrlat=lats[-1], llcrnrlon=lons[0], urcrnrlat=lats[0], urcrnrlon=lons[-1])
  40. m.readshapefile(r'E:\shpmap\bou2_4p', ' ', color='grey', linewidth=2)
  41. m.drawparallels(np.linspace(lats[0],lats[-1],10), labels=[1, 0, 0, 0], linewidth=0.5, fontsize=20)
  42. m.drawmeridians(np.linspace(lons[0],lons[-1],10), labels=[0, 0, 0, 1], linewidth=0.5, fontsize=20)
  43. x[~bool_grid_1d]=np.nan
  44. y[~bool_grid_1d]=np.nan
  45. m.scatter(x,y,s=10,color='red')

  46. ###########################################################################
  47. ###### Calculating area mean by np,nanmean#################################
  48. ###########################################################################
  49. z[~bool_grid_original_shape]=np.nan
  50. print np.nanmean(z)/9.8

  51. plt.show()
复制代码



点的位置的示意图如下:
Figure_1.png

最后得到的计算结果如下:
QQ截图20201011155857.png
可见,北京市、天津市和辽宁省三地的平均地形海拔高度为236m。

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

评分

参与人数 1金钱 +10 收起 理由
wet510 + 10 提到了西藏的边界,楼主很是细心,赞

查看全部评分

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

新浪微博达人勋

发表于 2020-10-11 16:28:15 | 显示全部楼层
楼主厉害,写的帖子都是我想要的,赞赞赞
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-10-12 07:42:07 | 显示全部楼层
有趣有趣 谢谢分享
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-10-12 09:13:00 | 显示全部楼层
谢谢分享,又有新收获
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-10-12 09:40:08 | 显示全部楼层
path函数的使用,matlab里的inpolygon函数应该也能做到。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-10-12 16:44:38 | 显示全部楼层
qzf688018 发表于 2020-10-12 09:40
path函数的使用,matlab里的inpolygon函数应该也能做到。

matlab是肯定可以的,毕竟matlab的科学计算的功能比python强大很多
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-3-31 15:09:22 | 显示全部楼层
棒棒哒,也是想要的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-12-8 12:55:15 | 显示全部楼层
大大,请教以下问题,1.此方法是否可以实现某个物理量等值线间的地理面积的计算?2.在用再分析资料时,往往某地海拔高度下仍有数值,以前的方法基本都是贴个狗皮膏(尤其是grads),很不美观,根据您的mask方法,似乎应该可以实现海拔高度下进行nan处理,这需要如何实现?望解答,谢谢!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-12-9 09:42:32 | 显示全部楼层
feiji158 发表于 2021-12-8 12:55
大大,请教以下问题,1.此方法是否可以实现某个物理量等值线间的地理面积的计算?2.在用再分析资料时,往往 ...

你好,第一个问题还真不了解。第二个问题,可以构建一个经纬度二维的mask_array的数组,判断一下,当某个经纬度的高度低于当地海拔时,设置mask_array该位置的元素为nan,否则为1,最后用mask_array与原始的数组相乘,nan乘出来就是nan,1乘出来就是原始值。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

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

厉害

                               
登录/注册后可看大图

                               
登录/注册后可看大图

                               
登录/注册后可看大图
,谢谢大佬解答.
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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