Optimización de la identificación de los límites de las masas arbóreas

Question

Estoy trabajando en un flujo de trabajo para identificar los límites de las masas arbóreas. Mi resultado deseado son las líneas azules - sé que nunca será así, pero necesito que sea lo mejor posible. Mi mejor resultado son las líneas rojas. El resultado final no deben ser polígonos, sino líneas.

Esto es lo que he hecho hasta ahora: Estoy usando un DSM normalizado con tamaño de píxel 2x2 como mi capa de partida. Este es mi flujo de trabajo:

Utilizo i.segment donde mis mejores resultados han sido con Umbral de diferencia=0,4 y Número mínimo de celdas = 150. Entonces utilizo r.to.vect porque Polygonize por alguna razón no funcionaba en este flujo de trabajo (?). Después de esto convierto a líneas, y borro los duplicados. Uso la simplificación de líneas con Tolerence 8.5 (mejor resultado) Y por último uso v.generalize con valor de tolerancia máxima=21 con el algoritmo Douglas.

Así que utilizo cuatro variables para mejorar mi resultado. He intentado utilizar algunas más de las otras variables, pero sinceramente no estoy seguro de cómo afectan al resultado.

Así que pregunto aquí por cualquier sugerencia para mejorar mi resultado.

Actualización: Me las arreglé para hacerlo un poco mejor por seguir cambiando las variables anteriores. Probé algunas otras variables:

• En segmentación cambié la cantidad de memoria a usar en MB, pero no tuvo efecto.
• En v.generalize He probado todos los algoritmos diferentes, y Douglas parece ser el mejor.
• También he probado a cambiar el parámetro Look-ahead en v.generalize pero no tuvo ningún efecto positivo.

Así que sigo esperando alguna sugerencia "fuera de lo común" para mejorarlo.

Actualización-2: Todavía estoy trabajando en esto. He probado diferentes variaciones para las entradas raster. He intentado lo siguiente:

1. DSM normalizado (resolución 2x2)
2. MDS normalizado (2x2) + ortofoto
3. DSM normalizado (2x2) + ortofoto + ortofoto CIR
4. Ortofoto + Ortofoto CIR
5. DSM normalizado (2x2) + DSM (2x2)

También he probado diferentes resoluciones, pero 2x2 da el mejor resultado. También he aprendido que la combinación de tres entradas de trama da un resultado peor. Mis mejores resultados son o bien utilizando sólo Normalized DSM (2x2) o la entrada combinada con Normalized DSM (2x2) + DSM (2x2).

Sigo esperando sugerencias. Estoy abierto a probar otro producto de código abierto.

Answer 1

He probado con Python. Espero que esto sigue siendo aceptable. Básicamente, hice un poco de desenfoque (suavizado gaussiano, puedes desenfocar la imagen más o menos para hacer las áreas homogéneas más homogéneas), y luego apliqué un algoritmo de segmentación. Mi script se basa en gran medida en algunos documentación sobre scikit-image y esto es lo que obtengo en términos de segmentos:

Por supuesto, puede ajustar los parámetros, y si desea guardar el resultado como un GeoTIFF, por ejemplo, para extraer polígonos, sólo tiene que guardar el archivo segments_slic como GeoTIFF con la geolocalización adecuada. He dejado esa parte del código comentada al final.

# Big import block
from PIL import Image
import matplotlib.pyplot as plt
import numpy as np

from skimage.color import rgb2gray
from skimage.segmentation import slic
from skimage.segmentation import mark_boundaries

# Open image
im_fp = Image.open("forest.png")
# Turn into numpy array, discard last band
# (alpha channel?)
img = np.asarray(im_fp)[:, :, :-1]

# Smooth spatially a bit
img=ndi.gaussian_filter(img, sigma=5)

# Calculate segments. 
# Feel free to tweak n_segments, compactness
# and sigma
segments_slic = slic(img, n_segments=25, compactness=50, 
                     sigma=1, start_label=1)

# Plot image and boundaries
plt.imshow(mark_boundaries(img, segments_slic))

# Export as GeoTIFF. Check
#drv = gdal.GetDriverByName('GTiff')     # create driver for writing geotiff file
#outRaster = drv.CreateCopy('forest_segments.tif', 'forest.tif, , 0 )   # create new copy of inut raster on disk
#outRaster.SetGeoTransform([top_left_x, hor_resolution, 0, top_left_y, 0, -vert_res])
#outRater.SeProjection()
#newBand = outRaster.GetRasterBand(1)                               # get the first (and only) band of the new copy
#newBand.WriteArray(segments_slic)                                           # write array data to this band 
#outRaster = None