Merging rasters with rasterio in blocks to avoid MemoryError?

Question

I am trying to merge two rasters with rasterio.

However, when I run dest_array, out_transform = merge(sources), I receive the following error:

MemoryError: Unable to allocate array with shape (1, 56160, 85196) and data type float64

I know rasterio can read and write files in blocks.

Unfortunately, I cannot seem to find a way to do it when calling merge.

Is there a way to use blocks with merge or another way I can merge rasters without using the merge method?

I know how to read and write data in blocks, but I have no idea how to do that when using the merge method.

Here is the code I am using until the line that calls merge and that gives me the error:

import rasterio
from rasterio.merge import merge
import glob
import os

def merge_raster_in_folder(src_dir, dest_dir):
    """Merge all the raster (TIFF extension) within the given directory"""

    # create list of raster files in source folder
    rasters = []
    os.chdir(src_dir)
    for file in glob.glob("*.tif"):
        rasters.append(file)

    # create list of raster objects from the list of raster filename
    sources = [rasterio.open(raster) for raster in rasters]

    # create array representing all source rasters mosaicked together
    dest_array, out_transform = merge(sources)

Answer 1

I finally solved the problem by using numpy.memmap to

create a memory-map to an array stored in a binary file on disk

and then processing the input rasters in windows and blocks. It might be slower and but it works and I'm happy with the result (need to thank user @Thomas that helped me in some steps).

The code I am using is taken and modified from the source code of rasterio merge.py, until the part where it creates a "huge" (in my case) numpy array of zeros that was hitting the MemoryError in my case.

Here is the final code:

from tempfile import mkdtemp
import rasterio
from rasterio import Affine
from rasterio import windows
import math
import numpy as np
import os

INPUT_FILES = [r'path/to/raster1.tif', r'path/to/rasterN.tif']

sources = [rasterio.open(raster) for raster in INPUT_FILES]

memmap_file = os.path.join(mkdtemp(), 'test.mymemmap')

# adapted from https://github.com/mapbox/rasterio/blob/master/rasterio/merge.py
first = sources[0]
first_res = first.res
dtype = first.dtypes[0]
# Determine output band count
output_count = first.count

# Extent of all inputs
# scan input files
xs = []
ys = []
for src in sources:
    left, bottom, right, top = src.bounds
    xs.extend([left, right])
    ys.extend([bottom, top])
dst_w, dst_s, dst_e, dst_n = min(xs), min(ys), max(xs), max(ys)

out_transform = Affine.translation(dst_w, dst_n)

# Resolution/pixel size
res = first_res
out_transform *= Affine.scale(res[0], -res[1])

# Compute output array shape. We guarantee it will cover the output
# bounds completely
output_width = int(math.ceil((dst_e - dst_w) / res[0]))
output_height = int(math.ceil((dst_n - dst_s) / res[1]))

# Adjust bounds to fit
dst_e, dst_s = out_transform * (output_width, output_height)

# create destination array
# destination array shape
shape = (output_height, output_width)
# dest = np.zeros((output_count, output_height, output_width), dtype=dtype)
# Using numpy.memmap to create arrays directly mapped into a file
dest_array = np.memmap(memmap_file, dtype=dtype,
                       mode='w+', shape=shape)

dest_profile = {
    "driver": 'GTiff',
    "height": dest_array.shape[0],
    "width": dest_array.shape[1],
    "count": output_count,
    "dtype": dest_array.dtype,
    "crs": '+proj=latlong',
    "transform": out_transform
}

# open output file in write/read mode and fill with destination mosaick array
with rasterio.open(
    os.path.join(r'out/dir', 'test.tif'),
    'w+',
    **dest_profile
) as mosaic_raster:
    for src in sources:
        for ji, src_window in src.block_windows(1):

            print(ji)
            r = src.read(1, window=src_window)
            # store raster nodata value
            nodata = src.nodatavals[0]
            # replace zeros with nan
            r[r == nodata] = np.nan
            # convert relative input window location to relative output window location
            # using real world coordinates (bounds)
            src_bounds = windows.bounds(
                src_window, transform=src.profile["transform"])
            dst_window = windows.from_bounds(
                *src_bounds, transform=mosaic_raster.profile["transform"])

            # round the values of dest_window as they can be float
            dst_window = windows.Window(round(dst_window.col_off), round(
                dst_window.row_off), round(dst_window.width), round(dst_window.height))
            # before writing the window, replace source nodata with dest nodata as it can already have been written (e.g. another adjacent country)
            # https://stackoverflow.com/a/43590909/1979665
            dest_pre = mosaic_raster.read(1, window=dst_window)
            mask = (np.isnan(r))
            r_mod = np.copy(r)
            r_mod[mask] = dest_pre[mask]
            mosaic_raster.write(r_mod, 1, window=dst_window)

os.remove(memmap_file)

Answer 2

I faced the same issue. The easiest solution, in this case, is not using rasterio. Instead, using gdal will be much easier, without any need to modify the source code. Follow the second step of this video carefully.

Basically, they are creating a virtual raster XML file using gdal.BuildVRT() function of gdal. The virtual raster XML stores all the individual rasters that are to be merged. And the gdal.Translate() function is used to convert the XML to a GTiff. The BuildVRT() function bypasses the use of numpy.zeros() function in rasterio.merge() which is the source of memory error in the case of large datasets.