Revisions to How to create an rgb GeoTIFF file raster from bands using the GDAL python module?

deleted 354 characters in body

Source Link

edited Jul 19, 2017 at 11:29

233
1
2
9

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):
        if band is 0:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
        if band is 1:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_GreenBand)
        if band is 2:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
        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

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):
        if band is 0:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
        if band is 1:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_GreenBand)
        if band is 2:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
        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

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

deleted 154 characters in body

Source Link

edited Jul 19, 2017 at 11:19

Školstvo

233
1
2
9

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):
        if band is 0:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
        if band is 1:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_GreenBand)
        if band is 2:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
        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 editedhave updated the scriptcode to reflect Lukes recomendations. The file is written but I think there is a bug in GDAL where it doesn't actually writeworking example to output a 16-bit rgb image to the RGD image channels. The output still remains as grey,alpha,alpha.

Gdal.info says that the image is formatted correctly but that is of little consolementraster. [1]: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

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)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    
    for band in range(bands):
        if band is 0:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
        if band is 1:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_GreenBand)
        if band is 2:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
        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: I edited the script to reflect Lukes recomendations. The file is written but I think there is a bug in GDAL where it doesn't actually write a 16-bit rgb image to the RGD image channels. The output still remains as grey,alpha,alpha.

Gdal.info says that the image is formatted correctly but that is of little consolement. [1]: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

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):
        if band is 0:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
        if band is 1:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_GreenBand)
        if band is 2:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
        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. [1]: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

Further clarification

Source Link

edited Jul 19, 2017 at 11:06

Školstvo

233
1
2
9

I'm trying to create a new raster from three bands (rgb). I got the script mostly from here.[here.][1] 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.

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

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,format,spatialRef):
    cols
    bands = array.shape[1]shape[0]
    rowscols = array.shape[0]shape[1]
    originXrows = rasterOrigin[0]array.shape[2]
    originY = rasterOrigin[1]
    driveroriginX = gdal.GetDriverByName('GTiff')rasterOrigin[0]
    #HereoriginY is= whererasterOrigin[1]
 I assign three bands 
 to the raster withdriver int= 3gdal.GetDriverByName(format)
    outRaster = driver.Create(newRasterfn, cols, rows, 3bands, gdal.GDT_UInt16)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband 
 = outRaster.GetRasterBand  for band in range(1bands):
    outband.WriteArray(array)    if band is 0:
    outRasterSRS = osr      outRaster.SpatialReferenceGetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
    outRasterSRS.ImportFromEPSG(6962)    if band is 1:
            outRaster.SetProjectionGetRasterBand(outRasterSRSband+1).ExportToWktSetColorInterpretation()gdal.GCI_GreenBand)
    outband    if band is 2:
            outRaster.FlushCacheGetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
    
#The source image is loadedoutRaster.GetRasterBand(band+1).WriteArray(array[band,:,:])
 from disk in this example for simplicity's sakeoutRaster.GetRasterBand(band+1).FlushCache()
 It is 16x16 pixelsoutRasterSRS big= osr.SpatialReference()
 It is also aoutRasterSRS.ImportFromEPSG(spatialRef)
 16-bit image, but thatoutRaster.SetProjection(outRasterSRS.ExportToWkt())
 is irrelevant  

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

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

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

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

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

exampleRGB = example_rgb_creationget_rgb_bands()

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

#We call thearray2raster(exampleOutput,rasterOrigin,1, actual1,exampleRGB,'GTiff',3879) rasterization

createdFile function= heregdal. Note that I onlyOpen("C:\\ExampleOutput.tif")

print send'old onwfile bandinfo:'+gdal. SendingInfo(source)
print three'new returnsfile theinfo:'+gdal.Info(createdFile)

createdFile above= error.None
array2raster( 
exampleOutput,rasterOrigin,1, 1,exampleRGB[:,:,1])= None

EDIT: I edited the script to reflect Lukes recomendations. The file is written but I think there is a bug in GDAL where it doesn't actually write a 16-bit rgb image to the RGD image channels. The output still remains as grey,alpha,alpha.

Gdal.info says that the image is formatted correctly but that is of little consolement. [1]: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

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.

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

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,format):
    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    driver = gdal.GetDriverByName('GTiff')
    #Here is where I assign three bands to the raster with int 3
    outRaster = driver.Create(newRasterfn, cols, rows, 3, gdal.GDT_UInt16)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(6962)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
    
#The source image is loaded from disk in this example for simplicity's sake. It is 16x16 pixels big. It is also a 16-bit image, but that is irrelevant
def example_rgb_creation():
    source = gdal.Open("C:\\RGBtest.tif")

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

    rgb_array = np.zeros((16,16,3), 'uint16')

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

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

exampleRGB = example_rgb_creation()

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

#We call the actual rasterization function here. Note that I only send onw band. Sending three returns the above error.
array2raster(exampleOutput,rasterOrigin,1, 1,exampleRGB[:,:,1])

I'm trying to create a new raster from three bands (rgb). I got the script mostly from [here.][1] 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.

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)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
     
    for band in range(bands):
        if band is 0:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_RedBand)
        if band is 1:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_GreenBand)
        if band is 2:
            outRaster.GetRasterBand(band+1).SetColorInterpretation(gdal.GCI_BlueBand)
        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: I edited the script to reflect Lukes recomendations. The file is written but I think there is a bug in GDAL where it doesn't actually write a 16-bit rgb image to the RGD image channels. The output still remains as grey,alpha,alpha.

Gdal.info says that the image is formatted correctly but that is of little consolement. [1]: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

Source Link

asked Jul 17, 2017 at 15:23

Školstvo

233
1
2
9

Loading

Stack Exchange Network

Return to Question