GDAL WriteArray issue

Question

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()

I originally posted an answer then I realized I misunderstood what you were doing. So questions: have you tried moving the `print` statements to after `outDs.SetProjection(inDs.GetProjection())` and `outDs.SetGeoTransform(inDs.GetGeoTransform())` to make sure that's not the problem? If you end the script after `outDs = None`, then look at the file, is it empty? Confirm the problem is the write, not how you're reading it back. — agf, Jul 22 '11 at 15:25
I figured out the issue: in order to work with the newly created file and store the info, I needed to assign the driver.Create() to a variable: i.e. "outFile = driver.Create("C:\\Documents and Settings\\patrick\\Desktop\\tiff elevation\\EBK1KM\\newfile.tif",cols,rows,3,GDT_UInt16)" after that, I worked with the outFile's variables and data members to get the desired results. Thanks for your help! — Pat, Aug 08 '11 at 16:41

Mike T · Answer 1 · 2020-02-27T20:50:42.133

23

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

edited Feb 27 '20 at 20:50

answered Jul 28 '11 at 20:15

Mike T

41,085
18
152
203

it's the first line you gave me - setting the Create() equal to a variable solved everything, and i just needed to work with that variable after. thanks! – Pat Aug 08 '11 at 16:42
I guess @Pat forgot to do those steps. :-/ – Richard Aug 19 '14 at 15:52

GDAL WriteArray issue

1 Answers1

Linked