Following the link in the comment from @kersten I found this https://svn.osgeo.org/gdal/trunk/autotest/gcore/tiff_write.py with lots of interesting snippets. I think what you want to do is this:
...
ct = gdal.ColorTable()
# Some examples
ct.SetColorEntry( 0, (0, 0, 0, 255) )
ct.SetColorEntry( 1, (0, 255, 0, 255) )
ct.SetColorEntry( 2, (255, 0, 0, 255) )
ct.SetColorEntry( 3, (255, 0, 255, 255) )
# Set the color table for your band
outband.SetColorTable(color_table)
...
Also, you may want to use CreateColorRamp
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
Best Answer
You are not defining a spatial reference for the shapefile in your code. In order to do so, you have to specify it as the second argument when using the
CreateLayer()method. For example, if your spatial reference would happen to be WGS84, you could write the following:I guess you want the spatial reference to be the same as the raster, so you could do the following as well:
As for the intervals, I am not familiar with the
ContourGenerate()function. However, you can easily create a list with intervals using Python'srange()function by specifying a step as the third argument. For example:Then, you can simply specify the intervals you just created as the argument.