I'm trying the learn the ropes of Remote Sensing image processing using Python GDAL bindings and numpy. As a first attempt, I'm reading a Landsat8 geotiff file, do a simple manipulation and write the result to a new file. The code below appears to work fine, except that the original raster is dumped in the output file, rather than the manipulated raster.
Any comments or suggestions are welcome, but particularly notes on why the manipulated raster does not show in the result.
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
I use Python 2.7.1 on a Windows 7 32 bit machine.