2

假设我们有 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 月的平均值）。

知道为什么吗？有没有其他优雅的方式来写这个？谢谢。

4

4 回答 4

3

我只想指出 Fir Nor 的解决方案是不正确的，因为在处理常规纬度/经度网格上的空间数据时，您不能简单地使用算术平均值 (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 数组

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

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

使用 xarray 进行权重计算和错误诊断的另一个示例是我的答案：

于 2018-02-22T09:00:33.227 回答
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 将计算每个网格单元面积加权的平均值。

于 2020-08-26T11:28:23.507 回答
1

过了一会儿，我又重新审视了这个问题，显然上面的方法几乎是正确的。经过一些调整，在单个数据文件上进行测试，并与 GrADS 解决方案进行交叉检查，我得到了这样的结果：

    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 "............."

这给了我这些结果：

    >>> execfile('find_percentile.py')
    0.7830327034
    .............
    1.56235361099
    .............

使用 GrADS 计算时的结果是相同的。

建议后编辑：

    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 "............."

回到最初的问题，在多个文件中执行此操作并在 txt/csv 文件中打印它仍在建设中（和测试）。

于 2014-09-22T18:24:48.860 回答
0

我相信它可以很容易地用easymore包来完成。

第一步是创建 shapefile。这可以是任何形状（例如点、盆地或矩形）。在您的情况下，它将是一个矩形 shapefile，其中一个形状定义了边界。这可以在 QGIS、ArcGIS 或 python 中完成：

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

接下来是调用easymore python包并将变量重新映射到感兴趣的shapefile，如下所示：

# 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- 使用这个包，您可以提供一个具有多个形状（感兴趣区域）的 shapefile，它可以一次性完成重新映射。例如，您可以改为提供世界各国的 shapefile。

2- 如果您的框小于网格（多边形）或点，则返回值将是小框或点所在的网格。

3- 重新映射和加权在等面积内完成，以计算 WGS84 中更高纬度的不同等面积网格。

4- 代码足够智能，因此您无需担心 0 到 360 lon 格式和 -180 到 180 lon 格式的目标 shapefile。例如，如果盒子是 NA，负 lon 值，shapefile 可以以负 lon 格式给出 -180 到 180，而 nc 文件具有非负 lon 值（0 到 360）。

包含更多示例的 GitHub 页面：

https://github.com/ShervanGharari/EASYMORE

于 2021-09-10T05:34:49.423 回答