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

Question

假设我们有 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 和 lat=4, lat=6.5 之间。这意味着，仅读取该区域中的变量 -:

--------------------
|经度:98.5 纬度:6.5|
| |
|纬度:4 经度: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

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

第一个区域=[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]

第二个区域=[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]

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

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

Answer 1

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

Here is a discussion on the python xarray pages that demonstrates the differences that occur if you do not apply a weighted mean.

我还制作了一个未装箱的气候 youtube video on this topic解释为什么未加权均值不正确以及如何使用 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 将其转换为二维数组

然后在二维数组上生成权重，并计算weighted average :

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

Another example of the weights calculation using xarray and a diagnostic of the error is found my my answer here :