Lire, modifier et écrire un géotiff avec GDAL en python

J'essaie d'apprendre les ficelles du traitement d'image de télédétection en utilisant les liaisons Python GDAL et numpy. Dans un premier temps, je lis un fichier géotiff Landsat8, je fais une manipulation simple et j'écris le résultat dans un nouveau fichier. Le code ci-dessous semble fonctionner correctement, sauf que le raster d'origine est sauvegardé dans le fichier de sortie, plutôt que le raster manipulé.

Tous les commentaires ou suggestions sont les bienvenus, mais en particulier les notes expliquant pourquoi le raster manipulé n'apparaît pas dans le résultat.

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

J'utilise Python 2.7.1 sur une machine Windows 7 32 bits.

Votre script ne contient pas la méthode ds.FlushCache, qui enregistre sur disque ce que vous avez en mémoire à la fin des modifications. Voir ci-dessous une version corrigée de votre exemple. Notez que j'ai également ajouté deux lignes pour définir la projection et la géotransformation en entrée

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