I have two large datasets,
RP2. I would like to subtract one from the other, however the datasets are extremely large and so exit out with a
MemoryError when I try a direct read:
RP1 = rasterio.open(big_file_path) RP1 = RP1.read(1, masked=True) RP1 = RP1_dataset.read(1, masked=True)
File "rasterio_io.pyx", line 348, in
Therefore I have tried to resort to windowed reading and writing as explained by the
rasterio docs. I wrote the following script:
import numpy import rasterio image = numpy.ones((150, 250), dtype=rasterio.ubyte) * 127 RP1_filename = r'RP1.tif' RP2_filename = r'RP2.tif' RP2 = rasterio.open(RP2_filename) profile = RP2.profile with rasterio.open(RP1_filename) as RP1: for block_index, window in RP1.block_windows(1): RP1_block = RP1.read(window=window,masked = True) RP2_block = RP2.read(window=window, masked = True) result_block = RP2_block - RP1_block print(result_block.shape) RP2.close()
RP2 have the same window schema I took the liberty of indexing into
RP2 by the current
RP1 block index. This seems to work fine and the
result_block variable contains a masked array as expected.
The problem is I don't know how to collect all these arrays properly such that I can write a new output raster with the result values.
Does anyone know how to do that?