Annonce
#1 Tue 28 March 2017 13:08
- alexelan
- Juste Inscrit !
- Date d'inscription: 3 Feb 2017
- Messages: 1
QGIS 2.14.6 python: MemoryError fusion de raster
Bonjour,
Je dois fusionner des 7 rasters en partant du moins précis au plus précis (les plus précis écrasant les valeurs des moins précis). Je les ai tous recalé à la même taille.
Pour faire cette opération j'ai un script python dont je ne comprend rien et que je n'ai pas écrit mais qui devrait marcher.
Le problème c'est que la personne qui m'a laissé ce script n'est plus joignable et que quand je le lance je me retrouve bloqué avec un Memory Error.
Voici le script:
Code:
---------------------------- DEBUT DU SCRIPT--------------------------
import os
import gdal
import numpy
import osr
import sys
def raster2array(raster):
band = raster.GetRasterBand(1)
return band.ReadAsArray()
def getOrigin(raster):
geotransform = raster.GetGeoTransform()
return [geotransform[0], geotransform[3]]
def getLayers():
return iface.legendInterface().layers()
def array2raster(rasterfn,newRasterfn,array, bitDepth = 7):
raster = gdal.Open(rasterfn)
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize
driver = gdal.GetDriverByName('HFA')
outRaster = driver.Create(newRasterfn, cols, rows, 1, bitDepth)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
gdal.UseExceptions()
# Get layers as they appear in the QGis Layers panel
layers = getLayers()
firstLayer = layers.pop()
firstRasterPath = firstLayer.dataProvider().dataSourceUri()
raster = gdal.Open(firstRasterPath)
print "Using "" + firstLayer.name() + "" as base layer"
origin = getOrigin(raster)
dest = raster2array(raster)
for layer in reversed(layers):
print "Processing "" + layer.name() + ""."
provider = layer.dataProvider()
path = provider.dataSourceUri()
raster = gdal.Open(path)
layerOrigin = getOrigin(raster)
layerData = raster2array(raster)
dest = numpy.where(numpy.logical_and(0 < layerData, layerDataHors ligne
