2

I'm trying to create a new raster from three bands (rgb). I got the script mostly from here. I can create a single band raster with no problems, but I'm having trouble inputting a tree band numpy array. I'm greeted with this error message when I attempt it.

ValueError: expected array of dim 2

If I try to create a three band GeoTIFF image I end up with a grey,alpha,alpha image instead of an rgb image, nor am I able to populate these with my three-dimensional numpy array. The following script works for creating a greyscale raster. with three bands as grey,alpha,alpha.

from osgeo import gdal, osr
import os
import io
import numpy as np

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

    bands = array.shape[0]
    cols = array.shape[1]
    rows = array.shape[2]

    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName(format)
    outRaster = driver.Create(newRasterfn, cols, rows, bands,gdal.GDT_UInt16,options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF',])
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

    for band in range(bands):
        outRaster.GetRasterBand(band+1).WriteArray(array[band,:,:])
        outRaster.GetRasterBand(band+1).FlushCache()
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(spatialRef)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())


def get_rgb_bands():
    source = gdal.Open("C:\\RGBtest.tif")

    red = source.GetRasterBand(1).ReadAsArray()
    green = source.GetRasterBand(2).ReadAsArray()
    blue = source.GetRasterBand(3).ReadAsArray()

    #rgbOutput = source.ReadAsArray() #Easier method
    rgbOutput = np.zeros((3,16,16), 'uint16')

    rgbOutput[0,...] = red
    rgbOutput[1,...] = green
    rgbOutput[2,...] = blue

    #Clear so file isn't locked
    source = None
    return rgbOutput

exampleRGB = get_rgb_bands()

exampleOutput = "C:\\ExampleOutput.tif"
rasterOrigin=[0,0]

array2raster(exampleOutput,rasterOrigin,1, 1,exampleRGB,'GTiff',3879) 

createdFile = gdal.Open("C:\\ExampleOutput.tif")

print 'old file info:'+gdal.Info(source)
print 'new file info:'+gdal.Info(createdFile)

createdFile = None

exampleOutput = None

EDIT: Luke answered this question perfectly I have updated the code to reflect a working example to output a 16-bit rgb raster.

2

You can set the colour interpretation for most rasters using Band.SetColorInterpretation, i.e

outband.SetColorInterpretation(gdal.GCI_RedBand)

However, this is only understood by GDAL based software. To specifically create a GeoTIFF with an RGB baseline PhotometricInterpretation tag (this is not a GDAL tag, but part of the TIFF spec) pass a PHOTOMETRIC=RGB option to the GeoTIFF driver when creating the raster:

driver = gdal.GetDriverByName('GeoTIFF')
options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']
outRaster = driver.Create(filepath, cols, rows, bands, datatype, options=options)

You're also trying to write 3 bands to a single band:

(outRaster.GetRasterBand(1).WriteArray(etc...))

Get used to having your arrays in [bands, rows, cols] order not [cols, rows, bands] as this is how gdal arranges them when reading into a numpy array. Then you need to loop through the array "bands" and write each of them to the raster bands:

bands = array.shape[0]
for band in range(bands):
    outRaster.GetRasterBand(band+1).WriteArray( array[band, :, :] )

Try something like:

from osgeo import gdal, osr
import os
import io
import numpy as np

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,format='GTiff'):

    bands = array.shape[0]
    rows = array.shape[1]
    cols = array.shape[2]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    driver = gdal.GetDriverByName(format)
    #Here is where I assign three bands to the raster with int 3
    options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']
    outRaster = driver.Create(newRasterfn, cols, rows, bands, gdal.GDT_UInt16, options=options)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    #outband = outRaster.GetRasterBand(1)
    #outband.WriteArray(array)
    for band in range(bands):
        outRaster.GetRasterBand(band+1).WriteArray( array[band, :, :] )

    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(6962)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())

    return outRaster


def example_rgb_creation():
    source = gdal.Open("C:\\RGBtest.tif")
    rgb_array = source.ReadAsArray()

    return rgb_array

exampleRGB = example_rgb_creation()
exampleOutput = "C:\\Temp\\ExampleOutput.tif"
rasterOrigin=[0,0]

outRaster = array2raster(exampleOutput,rasterOrigin,1, 1,exampleRGB)

outRaster.GetRasterBand(1).SetColorInterpretation(gdal.GCI_RedBand)
outRaster.GetRasterBand(2).SetColorInterpretation(gdal.GCI_GreenBand)
outRaster.GetRasterBand(3).SetColorInterpretation(gdal.GCI_BlueBand)

del outRaster
  • Thank you very much. I got to create a tiff file but it is unreadable now. I populated separate rgb bands and then populated them from separate rgb arrays. But I think that I'm not writing the actual driver correctly. I don't understand how the driver gets the bands. – Školstvo Jul 18 '17 at 13:14
  • I tried different ways but the bands still stay as grey, alpha, alpha instead of r,g,b – Školstvo Jul 18 '17 at 18:46
  • Your edit made everything clear, except that writing to a 16-bit GeoTIFF still creates bands as grey, alpha, alpha instead of r,g,b. If I write gdal.GDT_Byte thus writing an 8-bit raster, everything works beautifully (although clamped to 256 colors). The problem is that I absolutely need the extra dynamic range, since the raster will be used by a different program in the future. – Školstvo Jul 19 '17 at 6:59
  • The raster I'm outputting is a normalmap for another program. The channels matter very much, and I need to use georeferenced coordinates. I found this documentation about designating specific color-types but I'm unsure how to assign this data to the driver? – Školstvo Jul 19 '17 at 7:21
  • The other program I am using is already compiled to work with RGB color channels. Currently Gdal outputs three channels as greyscale, alpha, alpha. The program interprets this as greyscale = red, greyscale = green, greyscale = blue, alpha1 = transparency, alpha2 = nothing. I can see the channels assigned in Photoshop. – Školstvo Jul 19 '17 at 8:10

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.