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.
gdainfo -stats original.tiff
ygdal-config --version
también eso podría ayudar.Respuestas:
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
fuente
GetProjection()
entregar el EPSG correcto, pero parece que no se aplica. ¿Falta la deformación GDAL? ¡Gracias!outdata.GetRasterBand(1).WriteArray(arr_out)
para escribir una imagen multiespectral que tenga más de una banda?