1

I'm attempting to make a spatio-temporal statistics analysis on a few chunks of data using rasterstats, rasterio and pandas/geopandas on Python.

# create dataframes for areas of interest and specify coordinate systems

dornod_df = pd.DataFrame(mn_provinces.loc[1]).transpose()
dornod_polygon = gpd.GeoDataFrame(dornod_df, geometry="geometry")
dornod_polygon.crs = mn_provinces.crs
# creating a pandas dataframe to put the spatial data in it
data = pd.DataFrame('', columns=['Recorded_Month', 'Average_Fires_Dornod'], index=np.arange(1, 13))
# reading information from the rasters and extracting zonal stats

for rast in os.listdir(r'C:\Users\david\Desktop\2001'):
    if rast[-4:] == '.tif':
        rf = rasterio.open(r'C:\Users\david\Desktop\2001' + '\\' + rast)
        fire_array = rf.read(1)
        affine = rf.transform

        avg_fire_dornod = rasterstats.zonal_stats(dornod_polygon,
                                                  fire_array,
                                                  affine=affine,
                                                  stats=['mean'],
                                                  geojason_out=True)

        print(rast)

This returns an error:

Traceback (most recent call last):
  File "C:\anaconda3\envs\thesis_py38\lib\site-packages\IPython\core\interactiveshell.py", line 3437, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-f67c0b21d970>", line 1, in <module>
    runfile('D:/GIS/undergraduateThesis/programming/geodeltalabs_practice/spatial_temporal_statistics.py', wdir='D:/GIS/undergraduateThesis/programming/geodeltalabs_practice')
  File "C:\Program Files\JetBrains\PyCharm Community Edition 2020.3.5\plugins\python-ce\helpers\pydev\_pydev_bundle\pydev_umd.py", line 197, in runfile
    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script
  File "C:\Program Files\JetBrains\PyCharm Community Edition 2020.3.5\plugins\python-ce\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "D:/GIS/undergraduateThesis/programming/geodeltalabs_practice/spatial_temporal_statistics.py", line 28, in <module>
    avg_fire_dornod = rasterstats.zonal_stats(dornod_polygon,
  File "C:\anaconda3\envs\thesis_py38\lib\site-packages\rasterstats\main.py", line 31, in zonal_stats
    return list(gen_zonal_stats(*args, **kwargs))
  File "C:\anaconda3\envs\thesis_py38\lib\site-packages\rasterstats\main.py", line 156, in gen_zonal_stats
    fsrc = rast.read(bounds=geom_bounds)
  File "C:\anaconda3\envs\thesis_py38\lib\site-packages\rasterstats\io.py", line 305, in read
    new_array = boundless_array(
  File "C:\anaconda3\envs\thesis_py38\lib\site-packages\rasterstats\io.py", line 186, in boundless_array
    out = np.ones(shape=window_shape) * nodata
  File "C:\anaconda3\envs\thesis_py38\lib\site-packages\numpy\core\numeric.py", line 203, in ones
    a = empty(shape, dtype, order)
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 97.8 PiB for an array with shape (92978367, 148080780) and data type float64

The dataset I'm using is MODIS' 500 meter burnt area MCD64A1.

The raster data I'm working with is 16-bit signed integer, if that helps. I've looked on the internet, and many sources say to increase page file, but I'm not certain how I'll be able to allocate 97.8 PiB.

I have tried the same code on a bigger dataset, but it returned the exact same error: the PiB allocation requirement was the same 97.8 PiB. The size of the area of interest is approximately 120000 square kilometers.

How do I overcome this memory error?

EDIT:

The Python I'm using is the version 3.8.8 (64-bit) packaged by conda.

I have 8 GB of RAM.

gdalinfo of a sample data returns this:

Size is 1793, 910
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (112.056152343750000,50.286621093750000)
Pixel Size = (0.004394531250000,-0.004394531250000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 112.0561523,  50.2866211) (112d 3'22.15"E, 50d17'11.84"N)
Lower Left  ( 112.0561523,  46.2875977) (112d 3'22.15"E, 46d17'15.35"N)
Upper Right ( 119.9355469,  50.2866211) (119d56' 7.97"E, 50d17'11.84"N)
Lower Right ( 119.9355469,  46.2875977) (119d56' 7.97"E, 46d17'15.35"N)
Center      ( 115.9958496,  48.2871094) (115d59'45.06"E, 48d17'13.59"N)
Band 1 Block=1793x2 Type=Int16, ColorInterp=Gray
  NoData Value=32767

ogrinfo of the shapefile representing the area of interest is as follows:

Layer name: dornod1

Metadata:

DBF_DATE_LAST_UPDATE=2021-04-10

Geometry: Polygon

Feature Count: 1

Extent: (1445925.343002, 5277981.646306) - (2096670.952278, 5686577.978462)

Layer SRS WKT:

PROJCRS["WGS 84 / UTM zone 47N",

BASEGEOGCRS["WGS 84",

DATUM["World Geodetic System 1984",

ELLIPSOID["WGS 84",6378137,298.257223563,

LENGTHUNIT["metre",1]]],

PRIMEM["Greenwich",0,

ANGLEUNIT["degree",0.0174532925199433]],

ID["EPSG",4326]],

CONVERSION["UTM zone 47N",

METHOD["Transverse Mercator",

ID["EPSG",9807]],

PARAMETER["Latitude of natural origin",0,

ANGLEUNIT["Degree",0.0174532925199433],

ID["EPSG",8801]],

PARAMETER["Longitude of natural origin",99,

ANGLEUNIT["Degree",0.0174532925199433],

ID["EPSG",8802]],

PARAMETER["Scale factor at natural origin",0.9996,

SCALEUNIT["unity",1],

ID["EPSG",8805]],

PARAMETER["False easting",500000,

LENGTHUNIT["metre",1],

ID["EPSG",8806]],

PARAMETER["False northing",0,

LENGTHUNIT["metre",1],

ID["EPSG",8807]]],

CS[Cartesian,2],

AXIS["(E)",east,

ORDER[1],

LENGTHUNIT["metre",1]],

AXIS["(N)",north,

ORDER[2],

LENGTHUNIT["metre",1]],

ID["EPSG",32647]]

Data axis to CRS axis mapping: 1,2

ID: Integer64 (10.0)

OBJECTID: Integer64 (10.0)

AREA_M2: Real (19.8)

LENGTH_M: Real (19.8)

TYPE: Integer (5.0)

NAME: String (50.0)

NAME_L1: String (50.0)

PARENT_ID: Integer64 (10.0)

UNIT_STATU: Integer (5.0)

SHAPE_AREA: Real (19.11)

SHAPE_LEN: Real (19.11)
3
  • Is it python 32 or 64 bit?
    – Taras
    Apr 10, 2021 at 7:58
  • Okay, I have included the output of gdalinfo and ogrinfo in the question.
    – sokhymi
    Apr 10, 2021 at 8:11
  • Thank you! Learned something new :)
    – sokhymi
    Apr 10, 2021 at 8:49

1 Answer 1

5

Your vector is in a different CRS to your raster so rasterstats is trying to create an enormous (97 petabyte) raster.

Either add a LOT more RAM :) or reproject your vector to match your raster:

dornod_df.to_crs(4326)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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