Python读取可变大小的多个NetCDF降雨文件

Question

我遇到的问题是澳大利亚气象局向我提供了降雨数据文件，其中包含所有活动仪表每 30 分钟记录一次的降雨记录。问题是 1 天有 48 个 30Minute 文件。我想创建特定仪表的时间序列。这意味着读取所有 48 个文件并搜索仪表 ID，确保在 1 30 分钟内仪表没有记录任何内容时它不会失败？ 这是文件格式的链接：

https://dl.dropboxusercontent.com/u/15223371/14/gauge_30min_20100214_000000.nc

https://dl.dropboxusercontent.com/u/15223371/14/gauge_30min_20100214_003000.nc

https://dl.dropboxusercontent.com/u/15223371/14/gauge_30min_20100214_010000.nc

这是我迄今为止尝试过的：

"""
This script is used to read a directory of raingauge data from a Data Directory





"""
from anuga.file.netcdf import NetCDFFile
from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
                            netcdf_float
import os
import glob
from easygui import *
import string
import numpy
"""
print 'Default file Extension...'
msg="Enter 3 letter extension."
title = "Enter the 3 letter file extension to search for in DIR "
default = "csv"
file_extension = enterbox(msg,title,default)
"""


print 'Present Directory Open...'
title = "Select Directory to Read Multiple rainfall .nc files"
msg = "This is a test of the diropenbox.\n\nPick the directory that you wish to open."
d = diropenbox(msg, title)
fromdir = d

filtered_list = glob.glob(os.path.join(fromdir, '*.nc'))
filtered_list.sort()

nf = len(filtered_list)
print nf

import numpy

rain = numpy.zeros(nf,'float')
t = numpy.arange(nf)

Stn_Loc_File='Station_Location.csv'
outfid = open(Stn_Loc_File, 'w')

prec = numpy.zeros((nf,1752),numpy.float)

gauge_id_list = ['570002','570021','570025','570028','570030','570032','570031','570035','570036',
                 '570047','570772','570781','570910','570903','570916','570931','570943','570965',
                 '570968','570983','570986','70214','70217','70349','70351']
"""
title = "Select Gauge to plot"
msg = "Select Gauge"
gauge_id = int(choicebox(msg=msg,title=title, choices=gauge_id_list))
"""
#for gauge_id in gauge_id_list:
#    gauge_id = int(gauge_id)
try:    

    for i, infile in enumerate(filtered_list):

        infilenet = NetCDFFile(infile, netcdf_mode_r)
        print infilenet.variables
        raw_input('Hold.... check variables...')
        stn_lats = infilenet.variables['latitude']
        stn_longs = infilenet.variables['longitude']
        stn_ids = infilenet.variables['station_id']
        stn_rain = infilenet.variables['precipitation']

        print stn_ids.shape
        #print stn_lats.shape
        #print stn_longs.shape
        #print infile.dimensions
        stn_ids = numpy.array(stn_ids)

        l_id = numpy.where(stn_ids == gauge_id)
        if stn_ids in gauge_id_list:
            try:
                l_id = l_id[0][0]
                rain[i] = stn_rain[l_id]
            except:
                rain[i] = numpy.nan
    print 'End for i...'            
    #print rain

    import pylab as pl

    pl.bar(t,rain)
    pl.title('Rain Gauge data')
    pl.xlabel('time steps')
    pl.ylabel('rainfall (mm)')
    pl.show()
except:
    pass 
raw_input('END....')

Answer 1

好的，你得到的数据格式比它需要的更复杂。他们可以很容易地将一整天的时间塞进一个 netCDF 文件中。实际上，解决此问题的一种选择是将所有文件合并为一个具有时间维度的文件，例如使用 NCO 命令行工具。

但这里有一个使用 scipy netcdf 模块的解决方案。我相信它已被弃用——我自己，我更喜欢 NetCDF4 库。主要做法是：用np.nan值预设你的输出数据结构；循环遍历您的输入文件并检索降水和站点 ID；对于您感兴趣的每个站点，检索索引，然后在该索引处进行降水；添加到输出结构。 （我没有做提取时间戳的工作——这取决于你。）

import glob
import numpy as np
from scipy.io import netcdf

# load data file names 
stationdata = glob.glob('gauge*.nc')
stationdata.sort()
# initialize np arrays of integer gauging station ids
gauge_id_list = ['570002','570021','570025','570028','570030','570032','570031','570035','570036',
                 '570047','570772','570781','570910','570903','570916','570931','570943','570965',
                 '570968','570983','570986','70214','70217','70349','70351']
gauge_ids = np.array(gauge_id_list).astype('int32')
ngauges = len(gauge_ids)
ntimesteps = 48
# initialize output dictionary
dtypes = zip(gauge_id_list, ['float32']*ngauges)
timeseries_per_station = np.empty((ntimesteps,))
timeseries_per_station.fill(np.nan)
timeseries_per_station = timeseries_per_station.astype(dtypes)

# Instead of using the index, you could extract the datetime stamp 
for timestep, datafile in enumerate(stationdata):
    data = netcdf.NetCDFFile(datafile, 'r')
    precip = data.variables['precip'].data
    stid = data.variables['stid'].data
    # create np array of indices of the gaugeid present in file
    idx = np.where(np.in1d(stid, gauge_ids))[0]
    for i in idx:
        timeseries_per_station[str(stid[i])][timestep] = precip[i]
    data.close()

np.set_printoptions(precision=1)
for gauge_id in gauge_id_list:
    print "Station %s:" % gauge_id
    print timeseries_per_station[gauge_id]

输出如下：

Station 570002:
[ 1.9  0.3  0.   nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan
  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan
  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan
  nan  nan  nan]
Station 570021:
[  0.   0.   0.  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan
  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan
  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan  nan
  nan  nan  nan]
...

（显然，只有三个文件。）

编辑： OP 指出，由于他的变量名称是“precipitation”和“station_id”，因此代码运行时不会没有错误。该代码在他发布的文件上为我运行。显然，他应该使用提供给他的文件中使用的任何变量名。由于它们似乎是供他使用的定制文件，因此可以想象作者在变量命名上可能不一致。