python - 将大型栅格数据输入 PyTables 的有效方法

Question

我正在寻找将 20GB 大小的栅格数据文件 (GeoTiff) 提供给 PyTables 以进一步进行核心计算的有效方法。

目前我正在使用 Gdal 将其读取为 numpy 数组，并将 numpy 数组写入pytables 使用下面的代码:

import gdal, numpy as np, tables as tb

inraster = gdal.Open('infile.tif').ReadAsArray().astype(np.float32)
f = tb.openFile('myhdf.h5','w')
dataset = f.createCArray(f.root, 'mydata', atom=tb.Float32Atom(),shape=np.shape(inraster)
dataset[:] = inraster
dataset.flush()
dataset.close()
f.close()
inraster = None

不幸的是，由于我的输入文件非常大，在将其读取为 numpy 错误时，我的 PC 显示内存错误。有没有其他方法可以将数据输入 PyTables 或有任何改进我的代码的建议？

Answer 1

我没有 geotiff 文件，所以我摆弄了一个普通的 tif 文件。如果将数据写入 pytables 文件，您可能必须省略形状中的 3 和写入中的切片。本质上，我遍历数组而不是一次将所有内容读入内存。您必须调整 n_chunks 以便一次性读取的 block 大小不会超过您的系统内存。

ds=gdal.Open('infile.tif')

x_total,y_total=ds.RasterXSize,ds.RasterYSize

n_chunks=100

f = tb.openFile('myhdf.h5','w')
dataset = f.createCArray(f.root, 'mydata', atom=tb.Float32Atom(),shape=(3,y_total,x_total)


#prepare the chunk indices
x_offsets=linspace(0,x_total,n_chunks).astype(int)
x_offsets=zip(x_offsets[:-1],x_offsets[1:])
y_offsets=linspace(0,y_total,n_chunks).astype(int)
y_offsets=zip(y_offsets[:-1],y_offsets[1:])

for x1,x2 in x_offsets:
    for y1,y2 in y_offsets:
        dataset[:,y1:y2,x1:x2]=ds.ReadAsArray(xoff=x1,yoff=y1,xsize=x2-x1, ysize=y2-y1)