在网格化的netCDF-fi中计算一个选择区域内变量的平均值

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 and lat=4, lat=6.5之间。这意味着，要只读取此区域中的变量-：

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

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

^{pr2}$

我试着用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月的平均值）。在

知道为什么吗？还有其他优雅的写作方式吗？ 谢谢。在

Answer 1

我只想指出，Fir-Nor的解是不正确的，因为你不能简单地使用算术平均值(np平均值)当在常规的lat/lon网格上处理空间数据时，就像这里的情况一样，因为网格单元的大小随着你向两极移动而改变！在

最好不要担心这个，用CDO做手术：

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

Answer 2

过了一段时间，我设法再次研究这个问题，显然上面的方法几乎是正确的。经过几次调整，在单个数据文件上进行了测试，并与GrADS solution交叉检查，我得到了如下结果：

    f = nc.Dataset('~/data/TRMM3H/3B42.19980101.12.7A.nc')
    pcpt = f.variables['pcp'][:]
    lon = f.variables['longitude'][:]
    lat = f.variables['latitude'][:]

    #select two regions
    latidx1 = (lat >=4. ) & (lat <=6.5 ) 
    lonidx1 = (lon >=100.5 ) & (lon <=101.5 ) 
    latidx2 = (lat >=2.5 ) & (lat <=5.0 ) 
    lonidx2 = (lon >=101. ) & (lon <=102. ) 

    rainf = pcpt[:]
    #these basically listing the values in an array (2 in this case)
    rainf1 = rainf[:, latidx1][..., lonidx1]
    rainf2 = rainf[:, latidx2][..., lonidx2]
    rainf_1 = rainf1
    rainf_2 = rainf2

    #time to get the mean values
    print np.mean(rainf_1)
    print "............."
    print np.mean(rainf_2)
    print "............."

结果是：

^{pr2}$

用梯度法计算结果是一致的。在

建议后编辑：

    f = nc.Dataset('~/data/TRMM3H/3B42.19980101.12.7A.nc')
    pcpt = f.variables['pcp'][:]
    lon = f.variables['longitude'][:]
    lat = f.variables['latitude'][:]

    #select two regions
    latidx1 = (lat >=4. ) & (lat <=6.5 ) 
    lonidx1 = (lon >=100.5 ) & (lon <=101.5 ) 
    latidx2 = (lat >=2.5 ) & (lat <=5.0 ) 
    lonidx2 = (lon >=101. ) & (lon <=102. ) 

    #these basically listing the values in an array (2 in this case)
    rainf1 = pcpt[:, latidx1][..., lonidx1]
    rainf2 = pcpt[:, latidx2][..., lonidx2]
    rainf_1 = rainf1
    rainf_2 = rainf2

    #time to get the mean values
    print np.mean(rainf_1)
    print "............."
    print np.mean(rainf_2)
    print "............."

仍在测试文件中的原始文件。在