Python/GDAL--¿Manejar grandes rastreos y evitar MemoryErrors?

Question

Estoy usando python 2.7.10. Estoy tratando de leer un gran Raster con gdal y numpy array. Pero se muestra en MemoreErrors cuando uso esto: banda. ReadAsArray()

Nota: este error se produce cuando la trama de entrada es muy grande, no para todas las tramas

Aquí está mi código:

    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"

¿Cómo puedo resolver este problema? De lo contrario, no puedo hacer más análisis.

Answer 1

He utilizado estos ejemplos de Optimización de Python GDAL ReadAsArray aquí y utilizar esto para hacer frente a la lectura y escritura de raster y ha funcionado. También me gustaría compartir mis modificaciones:

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

Puedes utilizar el módulo de python RIOS (Raster I/O Simplification). Lee una imagen en la memoria en un bloque de 400 x 400 x nbands y escribe la salida en una nueva imagen. Maneja la apertura de la imagen existente y la creación de la nueva imagen para que el usuario pueda centrarse en el procesamiento. RIOS está aquí: https://bitbucket.org/chchrsc/rios/overview .