在网格化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.

这些文件的维度如下：X大小=1440，Y大小=400，Z大小=1，T大小=1。经度范围是0到360，纬度范围是-50到50。 我想计算某个区域的降水量，比如在lon=98.5, lon=100 和 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])

这个脚本生成了一个列表（在我的情况下）包含15个值，保存在CSV文件中。 但是当我尝试获取另一个区域的值，并调整我认为必要的部分，比如：

     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

我觉得用easymore这个包可以很简单地做到这一点。

第一步是创建一个形状文件。这个形状可以是任何形式的，比如点、子流域或者矩形。在你的情况下，它将是一个矩形的形状文件，里面有一个形状来定义边界。你可以在QGIS、ArcGIS或者用Python来完成这一步：

从边界框坐标列表创建形状文件

接下来就是调用easymore这个Python包，把变量映射到你感兴趣的形状文件上，操作非常简单，如下所示：

# loading EASYMORE
from easymore.easymore import easymore

# initializing EASYMORE object
esmr = easymore()

# specifying EASYMORE objects
# name of the case
esmr.case_name                = 'TRMM_3B43'              
# temporary path that the EASYMORE generated GIS files and remapped file will be saved
esmr.temp_dir                 = 'path/temporary/'
# name of target shapefile that the source netcdf files should be remapped to;
# it was created in the first step
esmr.target_shp               = 'path/target_shapefiles/box.shp'
# name of netCDF file(s); multiple files can be specified with *
esmr.source_nc                = ' d:/data/trmm/3B43*.nc'
# name of variables from source netCDF file(s) to be remapped
esmr.var_names                = ['pcp']
# name of variable longitude in source netCDF files
esmr.var_lon                  = 'longitude'
# name of variable latitude in source netCDF files
esmr.var_lat                  = 'latitude'
# name of variable time in source netCDF file; should be always time
esmr.var_time                 = 'time'
# location where the remapped netCDF, csv file will be saved
esmr.output_dir               = 'path/output/'
# if required that the remapped values to be saved as csv as well
esmr.save_csv                 = True

# execute EASYMORE nc remapper
esmr.nc_remapper()

这段代码会为每个原始的nc文件生成重新映射的nc文件和它的csv版本，存放在输出目录里。重新映射的文件将是你感兴趣的形状区域内的降水量的面积平均值，时间分辨率保持原样（比如说按天）。然后你可以轻松地将它们转换为按月的时间步长进行比较。

这个方法的优点有：

1- 使用这个包，你可以提供一个包含多个形状（感兴趣区域）的形状文件，它会一次性完成重新映射。比如，你可以直接提供世界各国的形状文件。

2- 如果你的框比网格（多边形）或点要小，返回的值将是这个小框或点所在的网格。

3- 重新映射和加权是按等面积进行的，以考虑在WGS84坐标系中高纬度地区不同的等面积网格。

4- 这段代码很聪明，所以你不需要担心经度格式的问题，比如0到360的经度格式和目标形状文件的-180到180的经度格式。例如，如果框在北美，负的经度值，形状文件可以用负经度格式-180到180，而nc文件则用非负经度值（0到360）。

更多示例可以在GitHub页面找到：

https://github.com/ShervanGharari/EASYMORE

Answer 2

如果你在使用Linux系统，可以通过nctoolkit来解决这个问题（nctoolkit.readthedocs.io/en/latest/）。下面的代码应该能完成所有的工作：

import nctoolkit as nc
ff = '~/data/TRMM3H/3B42.19980101.12.7A.nc'
data = nc.open_data(ff)
data.crop(lon = [98.5, 100], lat = [4, 6.5])
data.spatial_mean()

注意：这个方法使用了CDO作为后台工具，而spatial_mean会根据每个网格单元的面积来计算加权平均值。

Answer 3

我想指出，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）。

如果你的纬度是 1D 的话，可以使用 numpy.meshgrid 将其转换为 2D 数组。

然后在这个 2D 数组上生成权重，并计算 加权平均数：

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

这里有另一个使用 xarray 计算权重的例子，以及我在这里的回答中对错误的诊断。