Leer, modificar y escribir un geotiff con GDAL en python

11

Estoy tratando de aprender las cuerdas del procesamiento de imágenes de detección remota utilizando enlaces Python GDAL y numpy. Como primer intento, estoy leyendo un archivo geotiff Landsat8, hago una simple manipulación y escribo el resultado en un nuevo archivo. El siguiente código parece funcionar bien, excepto que el ráster original se volca en el archivo de salida, en lugar del ráster manipulado.

Cualquier comentario o sugerencia son bienvenidos, pero particularmente notas sobre por qué el ráster manipulado no se muestra en el resultado.

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

Yo uso Python 2.7.1 en una máquina con Windows 7 de 32 bits.

HDR
fuente
Lo hice funcionar en un DEM (Ubuntu, python 2.7.1) y produjo el resultado esperado, con todo por debajo del valor medio establecido en 10000 y escrito en un nuevo tiff. No está copiando la geotransformación a la nueva imagen para que no se proyecte, por lo que es posible que tenga que tener esto en cuenta al intentar verla (hay una línea para hacer esto, pero tendré que desenterrarla). Si puede editar su pregunta con la salida de gdainfo -stats original.tiffy gdal-config --versiontambién eso podría ayudar.
Steven Kay
Hola, gracias por investigar esto! Sé que descuidé la geotransformación, pensé en masticar eso más tarde. Sin embargo, veo toda la imagen de salida (usando Irfanview), así que creo que no puede ser. Generaré la información que solicitaste cuando vuelva al asiento esta noche.
HDR
Hola, me cuesta proporcionar la información que solicitaste. Estoy usando el enlace Python GDAL y no estoy seguro de cómo los comandos que especifique corresponden a un comando de Python. En cualquier caso, estoy usando GDAL-1.11.2-cp27-none-win32, tal como lo adquirí desde aquí . Actualizaré mi publicación con algunas estadísticas sobre el .tiff original.
HDR
¿Qué sería arr_min?
fluidmotion
arr_min = 0. He actualizado la publicación para mostrar esto. ¡Gracias!
HDR

Respuestas:

13

A su script le falta el método ds.FlushCache, que guarda en el disco lo que tiene en la memoria al final de las modificaciones. Vea a continuación una versión corregida de su ejemplo. Tenga en cuenta que también agregué dos líneas para establecer la proyección y la geotransformación como entrada

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
Andrea Massetti
fuente
El archivo no se proyecta. Estoy leyendo un archivo HDF5 y selecciono la proyección de la banda que quiero exportar para GetProjection()entregar el EPSG correcto, pero parece que no se aplica. ¿Falta la deformación GDAL? ¡Gracias!
Michael
¿con qué debo reemplazar outdata.GetRasterBand(1).WriteArray(arr_out)para escribir una imagen multiespectral que tenga más de una banda?
Sai Kiran
El "1" en driver.Create () indica el número de bandas. Luego puede escribir en cada banda con outdata.GetRasterBand (band_number). Comienza en 1, no en cero.
Andrea Massetti hace