2

I have 9 vector layers that I rasterized. All of these have the same:

  • x;y origin
  • pixel size (1 pixel=1 meter)
  • SRS

When displayed separately in QGIS these 9 rasters are superimposed correctly but when merge using the Python script there is one or more rasters that are offset approximately 100 meters to the North.

Can you check if the script is correct or if the problem come from my data?

import os
import gdal
import struct
import numpy
import osr
import sys

def raster2array(raster):
    band = raster.GetRasterBand(1)
    return band.ReadAsArray()

def getOrigin(raster):
    geotransform = raster.GetGeoTransform()
#    print geotransform
    return [geotransform[0], geotransform[3]]

def getLayers():
    return iface.legendInterface().layers()

def array2raster(rasterfn,newRasterfn,array, xPadding = 0, yPadding = 0):
    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('GTiff')
    outRaster = driver.Create(newRasterfn, cols + xPadding, rows + yPadding, 1, gdal.GDT_Float64)
    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()

layers = getLayers()
firstLayer = layers.pop()
path = firstLayer.dataProvider().dataSourceUri()
raster = gdal.Open(path)

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)

    for x in range(0, len(layerData)):
        for y in range(0, len(layerData[x])):
            if layerData[x][y] > 0:
                dest[x][y] = layerData[x][y]

rasterfn = '/Users/sbooob/Desktop/QGis/SIG/1CLC12_Raster.img'
newValue = 0
newRasterfn = '/Users/sbooob/Desktop/QGis/SIG/out.tif'

array2raster(rasterfn, newRasterfn, dest)    

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.