Look at the documentation for driver.Create. Your current code is trying to write a 3D array to a single band. That is what raise ValueError("expected array of dim 2")
is telling you. You need to create a tiff with bands equal to the number of "stacked images". Then write a single image to each band. You cannot write all the bands at once.
Here is the same function you wrote, updated to reflect this change. There are also modifications because NewFileName isn't the parameter name....
def CreateGeoTiff(Name, Array, driver, NDV,
xsize, ysize, GeoT, Projection, DataType):
Array[np.isnan(Array)] = NDV
DataSet = driver.Create(Name, xsize, ysize, Array.shape[0], DataType)
DataSet.SetGeoTransform(GeoT)
DataSet.SetProjection( Projection.ExportToWkt() )
for i, image in enumerate(Array, 1):
DataSet.GetRasterBand(i).WriteArray( image )
DataSet.GetRasterBand(i).SetNoDataValue(NDV)
DataSet.FlushCache()
return Name
You can make this less verbose because your array will tell you the xsize and ysize.
def CreateGeoTiff(Name, Array, driver, NDV,
GeoT, Projection, DataType):
Array[np.isnan(Array)] = NDV
DataSet = driver.Create(Name, Array.shape[2], Array.shape[1], Array.shape[0], DataType)
DataSet.SetGeoTransform(GeoT)
DataSet.SetProjection( Projection.ExportToWkt() )
for i, image in enumerate(Array, 1):
DataSet.GetRasterBand(i).WriteArray( image )
DataSet.GetRasterBand(i).SetNoDataValue(NDV)
DataSet.FlushCache()
return Name
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
Why are you using for saving Scipy? I think you can use function imwrite from CV2.
For example:
cv2.imwrite('D:/trial/blur_55.tif', myIMG)
or you can use also another libraries. For exammple PIL or Matplotlib.
Example for Matplotlib: