Lesen, Ändern und Schreiben eines Geotiffs mit GDAL in Python

Ich versuche, die Seile der Fernerkundungs-Bildverarbeitung mit Python GDAL-Bindungen und Numpy zu lernen. Als ersten Versuch lese ich eine Landsat8-Geotiff-Datei, führe eine einfache Manipulation durch und schreibe das Ergebnis in eine neue Datei. Der folgende Code scheint einwandfrei zu funktionieren, außer dass das ursprüngliche Raster in der Ausgabedatei und nicht im manipulierten Raster gespeichert wird.

Kommentare oder Vorschläge sind willkommen, insbesondere aber Hinweise, warum das manipulierte Raster im Ergebnis nicht angezeigt wird.

import os
import gdal

gdal.AllRegister()

file = "c:\~\LC81980242015071LGN00.tiff"
(fileRoot, fileExt) = os.path.splitext(file)
outFileName = fileRoot + "_mod" + fileExt

ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()

[cols, rows] = arr.shape
arr_min = arr.Min()
arr_max = arr.Max()
arr_mean = int(arr.mean())

arr_out = numpy.where((arr < arr_mean), 10000, arr)

driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(outFileName, rows, cols, 1, gdal.GDT_UInt16)
outband = outdata.GetRasterBand(1)
outband.WriteArray(arr_out)
outdata = None

print arr_min
> 0
print arr_max
> 65535
print arr_mean
> 4856

Ich verwende Python 2.7.1 auf einem Windows 7 32-Bit-Computer.

In Ihrem Skript fehlt die Methode ds.FlushCache, mit der gespeichert wird, was Sie am Ende der Änderungen im Speicher haben. Unten sehen Sie eine korrigierte Version Ihres Beispiels. Beachten Sie, dass ich auch zwei Zeilen hinzugefügt habe, um Projektion und Geotransform als Eingabe festzulegen

import os
import gdal

file = "path+filename"
ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
[cols, rows] = arr.shape
arr_min = arr.min()
arr_max = arr.max()
arr_mean = int(arr.mean())
arr_out = numpy.where((arr < arr_mean), 10000, arr)
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(outFileName, rows, cols, 1, gdal.GDT_UInt16)
outdata.SetGeoTransform(ds.GetGeoTransform())##sets same geotransform as input
outdata.SetProjection(ds.GetProjection())##sets same projection as input
outdata.GetRasterBand(1).WriteArray(arr_out)
outdata.GetRasterBand(1).SetNoDataValue(10000)##if you want these values transparent
outdata.FlushCache() ##saves to disk!!
outdata = None
band=None
ds=None