Memory error: unable to allocate using zonal statistics

Question

I am doing a zonal statistics between a shapefile of cells around the African continent and a raster of population around the world. Both of them are in an equal area projection: World Cylindrical Equal area projection, CRS: ESRI:54034.

I manage python in a jupyter notebook which I access through miniconda prompt in Windows.

I am doing the following:

import numpy as np
import pandas as pd
import geopandas as gpd
from pyproj import CRS
import rasterio
from rasterstats import zonal_stats
from rasterio.mask import mask

Africa = gpd.read_file(data_directory + 'GridCells/ccode_GID.shp')
population = rasterio.open(pop_ag_1990)

#I do the zonal stats in the following way:

population_array, population_transform = mask(population, shapes=Africa.geometry, crop=True, nodata=np.nan)

zs1 = zonal_stats(Africa, population_array[0], affine=population_transform, stats=['sum'], nodata=np.nan)
Africa['population'] = [x['sum'] for x in zs1]

However, I received the following error:

MemoryError: Unable to allocate 1.61 TiB for an array with shape (662407, 667918) and data type float32

It is a bit strange because when I change the projection to WGS 84 EPSG:4326, the error does not appear, it works well. I have tried different solutions as restoring the kernel, closed and open the notebook, removing several files, etc. One solution that I have not tried is to increase the memory of the notebook changing the configuration as explained here: https://stackoverflow.com/questions/57948003/how-to-increase-jupyter-notebook-memory-limit. However, I work through a virtual environment, and the command "jupyter notebook" does not work for me.

Why my code does not work when I am using the cylindrical area projection? And how can I solve it?

[Edit: I add more information according to the comments I have received]

This is how I convert my shapefile in an equal area projection:

wcea_crs = CRS.from_string('esri:54034')

africa = africa.to_crs(crs=wcea_crs)

Africa.crs 

<Derived Projected CRS: ESRI:54034>
Name: World_Cylindrical_Equal_Area
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Coordinate Operation:
- name: World_Cylindrical_Equal_Area
- method: Lambert Cylindrical Equal Area
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

This is how I convert my raster in an equal area projection:

population = rasterio.open(population)

population.profile

{'driver': 'AAIGrid', 'dtype': 'float32', 'nodata': -9999.0, 'width': 2160, 'height': 1920, 'count': 1, 'crs': None, 'transform': Affine(0.0416666666667, 0.0, -30.0,
       0.0, -0.0416666666667, 40.00000000006399), 'tiled': False}

profile = population.profile
profile.update(driver='GTiff', crs=wcea_crs)
population_array = population.read()

population_output = "{}/population.tif".format(junkpath2)

with rasterio.open(population_output, 'w', **profile) as dst:
    dst.write(population_array) 

population = rasterio.open(population_output)

population.profile

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -9999.0, 'width': 2160, 'height': 1920, 'count': 1, 'crs': CRS.from_wkt('PROJCS["World_Cylindrical_Equal_Area",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Cylindrical_Equal_Area"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'), 'transform': Affine(0.0416666666667, 0.0, -30.0,
       0.0, -0.0416666666667, 40.00000000006399), 'tiled': False, 'interleave': 'band'}

Answer 1

rasterstats MemoryError issues are typically because of projection differences and in this case there is definitely a projection difference.

profile.update(driver='GTiff', crs=wcea_crs)

By doing the above, you are changing the projection definition, but not actually reprojecting the data so the coordinates do not change. It's like changing the file extension (projection definition) of a Word doc from .docx to .pdf and expecting Acrobat to read the file. It won't because it's a Word doc with the wrong file extension not a PDF, you need to actually export/save (reproject) your Word doc to a different format.

Your raster is still in EPSG:4326 with geographic lon lat coordinates, but rasterio thinks it's in 54034 so tries to rasterize the vectors to a ridiculously big array (662407, 667918) that matches the rasters pixel resolution. The size required for this array is (662407x667918x32)/8/1024^3 = 1648 Gb or 1.6 TB (terabytes!)* and you're very unlikely to have that much RAM available.

* Actually GiB/TiB but lets not get pedantic :)

If you want to reproject your raster, you can follow the rasterio reprojection example, but note that reprojecting will resample the raster and introduce error so it's best to avoid that if you don't really need to reproject.

However, since both your raster and vector are in EPSG:4326 and you are not doing distance or area calculations that require an equal area projection, there's no need to reproject anyway.