1

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.

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

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://bitbucket.org/chchrsc/rios/overview.

| improve this answer | |
2

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................."
| improve this answer | |

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.