Python/GDAL--Handling Big Rasters and avoid MemoryErrors?

Question

I am using python 2.7.10. I am trying to read a big Raster with gdal and numpy array. But it show In MemoreErrors when i use this: band. ReadAsArray()

Note: this error occurs when input raster is very big not for all of the rasters

Here is my Code:

 import sys, os, gdal, struct
from osgeo import osr
from gdalconst import *
import numpy
import time
import subprocess
startTime=time.time()
gdal.UseExceptions()
##OutCellSize = 1000
##OutFormat = "GTiff"
# Command line arguments
# change to known values if you prefer
InRas= r"F:\Geodaten\Eigene\Luftbild\Orthobilder20120502円\Aschhorn_mosaik.jpg"
src_ds=gdal.Open(InRas)
if src_ds is None:
 print 'Unable to open %s' % InRas
 sys.exit(1)
print "gdal can open the file"
srcband = src_ds.GetRasterBand(1)
try:
 bandArray = srcband.ReadAsArray()
 print bandArray
except:
 print "cannot convert raster to array.\n" "Need to compress raster"
# obtaining reference info from input raster
geotransform = src_ds.GetGeoTransform()
originX = geotransform[0]
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = src_ds.RasterXSize
rows = src_ds.RasterYSize ###we are here
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(src_ds.GetProjectionRef())
# creating raster in mepory as outraster from input
try:
 dst_ds = gdal.GetDriverByName("MEM").Create('', src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Byte)
 dst_ds.SetGeoTransform(geotransform)
 # Create for target raster the same projection as for the value raster
 dst_ds.SetProjection(raster_srs.ExportToWkt())
except:
 print "Unable to create in memory raster"

How could I solve this issue. Other wise I am not able to do further analysis.

• user30184: I would check the parameters in your geotransform = src_ds.GetGeoTransform()
• Shiuli Pervin: I have the following parameters in my geotransform: [518500.0, 0.1, 0.0, 5952000.0, 0.0, -0.1]
• user30184: Looks good. What values you have in "cols" and "rows"?
• Shiuli Pervin: cols:50000 and rows: 55000 file size: 1,04 GB
• user30184: In memory it will be uncompressed, three bands: 50*55*3=7.25 GB. Do you have that much free?

Answer 1

I have used this examples from Optimizing Python GDAL ReadAsArray here and use this for dealing with reading and writing raster and it has worked out. I would also like to to share my modifications:

import numpy
from numpy import zeros
from numpy import logical_and
from osgeo import gdal, osr
#import struct
'''Reading_Large_Raster_Blockwise'''
in_file = r"F:\Geodaten\Eigene\Luftbild\Orthobilder20120502円\Aschhorn_mosaik.jpg"
ds = gdal.Open(in_file)
#get projection reference
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(ds.GetProjectionRef())
geoT=ds.GetGeoTransform()
band = ds.GetRasterBand(1)
# Get "natural" block size, and total raster XY size. 
block_sizes = band.GetBlockSize()
x_block_size = block_sizes[0]
y_block_size = block_sizes[1]
xsize = band.XSize
ysize = band.YSize
print x_block_size, y_block_size
format = "GTiff"
driver = gdal.GetDriverByName( format )
dst_ds = driver.Create("original_blocks.tif", xsize, ysize, 1, band.DataType )
dst_ds.SetGeoTransform(geoT)
dst_ds.SetProjection(raster_srs.ExportToWkt())
print"working..................."
# Function to read the raster as arrays for the chosen block size.
def read_raster(x_block_size, y_block_size):
 raster = r"F:\Geodaten\Eigene\Luftbild\Orthobilder20120502円\Aschhorn_mosaik.jpg"
 ds = gdal.Open(raster)
 band = ds.GetRasterBand(1)
 xsize = band.XSize
 ysize = band.YSize
 blocks = 0
 for y in xrange(0, ysize, y_block_size):
 #print blocks
 if y + y_block_size < ysize:
 rows = y_block_size
 else:
 rows = ysize - y
 for x in xrange(0, xsize, x_block_size):
 if x + x_block_size < xsize:
 cols = x_block_size
 else:
 cols = xsize - x
 array = band.ReadAsArray(x, y, cols, rows)
 #print array.shape
 try:
 array[array>0]=1
 #print "we got them"
 except:
 print "could not find them"
 dst_ds.GetRasterBand(1).WriteArray(array, x, y)
 del array
 blocks += 1
#dst_ds.GetRasterBand(1).WriteArray(data)
read_raster(x_block_size, y_block_size)
band = None
ds = None
dst_ds = None
print"............... Done................."

• matandked: I've tried to apply your solution, but it doesn't work for me. I asked separate question: gis.stackexchange.com/questions/216698/…
• kory: his version only works for python 2
• kory: just add: you need to add "()" to prints, change the band.YSize to ds.RasterYSize, etc.

Answer 2

You could use the python module RIOS (Raster I/O Simplification). It reads an image into memory in a block of 400 x 400 x nbands and writes the output to a new image. It handles the opening of the existing image and the creation of the new image so the user can focus on the processing. RIOS is here: https://github.com/ubarsc/rios.