11

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.

  • 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 '15 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. – HDR Oct 1 '15 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. – HDR Oct 1 '15 at 20:53
  • what would arr_min be? – fluidmotion Oct 2 '15 at 2:08
  • arr_min = 0. I've updated the post to show this. Thanks! – HDR Oct 2 '15 at 18:11
13

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()
[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
  • 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 '17 at 12:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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