Le portail francophone de la géomatique

Bonjour,

y a t-il une raison sur le fait que vous passiez par le module gdal pour fusionner vos raster ?
Le plus simple vu que vous passez par la console python serait de passer par le module processing
et par l'algorithme 'gdalogr:merge'.
Si l'outil fusionner fonctionne je vous propose plutôt ce code :

import processing
#pour avoir de l'aide sur l'algo vous tapez 
#dans la console processing.alghelp('gdalorg:merge') 
#afin de voir les options qui vous conviennent
#paramètre de départ le chemin des rasters séparés par un ";"

input_layers = ';'.join([layer.source() for layer
    in processing.getAllLayers()])
newRasterfn = '/Users/sbooob/Desktop/QGis/SIG/out.tif'
processing.runalg('gdalogr:merge', input_rasters, 
    False, False,4,newRasterfn)

Dernière modification par lejedi76 (Tue 19 July 2016 13:46)

Il peut être , bien entendu, utile de poser la même question sur plusieurs forums, mais la politesse voudrait  qu'on le signale et qu'on assume en coordonnant les réponses.

http://gis.stackexchange.com/questions/ … ing-python
http://www.forumsig.org/showthread.php/ … post344283

Dernière modification par gene (Tue 19 July 2016 16:57)

Bonjour,

Pourquoi ne pas tout simplement utiliser une calculatrice raster... Je ne crois pas que gdal merge soit la solution c'est un outil dédié à la réalisation de mosaic (http://www.gdal.org/gdal_merge.html) et ce n'est clairement pas ce que vous essayez de faire.

Quelques remarques concernant votre script, il faut noté que la fonction ReadAsArray() retourne un tableau numpy et que dans un tableau numpy en 2D la première dimension représente toujours les lignes et la seconde les colonnes. On accède donc aux valeurs avec l'index y en 1er

data[row, column] == data[y, x]

Votre double itération devrait donc plutôt ressembler à ça:

ny, nx = layerData.shape
for x in range(nx):
    for y in range(ny):
        dest[y][x] = layerData[y][x]

Néanmoins affecter les valeurs une à une est une assez mauvaise pratique car verbeux et lent, il est préférable de profiter des fonctionnalités de numpy avec notamment la possibilité de tout simplement additionner des tableaux. Ainsi,

dest += layerData

remplacera avantageusement votre double itération.

Quoi qu'il en soit je vous engage à explorer l'utilisation des calculatrices raster. Il y en a plusieurs : celle de QGIS accessible depuis le menu raster et celles de SAGA, GDAL et GRASS accessible depuis processing.

Dernière modification par dominique.lys (Wed 20 July 2016 11:50)

Bonjour a vous deux!

Merci de vos réponse et pas de souci Lejedi76

Juste pour vous expliquer le contexte. Je suis en stage et je dois faire du SIG à un niveau qui me dépasse de beaucoup. Mon frère qui est développeur web et m'aide pour tout ce qui est "langage informatique" et python. Mais à nous deux on ne fait pas un très bon SIGiste, loin de la...

Mon frère viens de se lancer sur ce que tu nous a dit dominique.lys et le resultat est le même on a toujours 100m de decaalge verticale

Par contre pour la calculatrice raster c'est que je veux utiliser depuis le début mais on a pas trouvé de documentation précise sur comment s'en servir du coup on l'a un peu délaissé mais la on but et on risque de devoir s'en servir pour de bon!

En gros j'ai des 9 couches raster dont une qui sera à la base et qui va être compléter par les 8 autres. Je veux que les valeurs des pixels des couches que je vais ajouter se substituent à celle d'en dessous.
Je veux donc qu'il me remplace la valeurs du pixel de la couche de base par celle d'ajout. Si le pixel de la couche d'ajout est =0 alors il faut garder la valeur de la couche initial mais si il est différent de 0 alors il faut remplacer la valeur du pixel de la couche de base par celui de la couche d'ajout.
Donc soit je peux faire cela 8 fois et reprendre le fichier raster créer en sortie soit apparemment comme dit lejedi76 on peut le faire jusqu'à 6 paramètres ce qui pourrait nous faire gagner du temps.

Savez-vous comment le formuler pour que la calculatrice raster tourne correctement,  avec l'une ou l'autre des méthodes je suis preneur!

Merci de votre contribution!

Tristan et Yo

Dernière modification par Tris (Wed 20 July 2016 17:54)

J'ai essayé et ça ne marche pas, cela additionne les valeurs au leu de les remplacer, j'aurais aimé que ce soit aussi simple...

Prenons un exemple simplifié. raster1 est la couche inférieure, raster2 est la couche supérieure.
Je souhaite remplacer les pixels de raster1 par ceux de raster2, sauf la valeur du pixel de raster2 est égale à 0, auquel cas je veux conserver la valeur de raster 1 :

