Trying to process rasters in python and received CPLE_OutOfMemoryError error

Question

I have a script that does some geoprocessing, that makes sense in a test run, but when I actually try to use it with my actual dataset, appears to totally overload my computer with memory needs, and I am very confused about about to actually fix this. I am using the rasterio and geocube packages.

I have a polygons shapefile here called polygons. This shapefile has three columns with unique values: Blue, Red, Green, and test. I want to, for each of these columns, rasterize each polygon, and then sum. So I would then have a raster of all Blue values summed, and then of all Red values summed, and of all Green values summed, and then of all test values summed. I then want to divide each of the summed color rasters by the summed test raster. to produce a final Blue, Red, and Green raster output. Here is my code:

list = ['Blue', 'Red', 'Green', `test']

for i in list:
    polygons = gpd.read_file('Colors/Polygons.shp')
    polygon_IDs = polygons['Polygon_ID'].tolist()

    # Make a GeoJSON string of the bounding box feature
    bbox = gpd.GeoSeries(box(*polygons.total_bounds), crs=polygons.crs)
    geom = bbox.__geo_interface__["features"][0]["geometry"]

    # Add CRS
    geom["crs"] = {"properties": {"name": f"EPSG:{polygons.crs.to_epsg()}"}}

    cubes = []
    for j in polygon_IDs:
        x = polygons.loc[polygons['Polygon_ID'] == j]
        vector_fn = x
        out_grid = make_geocube(
            vector_data=vector_fn,
            measurements=[i],
            resolution=(-25, 25),
            fill=0,
            geom=json.dumps(geom)
        )
        cubes.append(out_grid)
    out_grid = sum(cubes)
    out_grid[i].rio.to_raster(f"Colors/Output_Raster_{i}.tif", dtype=out_grid[i].dtype)

d = ['Blue', 'Red', 'Green']
test = 'Colors/Output_Raster_test.tif'
folder = Path('Colors')

src = rasterio.open(test) 

for f in folder.glob("*.tif"):
    if any(color in f.name for color in d):

        med = rasterio.open(f) 
        new_raster = med.read(1)/src.read(1)

        # write raster
        profile = src.profile
        profile.update(
            dtype=new_raster.dtype,
            count=1,
            compress='lzw')
        with rasterio.open(f"Colors/Output_Raster_{f.name}.tif", 'w', **profile) as dst:
            dst.write(new_raster, 1)

However, when I run this code I received this error: MemoryError: Unable to allocate 5.96 GiB for an array with shape (25631, 31233) and data type float64

I then delete files to make space and then receive this error: CPLE_OutOfMemoryError: memdataset.cpp, 1513: cannot allocate 1x6404264184 bytes

I am guessing this means that my code is just to intensive to run smoothly and just requires too much memory. But I am not sure which part is using up all of the memory. I am guessing it probably has to do with all of the raster summing, but I am not sure. Would anyone perhaps know how I might alter my code so that it will not require so much memory or what steps I might be going wrong in?

Answer 1

Here's a non geocube way using rasterio and fiona (could also use the higher level geopandas which uses fiona under the hood). This

• loops through small chunks (aka a rasterio.windows.Window) and the features in your polygon dataset,
• rasterizes the features just in the small chunk area,
• sums the chunks,
• writes each summed chunk to the appropriate part of the output raster

import fiona as fio
from itertools import product
import numpy as np
import rasterio as rio
from rasterio.features import rasterize
from rasterio import windows


#########################################
# Change these to suit your data
polys = 'your_data.shp'
poly_value_field = 'Value'
sum_raster = 'your_output.tif'
pixel_size = 25  # adjust appropriately if using a geographic CRS e.g. 0.0025
dtype = 'float32'
windows_shape = (1024, 1024)
all_touched = False  # only pixels with centre covered by polygon
#########################################


def get_windows(window_shape, image_shape):
    """Convenience generator for non-overlapping Windows"""
    win_rows, win_cols = window_shape
    img_rows, img_cols = image_shape
    offsets = product(range(0, img_cols, win_cols), range(0, img_rows, win_rows))
    image_window = windows.Window(col_off=0, row_off=0, width=img_cols, height=img_rows)

    for col_off, row_off in offsets:
        window = windows.Window(
            col_off=col_off,
            row_off=row_off,
            width=win_cols,
            height=win_rows)

        yield window.intersection(image_window)


with fio.open(polys) as features:
    crs = features.crs
    xmin, ymin, xmax, ymax = features.bounds
    transform = rio.Affine.from_gdal(xmin, pixel_size, 0, ymax, 0, -pixel_size)
    out_shape = (int((ymax - ymin)/pixel_size), int((xmax - xmin)/pixel_size))

    with rio.open(
            sum_raster, 'w',driver='GTiff',
            height=out_shape[0], width=out_shape[1], count=1,
            dtype=dtype, crs=crs, transform=transform,
            tiled=True, options=['COMPRESS=LZW']) as raster:

        for window in get_windows(windows_shape, out_shape):
            window_transform = windows.transform(window, transform)
            window_shape = (window.height, window.width)  # can be smaller than windows_shape at the edges
            window_data = np.zeros(window_shape)

            for feature in features:
                    value = feature['properties'][poly_value_field]
                    geom = feature['geometry']
                    d = rasterize(
                        [(geom, value)],
                        all_touched=all_touched, out_shape=window_shape, transform=window_transform
                    )
                    window_data += d  # sum values up

            raster.write(window_data, window=window, indexes=1)