1

我需要访问一些 grib 文件。我已经弄清楚了如何使用pygrib。但是，我想出如何做到这一点的唯一方法是非常缓慢。

我有 34 年的 3 小时数据，它们每年以约 36 个文件的形式组织起来（或多或少每 10 天一个）。总共大约 1000 个文件。

每个文件有大约 80 条“消息”（每天 8 个值，持续 10 天）。（它们是空间数据，因此它们具有 (x,y) 维度）。

要读取我写的所有数据：

grbfile = pygrib.index(filename, 'shortName', 'typeOfLevel', 'level') 
var1 = grbfile.select(typeOfLevel='pressureFromGroundLayer', level=180, shortName='unknown')
for it in np.arange(len(var1)):
    var_values, lat1, lon1 = var1[it].data()
    if (it==0):
        tot_var = np.expand_dims(var_values,axis=0)
    else:
        tot_var = np.append(tot_var, np.expand_dims(var_values,axis=0),axis=0)

并对 1000 个文件中的每一个重复此操作。

有更快的方法吗？像一次加载每个 grib 文件的所有 ~80 层？就像是：

var_values, lat1, lon1 = var1[:].data()
4

1 回答 1

1

如果我理解正确，您希望每个文件中所有 80 条消息的数据都堆叠在一个数组中。

我必须警告您，该数组会变得非常大，并且可能会导致 NumPyMemoryError根据您的网格大小等抛出一个（之前发生在我身上）。

话虽如此，您可以执行以下操作：

# substitute with a list of your file names
# glob is a builtin  library that can help accomplish this
files = list_of_files

grib = pygrib.open(files[0]) # start with the first one

# grib message numbering starts at 1
data, lats, lons = grib.message(1).data()

# while np.expand_dims works, the following is shorter
# syntax wise and will accomplish the same thing
data = data[None,...] # add an empty dimension as axis 0

for m in xrange(2, grib.messages + 1):
    data = np.vstack((data, grib.message(m).values[None,...]))

grib.close()  # good practice

# now data has all the values from each message in the first file stacked up 
# time to stack the rest on there
for file_ in files[1:]:  # all except the first file which we've done
    grib = pygrib.open(file_)
    for msg in grib:
       data = np.vstack((data, msg.values[None,...]))

   grib.close()
print data.shape # should be (80 * len(files), nlats, nlons)

这可能会让你加快速度。pygrib.open对象的作用类似于生成器，因此它们会在调用每个对象时将其传递给您，而不是像a 的方法pygrib.gribmessage那样构建它们的列表。如果您需要特定文件中的所有消息，那么这就是我访问它们的方式。select()pygrib.index

希望能帮助到你！

于 2016-11-14T17:10:13.830 回答