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


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)
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.

  • I got it to work on a DEM (Ubuntu, python 2.7.1) and it produced the expected result, with everything below the mean value set to 10000 and written to a new tiff. You're not copying the geotransform over to the new image so it's unprojected, so you might need to factor that in when trying to view it (there's a one-liner to do this but I'll need to dig it out). If you can edit your question with the output from gdainfo -stats original.tiff and gdal-config --version too that could help.
    – Steven Kay
    Sep 30, 2015 at 21:48
  • Hi, thanks for looking into this! I know I neglected the geotransform, thought to chew on that later. I do see the entire output image though (using Irfanview), so that can't be it I think. I'll generate the info you requested when I'm back on seat tonight. Oct 1, 2015 at 7:15
  • Hi, I'm struggling to provide the info you asked. I'm using Python GDAL binding and I'm unsure how the commands you specifiy correspond to a Python command. In any case, I'm using GDAL-1.11.2-cp27-none-win32, as acquired from here. I'll update my post with some stats on the original .tiff. Oct 1, 2015 at 20:53
  • what would arr_min be? Oct 2, 2015 at 2:08
  • arr_min = 0. I've updated the post to show this. Thanks! Oct 2, 2015 at 18:11

1 Answer 1


Your script is missing the ds.FlushCache method, that saves to disk what you have in memory at the end of the modifications. See below a corrected version of your example. Notice that I also added two lines to set projection and geotransform as input

import os
import gdal

file = "path+filename"
ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
[rows, cols] = 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, cols, rows, 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).SetNoDataValue(10000)##if you want these values transparent
outdata.FlushCache() ##saves to disk!!
outdata = None
  • The outfile is not projected. I am reading an HDF5 file and select the projection from the band which I want to export GetProjection() delivers the correct EPSG, but it seems to be not applied. GDAL warp missing? Thanks!
    – Michael
    Oct 28, 2017 at 12:41
  • what should I replace with outdata.GetRasterBand(1).WriteArray(arr_out) to write a multispectral image that has more than one band?
    – mArk
    Apr 5, 2020 at 19:07
  • The "1" in driver.Create() indicates the number of bands. Then you can write on each band with outdata.GetRasterBand(band_number). It starts at 1, not zero. Apr 6, 2020 at 9:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.