Procesamiento de imágenes con Python, GDAL y Scikit-Image

Question

Estoy luchando con un procesamiento y espero poder resolverlo aquí.

Trabajo con la teledetección aplicada a la silvicultura, especialmente trabajando con datos LiDAR. La idea es utilizar Scikit-image para la detección de la copa de los árboles. Como soy nuevo en Python, consideré un gran triunfo personal hacer lo siguiente:

1. Importar un CHM (con matplotlib);
2. Ejecutar un filtro gaussiano (con el paquete scikit-image);
3. Ejecutar un filtro maxima (con el paquete scikit-image);
4. Ejecute el peak_local_max (con el paquete scikit-image);
5. Mostrar el CHM con los máximos locales (con matplotlib);

Ahora mi problema. Cuando importo con matplot, la imagen pierde sus coordenadas geográficas. Así que las coordenadas que tengo son sólo coordenadas básicas de la imagen (es decir, 250,312). Lo que necesito es obtener el valor del pixel bajo el punto máximo local de la imagen (puntos rojos en la imagen). Aquí en el foro he visto a un chico preguntando lo mismo ( ¿Obtener el valor del píxel de la trama GDAL bajo el punto OGR sin NumPy? ), pero ya tenía los puntos en un shapefile. En mi caso los puntos fueron calculados con scikit-image (Es un array con las coordenadas de cada copa de árbol). Así que no tengo el shapefile.

En conclusión, lo que quiero al final es un archivo txt con las coordenadas de cada máximo local en coordenadas geográficas, por ejemplo:

525412 62980123 1150 ...

Answer 1

En primer lugar, ¡bienvenido al sitio!

Los arrays Numpy no tienen un concepto de sistemas de coordenadas incorporado en el array. En el caso de una trama 2D, se indexan por columna y fila.

Nota Estoy haciendo la suposición de que usted está leyendo un formato de trama que es con el apoyo de GDAL .

En Python, la mejor manera de importar datos espaciales rasterizados es con la función rasterio paquete. Los datos brutos importados por rasterio siguen siendo un array de numpy sin acceso a los sistemas de coordenadas, pero rasterio también te da acceso a un método affine en el array fuente que puedes utilizar para transformar las columnas y filas del raster a coordenadas proyectadas. Por ejemplo:

import rasterio

# The best way to open a raster with rasterio is through the context manager
# so that it closes automatically

with rasterio.open(path_to_raster) as source:

    data = source.read(1) # Read raster band 1 as a numpy array
    affine = source.affine

# ... do some work with scikit-image and get an array of local maxima locations
# e.g.
# maxima = numpy.array([[0, 0], [1, 1], [2, 2]])
# Also note that convention in a numy array for a 2d array is rows (y), columns (x)

for point in maxima: #Loop over each pair of coordinates
    column = point[1]
    row = point[0]
    x, y = affine * (column, row)
    print x, y

# Or you can do it all at once:

columns = maxima[:, 1]
rows = maxima[:, 0]

xs, ys = affine * (columns, rows)

Y a partir de ahí puedes escribir tus resultados en un archivo de texto como quieras (te sugiero que eches un vistazo al incorporado csv módulo por ejemplo).

Answer 2

De un vistazo rápido a matplotlib, diría que tienes que alterar las escalas de los ejes después de la importación.

Answer 3

Por favor, inténtelo con el siguiente fragmento de código. Esto se puede utilizar para leer los datos de la imagen de la trama y escribir los datos procesados a la trama (archivo geotiff).

from PIL import Image,ImageOps
import numpy as np
from osgeo import gdal
#from osgeo import gdal_array
#from osgeo import osr
#from osgeo.gdalconst import *
#import matplotlib.pylab as plt

#from PIL import Image, ImageOps
#import gdal
#from PIL import Image
gdal.AllRegister()

################## Read Raster #################
inRaster='C:\python\Results\Database\Risat1CRS\CRS_LEVEL2_GEOTIFF\scene_HH\imagery_HH.tif'

inDS=gdal.Open(inRaster,1)
geoTransform = inDS.GetGeoTransform()
band=inDS.GetRasterBand(1)
datatype=band.DataType
proj = inDS.GetProjection()
rows = inDS.RasterYSize
cols=inDS.RasterXSize
data=band.ReadAsArray(0,0,cols,rows)#extraction of data to be processed#
############write raster##########
driver=inDS.GetDriver()
outRaster='C:\\python\\Results\\Database\\Temporary data base\\clipped_26July2017\\demo11.tif'
outDS = driver.Create(outRaster, cols,rows, 1,datatype)
geoTransform = inDS.GetGeoTransform()
outDS.SetGeoTransform(geoTransform)
proj = inDS.GetProjection()
outDS.SetProjection(proj)
outBand = outDS.GetRasterBand(1)
outBand.WriteArray(data1,0,0)
#data is the output array to written in tiff file
outDS=None 
im2=Image.open('C:\\python\\Results\\Database\\Temporary data base\\clipped_26July2017\\demo11.tif');
im2.show()