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
With some help from Jon, here is an answer:
set pathin=D:\procbox\In
set pathout=D:\procbox\Out
FOR %%N in (%pathin%\*.tif) DO (^
gdal_calc.py -A %%N --A_band=1 --outfile=A.tif --calc="A" --overwrite --type="UInt16" & ^
gdal_calc.py -A %%N --A_band=2 --outfile=B.tif --calc="A" --overwrite --type="UInt16" & ^
gdal_calc.py -A %%N --A_band=3 --outfile=C.tif --calc="A" --overwrite --type="UInt16" & ^
gdal_calc.py -A A.tif --A_band=1 -B B.tif --B_band=1 -C C.tif --C_band=1 ^
--outfile=%pathout%\%%~nN-cov.tif --calc="255*((A+B+C)==765)" --overwrite --type="UInt16")
FOR %%N in (%pathout%\*-cov.tif) DO gdal_polygonize.py %%N -f "ESRI Shapefile" -mask %%N %pathout%\%%~nN.shp
What this script does is it sequentially processes all .tif files in a directory, breaks them down into individual RGB bands (A,B,C .tif) and then calculates a new raster with 255 whenever all these three rasters add up to white (765). What was important was that I was working with a UInt8 tiff, and the input tiffs need to be at least UInt16 for the calculation not to loop over 255.
Then it creates a shapefile around all the white areas. I use this to check survey coverage. If you use the script, please note that in Windows BAT files you break rows with the ^ symbol, rather than with a backslash escape (\) like with shell scripting in Linux.
Best Answer
here's a better answer, use gdalbuildvrt with either srcnodata or vrtnodata flag:
If the next application in line doesn't understand .vrt, translate to a new tif: