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)


<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)


{'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:

population = rasterio.open(population_output)


{'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'}

  • 2
    rasterstats MemoryError issues are typically because of projection differences. Please include details (as an edit not a comment) about the vector and raster data - particularly extent coordinates of both and raster dimensions (rows/columns/cellsize) - and how you are changing the projection (in case you are setting a CRS instead of projecting to a CRS).
    – user2856
    Aug 29 at 22:49
  • 1
    Is your python 32 bit? If it is then it is restricted to 4 TiB of addressable memory, making a call for 1.61 TiB could be exceeding the maximum possible or that there is not a continuous block of 1.61 TiB to allocate. Aug 30 at 2:11
  • 2
    You aren't reprojecting the raster to new CRS, you are just changing the CRS definition. This is how to reproject a raster. But don't bother reprojecting the raster, just leave it in EPSG:4326 and if you need to, reproject the vectors to EPSG:4326.
    – user2856
    Aug 30 at 6:23
  • 2
    I don't believe you need an equal area CRS to do zonal stats. Particularly as by reprojecting your raster, you are resampling it and altering the values. You need an equal area CRS when you are using distance or area calcs.
    – user2856
    Aug 30 at 7:02
  • 2
    The error is because of the projection difference. Your raster is still in 4326 with geographic lon lat coordinates, but rasterio thinks it's in 54034 so tries to rasterize the vectors to an ridiculously big array that covers both extents.
    – user2856
    Aug 30 at 7:58

1 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.

Your Answer

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

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