10

我正在利用 python GDAL 将栅格数据写入 .tif 文件。这是代码：

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)

输出数据集关闭和重新打开之间的结果输出不同：

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]]

在将变量设置为无之前，我是否缺少一些命令来确保写入和保存文件？我尝试添加以下两个都没有运气：

outband1.FlushCache()
outDs.FlushCache()
4

1 回答 1

22

你不需要Create然后Open一个光栅（你正在阅读GA_ReadOnly）。您一开始也不需要gdal.AllRegister()，因为在将 GDAL 加载到 Python 时已经调用了它（请参阅Raster API 教程）。

在上面某处拾取（有修改）：

# 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

有时我添加它以确保文件被完全释放，并防止遇到一些陷阱：

del outDs, outBand
于 2011-07-28T20:15:35.310 回答