Membaca, memodifikasi, dan menulis geotiff dengan GDAL dengan python

Saya mencoba mempelajari tali pemrosesan gambar Penginderaan Jauh menggunakan binding Python GDAL dan numpy. Sebagai upaya pertama, saya membaca file geotiff Landsat8, melakukan manipulasi sederhana dan menulis hasilnya ke file baru. Kode di bawah ini tampaknya berfungsi dengan baik, kecuali bahwa raster asli dibuang di file output, daripada raster yang dimanipulasi.

Setiap komentar atau saran dipersilahkan, tetapi terutama catatan tentang mengapa raster yang dimanipulasi tidak muncul dalam hasilnya.

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

Saya menggunakan Python 2.7.1 pada mesin Windows 7 32 bit.

Skrip Anda tidak memiliki metode ds.FlushCache, yang menyimpan ke disk apa yang Anda miliki di memori di akhir modifikasi. Lihat di bawah versi contoh Anda yang sudah diperbaiki. Perhatikan bahwa saya juga menambahkan dua baris untuk mengatur proyeksi dan geotransform sebagai input

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