raster1 = [[ 0, 0, 1 ],
           [ 0, 1, 0 ],
           [ 1, 0, 0 ]]
           
raster2 = [[ 1, 0, 0 ],
           [ 0, 1, 0 ],
           [ 0, 0, 1 ]]

Je souhaite obtenir le rasterFusionné suivant :

rasterFusionné = [[ 1, 0, 1 ],
                  [ 0, 1, 0 ],
                  [ 1, 0, 1 ]]

La solution que vous proposez ne produit pas le même raster :

rasterFusionné = [[ 1, 0, 1 ],
                  [ 0, 2, 0 ],
                  [ 1, 0, 1 ]]

C'est pour effectuer un traitement différent si la valeur du pixel est égale à 0 que je me suis tourné vers Python. Peut-être qu'une des calculatrices raster peut aussi me permettre d'y arriver mais la formule doit être un peu plus compliquée qu'une simple addition.

Je ne pense pas que cela produise l'effet recherché.

Prenons un exemple simplifié. raster1 est la couche inférieure, raster2 est la couche supérieure.
Je souhaite remplacer les pixels de raster1 par ceux de raster2, sauf la valeur du pixel de raster2 est égale à 0, auquel cas je veux conserver la valeur de raster 1 :

raster1 = [[ 0, 0, 1 ],
           [ 0, 1, 0 ],
           [ 1, 0, 0 ]]
           
raster2 = [[ 1, 0, 0 ],
           [ 0, 1, 0 ],
           [ 0, 0, 1 ]]

Je souhaite obtenir le rasterFusionné suivant :

rasterFusionné = [[ 1, 0, 1 ],
                  [ 0, 1, 0 ],
                  [ 1, 0, 1 ]]

La solution que vous proposez ne produit pas le même raster :

rasterFusionné = [[ 1, 0, 1 ],
                  [ 0, 2, 0 ],
                  [ 1, 0, 1 ]]

C'est pour effectuer un traitement différent si la valeur du pixel est égale à 0 que je me suis tourné vers Python. Peut-être qu'une des calculatrices raster peut aussi me permettre d'y arriver mais la formule est un peu plus compliquée qu'une simple addition du coup, dommage... Mais merci quand même!

Problème résolu.

Le code de mon premier post était correct bien que loin d'être optimisé.

Le problème venait en fait du fait que deux de mes rasters n'avaient pas la même origine que les autres (les coordonnées du point supérieur gauche étaient différentes pour ces deux rasters et pour les 7 autres).

Mon script Pyhton ne prend pas en compte les coordonnées géographiques et ne manipule que les pixels des rasters. Il est donc important que les rasters aient bien tous la même origine.

Après avoir pris soin de définir une origine commune à tous mes rasters, et optimisé le script en prenant en compte les suggestions de dominique.lys, j'ai relancé le script. Le gain de performance est impressionnant, l'opération qui prenait plus de 45 minutes s'effectue désormais en moins d'une minute !

Au cas où ça puisse aider quelqu'un, voici comment j'ai défini une origine commune à mes rasters (cela créé une copie de chaque raster) :

* Dans la liste des couches, faire un clic droit sur un raster et cliquer sur "Enregistrer sous...".
* Dans la boîte de dialogue qui s'ouvre, sélectionner où enregistrer le fichier en cliquant sur "Parcourir".
* Dans la section "emprise", entrer les coordonnées des limites géographiques de chacun des bords du raster (Nord, Sud, Est et Ouest). S'assurer de choisir des valeurs qui englobent la partie visible de tous les rasters et noter ces valeurs dans un coin.
* Valider en cliquant sur "OK".
* Répéter cette opération pour les autres rasters à traiter, en utilisant les mêmes valeurs pour l'emprise.
* Importer ces calques dans QGis pour remplacer les raster n'ayant pas une origine commune.


Et voici la nouvelle version du script. Ce script n'est pas très flexible (il a été écrit pour résoudre un problème précis, si vous souhaitez le réutiliser, des ajustements seront peut-être nécessaires). Il travaille sur tous les calques présents dans la liste des couches, qu'ils soient masqués ou non, et dans l'ordre inverse où ils apparaissent dans cette liste (Le premier raster utilisé celui en bas de la liste, puis on écrira par dessus les valeurs du deuxième en partant du bas et ainsi de suite).

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, layerData <= 10000), layerData, dest)


# Write output to a file
newRasterfn = '/Users/sbooob/Desktop/QGis/SIG/merged.img'
array2raster(firstRasterPath, newRasterfn, dest, gdal.GDT_UInt16)

En tout cas merci à tous pour votre aide et vos conseil.

Bonne fin de journée

Tristan et Yo