计算网格 netCDF 文件中选定区域中的变量平均值

假设我们有 TRMM 降水数据，每个文件代表每个月的数据。例如，文件夹中的文件是：

     3B42.1998.01.01.7A.nc,
     3B42.1998.02.01.7A.nc, 
     3B42.1998.03.01.7A.nc, 
     3B42.1998.04.01.7A.nc, 
     3B42.1998.05.01.7A.nc, 
     ......
     ......
     3B42.2010.11.01.7A.nc,         
     3B42.2010.12.01.7A.nc.

这些文件的尺寸如下：Xsize=1440、Ysize=400、Zsize=1、Tsize=1。经度设置为 0 到 360，纬度设置为 -50 到 50。 我想计算某个区域的降水量，比如说介于两者之间的降水量lon=98.5, lon=100 and lat=4, lat=6.5。这意味着，仅读取该区域中的变量 -：

-------------------- |lon:98.5 lat:6.5| | | |lat:4 lon:100 | ---------------------

我曾经在 GrADS（网格分析和显示系统）中这样做过。在 GrADS 中，可以这样做：（简化版本）

      yy=1998
      while yr < 2011
        'sdfopen f:\data\trmm\3B42.'yy'.12.01.7A.nc'
        'd aave(pcp,lon=98.5,lon=100.0,lat=4.0,lat=6.5)'
         res=subwrd(result,4)
         rec=write('d:\precip.sp.TRMM3B42.1.'yy'.csv',res,append)   
         yy = yy+1
      endwhile

我尝试在 Python 中做同样的事情，但出了问题。 经过一些建议后，我现在在这里：

     import csv
     import netCDF4 as nc 
     import numpy as np

     #calculating december only
     f = nc.MFDataset('d:/data/trmm/3B43.????.12.01.7A.nc')#maybe I shouldn't do MFDataset?
     pcpt = f.variables['pcp']
     lon = f.variables['longitude']
     lat = f.variables['latitude']
     # Determine which longitudes
     latidx1 = (lat >=4.0 ) & (lat <=6.5 ) 
     lonidx1 = (lon >=98.5 ) & (lon <=100.0 ) 

     rainf1 = pcpt[:]
     rainf1 = rainf1[:, latidx1][..., lonidx1]
     rainf_1 = rainf1

     with open('d:/trmmtest.csv', 'wb') as fp:
          a = csv.writer(fp)
          for i in rainf_1:
              a.writerow([i])

该脚本在 CSV 文件中生成（在我的例子中）15 个值的列表。 但是，当我尝试获取另一个区域的值并调整我认为必要的值时，可以说：

     latidx2 = (lat >=1.0 ) & (lat <=1.5 ) 
     lonidx2 = (lon >=102.75 ) & (lon <=103.25 ) 

     rainf2 = pcpt[:]
     rainf2 = rainf2[:, latidx2][..., lonidx2]
     rainf_2 = rainf2

我得到与第一个相同的值。

firstarea=[0.511935,1.0771,0.613548,1.48839,0.445161,1.39161,1.03548,0.452903, 3.07725,2.84613 0.701613,2.10581,2.47839,3.84097,2.41065,1.38387]

secondarea=[0.511935,1.0771,0.613548,1.48839,0.445161,1.39161,1.03548,0.452903, 3.07725,2.84613,0.701613,2.10581,2.47839,3.84097,2.41065,1.38387]

我确实在单独的脚本上进行了测试，它仍然给我相同的值。我确实检查了地图（之前构建的），这两个区域的值不同（12 月的平均值）。

知道为什么吗？还有其他优雅的写法吗？ 谢谢。

我只是想指出 Fir Nor 的解决方案是不正确的（更新：fir Nor 的帖子已删除，它之前建议基于使用 np.mean 的解决方案），因为您不能简单地使用算术平均值（np .mean）在常规纬度/经度网格上处理空间数据时，就像这里的情况一样由于网格单元尺寸发生变化当你向两极移动时。

以下是 python xarray 页面上的讨论，演示了如果不应用加权平均值会出现的差异。

我还制作了一个未装箱的气候关于此主题的 YouTube 视频解释为什么未加权平均值不正确以及如何使用 CDO 计算空间统计数据。

1.CDO解决方案：

最好不要担心这个并使用 CDO 进行操作：

cdo fldmean -sellonlatbox,98.5,100,4.5,6 3B42.1998.05.01.7A.nc boxav.nc

2.Python解决方案

如果你想用Python来做到这一点，你需要为你的子区域生成权重，可以根据你的解决方案（或使用xarray.where）提取权重。

如果您的纬度是一维，您可以使用以下命令将其转换为二维数组numpy.meshgrid

然后在二维数组上生成权重，并计算加权平均:

 weights = np.cos(np.deg2rad(lat2d))
 meanrain = numpy.average(pcpt, weights=weights)

使用 xarray 进行权重计算和错误诊断的另一个示例可以在我的答案中找到: