I'm utilizing python GDAL to write a raster data into a .tif file. Here's the code:
import numpy, sys
from osgeo import gdal, utils
from osgeo.gdalconst import *
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
inDs = gdal.Open("C:\\Documents and Settings\\patrick\\Desktop\\tiff elevation\\EBK1KM\\color_a1.tif",GDT_UInt16)
if inDs is None:
print "couldn't open input dataset"
sys.exit(1)
else:
print "opening was successful!"
cols = inDs.RasterXSize
rows = inDs.RasterYSize
bands = inDs.RasterCount
driver = inDs.GetDriver()
driver.Create("C:\\Documents and Settings\\patrick\\Desktop\\tiff elevation\\EBK1KM\\newfile.tif",cols,rows,3,GDT_UInt16)
outDs = gdal.Open("C:\\Documents and Settings\\patrick\\Desktop\\tiff elevation\\EBK1KM\\newfile.tif")
if outDs is None:
print "failure to create new file"
sys.exit(1)
outBand1 = outDs.GetRasterBand(1)
outBand2 = outDs.GetRasterBand(2)
outBand3 = outDs.GetRasterBand(3)
data1 = inDs.GetRasterBand(1).ReadAsArray()
data2 = inDs.GetRasterBand(2).ReadAsArray()
data3 = inDs.GetRasterBand(3).ReadAsArray()
outBand1.WriteArray(data1,0,0)
outBand2.WriteArray(data2,0,0)
outBand3.WriteArray(data3,0,0)
print "before closing out the file"
print outDs.GetRasterBand(1).ReadAsArray(700,700,5,5)
print outDs.GetRasterBand(2).ReadAsArray(700,700,5,5)
print outDs.GetRasterBand(3).ReadAsArray(700,700,5,5)
outDs.SetProjection(inDs.GetProjection())
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs = None
outDs = gdal.Open("C:\\Documents and Settings\\patrick\\Desktop\\tiff elevation\\EBK1KM\\newfile.tif")
print "after reopening"
print outDs.GetRasterBand(1).ReadAsArray(700,700,5,5)
print outDs.GetRasterBand(2).ReadAsArray(700,700,5,5)
print outDs.GetRasterBand(3).ReadAsArray(700,700,5,5)
The resultant output between the closing and reopening of the output dataset are different:
before closing out the file
[[ 36 35 55 121 0]
[ 54 0 111 117 0]
[ 0 117 152 56 0]
[ 89 122 56 0 0]
[102 107 0 25 53]]
[[ 68 66 126 200 0]
[ 78 0 166 157 0]
[ 0 235 203 70 0]
[229 251 107 0 0]
[241 203 0 42 121]]
[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]
after reopening
[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]
[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]
[[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]
[0 0 0 0 0]]
is there some command I'm missing to ensure that the file is written and saved prior to setting the variable to None? I've tried adding both of the following with no luck:
outband1.FlushCache()
outDs.FlushCache()
Answer
You don't need to Create then Open a raster (which you were reading GA_ReadOnly). You also don't need gdal.AllRegister() at the beginning, as it has already been called when you load GDAL into Python (see the Raster API tutorial).
Picking up somewhere above (with modifications):
# Create a new raster data source
outDs = driver.Create(out_fname, cols, rows, 3, gdal.GDT_UInt16)
# Write metadata
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())
# Write raster data sets
for i in range(3):
outBand = outDs.GetRasterBand(i + 1)
outBand.WriteArray(data[i])
# Close raster file
outDs = None
Sometimes I add this to ensure the file is fully deallocated, and to prevent running into some gotchas:
del outDs, outBand