将大规模栅格数据高效输入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数组时，我的电脑出现了内存错误。有没有其他方法可以将数据导入PyTables，或者有什么建议可以改善我的代码？

Answer 1

我没有地理信息的geotiff文件，所以我试着用一个普通的tif文件来操作。如果你要把数据写入pytables文件，可能需要省略形状中的3和切片的部分。基本上，我是循环处理这个数组，而不是一次性把所有数据都读到内存里。你需要调整n_chunks，这样每次读取的数据块大小就不会超过你的系统内存。

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)