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\Orthobilder\20120502\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.

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\Orthobilder\20120502\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\Orthobilder\20120502\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................."

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.