0

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?

2
  • 1
    I'd explore using windowed read/writes: rasterio.readthedocs.io/en/latest/topics/…
    – mikewatt
    Commented Jan 5, 2022 at 23:39
  • Thanks for the tip. I tried the first part of my script, starting from list = ['Blue', 'Red', 'Green', test'] and down to out_grid[i].rio.to_raster(f"Colors/Output_Raster_{i}.tif", dtype=out_grid[i].dtype) just to see if those rasters actually would get made. They did not. After a good 10 or so minutes running, nothing happened, and eventually the memory just ran out again and I am left with no output rasters and this error message: CPLE_OutOfMemoryError: memdataset.cpp, 1513: cannot allocate 1x6404264184 bytes. I would like to use windows, but I'm not sure if they work with geocube. Commented Jan 6, 2022 at 2:38

1 Answer 1

3

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)
8
  • This works great! My output sum raster of all 500-some rasters was successfully produced! Thank you very much! Commented Jan 14, 2022 at 1:09
  • I do have a question though. I see at the top of the script you set both the vector and rasters to EPSG:5070. Does this mean that the output rasters will both match the input shapefile as EPSG:5070, such that when I open the input shapefile and the output rasters in QGIS they will directly overlap? I ask since I am not seeing shapefile and rasters are overlapping. Commented Jan 18, 2022 at 4:33
  • They do overlap for me, because the input shapefile is the same CRS, but you shouldn't use that CRS just because I did in my example, use whatever \coordinate system you want. See update.
    – user2856
    Commented Jan 18, 2022 at 4:57
  • My input shapefile is CRS: EPSG:2228 - NAD83 / California zone 4 (ftUS) - Projected and the output rasters are CRS: EPSG:5070 - NAD83 / Conus Albers - Projected, and they do not line up in QGIS. I tried assigning the rasters to the projection of the vector, but this is not working in QGIS. I am thinking that perhaps just restarting your script and setting the CRS to EPSG:2228 would work then. Also, in your updated script above for the line: pixel_size = 25 # adjust appropriately if using a geographic CRS e.g. 0.0025, how does pixel size determine the CRS? Commented Jan 18, 2022 at 5:21
  • That is why I said in the script comments # Change these to suit your data Just re-run the script (with my changes, copy the whole script, I made a change further down). Pixel size doesn't determine CRS. You choose an appropriate pixel size for your CRS.
    – user2856
    Commented Jan 18, 2022 at 5:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.