2

I have a numpy.array which I would like to save as GeoTiff using GDAL in python. The array was created using the following PCS in mind:

+proj=lcc +lat_1=48 +lat_2=48 +lat_0=48 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

where lat_0=48 and lon_0=0 are the standard latitude and longitude where there should be no distortion. I looked at this guide but with no success. The problem is that when I define the rasterOrigin in degrees and define pixelWidth and pixelHeight in metres this does not seem to work and when I load the raster in qgis the raster is misplaced. See the example below taken from the guide I pointed out but tweaked. The

outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

does not seem to work and the image is placed at the bottom left corner taken from the lat/lon details in the PCS and not from the SetGeoTransform which are different.

import gdal, ogr, os, osr
import numpy as np


def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):

    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.SetProjCS( "myproj" )
    outRasterSRS.SetWellKnownGeogCS( "WGS84")
    outRasterSRS.SetLCC(53.5,53.5,53.5,0,0,0)
    outRasterSRS.SetLinearUnits("Meter", pixelWidth)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    reversed_arr = array[::-1] # reverse array so the tif looks like the array
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster

if __name__ == "__main__":
rasterOrigin = (-20,40)
pixelWidth = 10
pixelHeight = 10
newRasterfn = 'test.tif'
array = np.array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])


main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)
  • pixel height will almost certainly be -10 (to do with image vs geographic origins) – Ian Turton Nov 2 '15 at 15:35
  • @iant pixel height is in metres - ie 1 pixel every 10 metres – Scott Nov 2 '15 at 15:57
  • But the image origin is at the top and the geographic origin is at the bottom, thus commonly we define rasters as having negative height to make the maths easier. So you pixel height will be -10m (unless the GDAL python bindings take care of that for you) – Ian Turton Nov 2 '15 at 15:59
  • @iant that is done by reversed_arr up in the main func but the problem I am having is to do with where the image is placed on the canvas rather that that. pixel height and width are fine. – Scott Nov 2 '15 at 16:01

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.