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 (Grid Analysis and Display System) でこれを行っていました。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]

別のスクリプトでテストしましたが、それでも同じ値が得られます。マップ (以前に作成) を確認しましたが、値はこれら 2 つの地域で異なります (12 月の平均)。

理由はありますか？これを書く他のエレガントな方法はありますか？どうも。

4

4 に答える 4

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 に答える