Extracción de áreas de copa de árbol de datos de detección remota (imágenes visuales y LiDAR)

13

Estoy buscando un método para procesar una imagen de detección remota y extraer las áreas de copa de los árboles individuales de la imagen.

Tengo imágenes de área visual de longitud de onda visual y datos LIDAR del área. La ubicación en cuestión es un área desértica, por lo que la cubierta arbórea no es tan densa como un área forestal. La resolución de las imágenes aéreas es de 0.5 pies por 0.5 pies. La resolución lidar es de aproximadamente 1 x 1 pie. Tanto los datos visuales como el lidar provienen de un conjunto de datos del condado de Pima, Arizona. Una muestra del tipo de imágenes aéreas que tengo está al final de esta publicación.

Esta pregunta ¿ Detección de árbol único en ArcMap? parece ser el mismo problema, pero no parece haber una buena respuesta allí.

Puedo obtener una clasificación razonable de los tipos de vegetación (e información sobre el porcentaje de cobertura general) en el área utilizando la clasificación Iso Cluster en Arcmap, pero esto proporciona poca información sobre árboles individuales. Lo más cercano que tengo a lo que quiero son los resultados de pasar la salida de la clasificación de isocluster a través de la característica Ráster a polígono en Arcmap. El problema es que este método se fusiona cerca de los árboles en un solo polígono.

Editar: Probablemente debería haber incluido más detalles sobre lo que tengo. Los conjuntos de datos en bruto que tengo son:

  • Los datos completos y un ráster tiff generado a partir de ellos.
  • Imágenes visuales (como la imagen de muestra que se muestra, pero que cubre un área mucho más amplia)
  • Mediciones directas manuales de un subconjunto de los árboles en el área.

De estos he generado:

  1. Las clasificaciones de suelo / vegetación.
  2. Los rásteres DEM / DSM.

muestra de imágenes aéreas

Theodore Jones
fuente
Tienes más datos que el enlace. ¿Tiene los archivos clasificados o solo el ráster DEM / DSM (¿cuál?) Es realmente no es fácil de hacer esto con longitudes de onda simplemente visuales con cierto grado de precisión.
Michael Stimson
Probablemente debería haber incluido más detalles sobre lo que tengo. Los conjuntos de datos sin procesar que tengo son: 1.Los datos completos y un ráster tiff generado a partir de ellos 2. Imágenes visuales (como la imagen de muestra que se muestra, pero que cubren un área mucho más amplia) 3. Mediciones directas manuales de un subconjunto de los árboles en la zona. De estos he generado: 1. las clasificaciones de suelo / vegetación 2. los rasters DEM / DSM
Theodore Jones
¿Tienes acceso a eCognition? Si no, ¿a qué software de procesamiento de imágenes o lenguajes de programación tiene acceso o conoce?
Aaron
No tengo una copia de eCognition, pero comprobaré si alguien que conozco en mi laboratorio / universidad la tiene porque parece popular para este tipo de cosas. Conozco Python, C y Java. Tengo una copia de Matlab pero soy un novato en eso. Tengo acceso a cualquiera de los software en esta lista softwarelicense.arizona.edu/students , además, por supuesto, ArcGIS.
Theodore Jones
Un poco más de detalle en las aplicaciones comerciales que tengo. Algunos de los que figuran en esa lista de software que he vinculado son Matlab, Mathematica, JMP y otras herramientas de estadísticas, y herramientas de desarrollo de software como Visual Studio.
Theodore Jones

Respuestas:

10

Existe una considerable cantidad de literatura sobre detección de corona individual en datos espectrales y lidar. Métodos sabios, tal vez comience con:

Falkowski, MJ, AMS Smith, PE Gessler, AT Hudak, LA Vierling y JS Evans. (2008) La influencia de la cubierta del dosel del bosque de coníferas en la precisión de dos algoritmos de medición de árboles individuales utilizando datos LIDAR. Canadian Journal of Remote Sensing 34 (2): 338-350.

Smith AMS, EK Strand, CM Steele, DB Hann, SR Garrity, MJ Falkowski, JS Evans (2008) Producción de mapas de estructura espacial de vegetación mediante análisis por objeto de la invasión de enebro en fotografías aéreas multitemporales. Canadian Journal Remote Sensing 34 (2): 268-285

Si está interesado en el método Wavelet (Smith et al., 2008), lo tengo codificado en Python, pero es muy lento. Si tiene experiencia con Matlab, aquí es donde se implementa en modo de producción. Tenemos dos documentos en los que identificamos ~ 6 millones de acres de invasión de enebro en el este de Oregón utilizando el método wavelet con imágenes NAIP RGB-NIR, por lo que está bien probado.

Baruch-Mordo, S., JS Evans, J. Severson, JD Naugle, J. Kiesecker, J. Maestas y MJ Falkowski (2013) Salvando el urogallo de los árboles: una solución proactiva para reducir una amenaza clave para un candidato Conservación biológica de especies 167: 233-241

Poznanovic, AJ, MJ Falkowski, AL Maclean y JS Evans (2014) Evaluación de precisión de algoritmos de detección de árboles en bosques de enebro. Ingeniería fotogramétrica y teledetección 80 (5): 627–637

Existen algunos enfoques interesantes, en general la descomposición de objetos, de la literatura de espacio de estado matemática aplicada que utiliza procesos gaussianos de resolución múltiple para descomponer las características del objeto a escala. Utilizo este tipo de modelos para describir procesos a múltiples escalas en modelos ecológicos, pero podría adaptarse para descomponer las características de los objetos de imagen. Divertido, pero un poco esotérico.

Gramacy, RB y HKH Lee (2008) treparon los modelos de proceso gaussianos bayesianos con una aplicación de modelado por computadora. Revista de la Asociación Americana de Estadística, 103 (483): 1119–1130

Kim, HM, BK Mallick y CC Holmes (2005) Análisis de datos espaciales no estacionarios utilizando procesos Gaussianos por partes. Revista de la Asociación Americana de Estadística, 100 (470): 653–668

Jeffrey Evans
fuente
+1 Especialmente para la opción 4; Dado que el OP tiene datos LIDAR, valdría la pena ejecutar el método wavelet en un modelo de superficie de dosel. Aunque, como saben, el método wavelet todavía no es corriente (o quizás nunca).
Aaron
En una oda al ideal de talla única, voy a comenzar a referirme al software comercial (p. Ej., ESRI, ERDAS) como software Big-box. A menudo, la mejor solución, o alguna, no está disponible en el "software Big-box". A menudo uno tiene que mirar hacia el desarrollo o las comunidades académicas en busca de respuestas a problemas analíticos espaciales complejos. Esto te saca de la corriente principal con mucha prisa. Afortunadamente, a estas comunidades les gusta compartir. Por eso también es importante que un analista no confíe en soluciones de botón.
Jeffrey Evans
2
Tiendo a estar de acuerdo sobre BBS para problemas espaciales complejos. Sin embargo, extraer un solo tipo de vegetación en un entorno árido, especialmente si tiene acceso a datos LIDAR, es bastante común. En este caso, no hay necesidad de reinventar la rueda desarrollando un nuevo enfoque para la identificación simple de árboles. Mis pensamientos son ¿por qué no utilizar un enfoque de botón preestablecido, especialmente en un paquete como eCognition, que es muy adecuado para la automatización?
Aaron
1
Debo agregar que eCognition tiene la capacidad de identificación de corona individual. Como ejemplo, puede encontrar un conjunto de reglas de muestra en la comunidad de eCog que utiliza un enfoque de crecimiento de semillas: busque "Conjunto de reglas de muestra de delineación de palmeras de aceite". La integración del nuevo algoritmo de coincidencia de plantillas de eCog y este enfoque de crecimiento de semillas podría ser un método muy poderoso.
Aaron
1
Estoy interesado en el código de Python que menciona para el método Wavelet de Smith (2008). ¿Está disponible en cualquier lugar?
Alpheus
3

Para crear un DHM, reste el DEM del DEM, esto se puede hacer en Esri Raster Calculator o GDAL_CALC . Esto pondrá todas tus elevaciones en un 'campo de juego nivelado'.

Sintaxis (Sustituya rutas completas para DEM, DSM y DHM):

GDAL_CALC.py -A DSM -B DEM --outfile=DHM --CALC "A-B"

El DHM será mayormente 0 (o lo suficientemente cerca), lo que hace que su valor de nodata. Con Raster Calculator o GDAL_CALC puede extraer valores más que un valor arbitrario en función de la cantidad de ruido que observe en el DHM. El objetivo de esto es reducir el ruido y resaltar solo las coronas de vegetación; en el caso de que dos 'árboles' estén adyacentes, esto debería dividirse en dos burbujas distintas.

Sintaxis (Sustituya las rutas completas para Binary y DHM y el valor observado para Value):

GDAL_CALC.py -A DHM --outfile=Binary --calc "A*(A>Value)"

Ahora, con GDAL_CALC o Esri IsNull, cree un ráster binario, que puede poligonizarse con GDAL_Polygonize o Esri Raster to Polygon .

Para refinar los polígonos, elimine los polígonos excesivamente pequeños y luego compárelos con las bandas RGB en busca de firmas, en Esri, la herramienta de Estadísticas zonales ayudará. Luego puede descartar los polígonos que claramente no tienen las estadísticas correctas (según la experimentación y sus datos, no puedo darle los valores).

Esto debería llevarlo a aproximadamente un 80% de precisión al trazar coronas individuales.

Michael Stimson
fuente
Gracias. Veré si obtengo buenos resultados con este método.
Theodore Jones
Tendrá que experimentar un poco para obtener los valores adecuados, sugiero recortar áreas pequeñas como muestras que son indicativas (similares a) las mejores / peores áreas en sus datos. Puede tomar media docena de ejecuciones de muestra para obtener sus parámetros, pero aún así es mejor que trazarlos manualmente.
Michael Stimson
3

eCognition es el mejor software para eso, lo hice usando otro software pero eCognition es mejor. Aquí está la referencia a la literatura sobre el tema:

Karlson, M., Reese, H. y Ostwald, M. (2014). Mapeo de copa de árbol en bosques gestionados (zonas verdes) de África occidental semiárida utilizando imágenes de WorldView-2 y análisis de imágenes basadas en objetos geográficos. Sensores, 14 (12), 22643-22669.

por ejemplo, http://www.mdpi.com/1424-8220/14/12/22643

Adicionalmente:

Zagalikis, G., Cameron, AD y Miller, DR (2005). La aplicación de la fotogrametría digital y las técnicas de análisis de imágenes para derivar las características de los árboles y los rodales. Revista canadiense de investigación forestal, 35 (5), 1224-1237.

por ejemplo, http://www.nrcresearchpress.com/doi/abs/10.1139/x05-030#.VJmMb14gAA

Giorgos Zagalikis
fuente
¿Podría explicar por qué eCognition es mejor? Las respuestas de solo enlace tienden a desaparecer cuando el enlace desaparece.
Aaron
1
eCognition es un software de análisis de imágenes basado en objetos que otros no son desde ahora. Usé un enfoque similar. La aplicación de técnicas de fotogrametría digital y análisis de imágenes para derivar las características de los árboles y los rodales G Zagalikis, AD Cameron, DR Miller Canadian Journal of Forest Research, 2005, 35 (5): 1224-1237, 10.1139 / x05-030 nrcresearchpress.com/doi /abs/10.1139/x05-030#.VJmMb14gAA
Giorgos Zagalikis
1
Gracias por la referencia Giorgos. Creo que estos comentarios funcionarían bien como una edición de su respuesta.
Aaron
3

Tuve el mismo problema hace un par de años. Tengo una solución que no requiere datos LAS filtrados u otros datos auxiliares. Si tiene acceso a datos LiDAR y puede generar DEM / DSM / DHM (DEM en adelante, no estoy debatiendo la semántica de la nomenclatura del modelo de superficie) a partir de diferentes retornos, el siguiente script puede ser útil.

El script arcpy ingiere 3 DEM y escupe un polígono forestal y archivos de formas de puntos de árbol. Los 3 DEM deben tener la misma resolución espacial (es decir, 1 metro) y extensiones, y representan los primeros retornos, los últimos retornos y la tierra desnuda. Tenía parámetros muy específicos para la extracción de verduras, pero los parámetros pueden modificarse para adaptarse a otras necesidades. Estoy seguro de que el proceso se puede mejorar, ya que este fue mi primer intento serio de secuencias de comandos de Python.

# Name:         Veg_Extractor.py
# Date:         2013-07-16
# Usage:        ArcMap 10.0; Spatial Analyst
# Input:        1 meter DEMs for first returns (DEM1), last returns (DEM2), and bare earth (BE)
# Output:       forest polygon (veg with height > 4m) shapefile with holes > 500m removed;
#               tree point (veg with height > 4m, crown radius of 9 cells) shapefile
# Notes:        Raises error if input raster cell sizes differ

import arcpy, os
from arcpy import env
from arcpy.sa import *

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Script arguments
dem1 = arcpy.GetParameterAsText(0) #input Raster Layer, First Return DEM
dem2 = arcpy.GetParameterAsText(1) #input Raster Layer, Last Return DEM
bare_earth = arcpy.GetParameterAsText(2) #input Raster Layer, Bare Earth DEM
outForest = arcpy.GetParameterAsText(3) #shapefile
outTree = arcpy.GetParameterAsText(4) #shapefile

# Make sure cell size of input rasters are same
arcpy.AddMessage("Checking cell sizes...")
dem1Xresult = arcpy.GetRasterProperties_management(dem1, "CELLSIZEX")
dem1Yresult = arcpy.GetRasterProperties_management(dem1, "CELLSIZEY")
dem2Xresult = arcpy.GetRasterProperties_management(dem2, "CELLSIZEX")
dem2Yresult = arcpy.GetRasterProperties_management(dem2, "CELLSIZEY")
beXresult = arcpy.GetRasterProperties_management(bare_earth, "CELLSIZEX")
beYresult = arcpy.GetRasterProperties_management(bare_earth, "CELLSIZEY")
dem1X = round(float(dem1Xresult.getOutput(0)),4)
dem1Y = round(float(dem1Yresult.getOutput(0)),4)
dem2X = round(float(dem2Xresult.getOutput(0)),4)
dem2Y = round(float(dem2Yresult.getOutput(0)),4)
beX = round(float(beXresult.getOutput(0)),4)
beY = round(float(beYresult.getOutput(0)),4)
if (dem1X == dem1Y == dem2X == dem2Y == beX == beY) == True:
    arcpy.AddMessage("Cell sizes match.")
else:
    arcpy.AddMessage("Input Raster Cell Sizes:")
    arcpy.AddMessage("DEM1: (" + str(dem1X) + "," + str(dem1Y) + ")")
    arcpy.AddMessage("DEM2: (" + str(dem2X) + "," + str(dem2Y) + ")")
    arcpy.AddMessage("  BE: (" + str(beX) + "," + str(beY) + ")")
    raise Exception("Cell sizes do not match.")

# Check map units
dem1_spatial_ref = arcpy.Describe(dem1).spatialReference
dem1_units = dem1_spatial_ref.linearUnitName
dem2_spatial_ref = arcpy.Describe(dem2).spatialReference
dem2_units = dem2_spatial_ref.linearUnitName
bare_earth_spatial_ref = arcpy.Describe(bare_earth).spatialReference
bare_earth_units = bare_earth_spatial_ref.linearUnitName
if (dem1_units == dem2_units == bare_earth_units) == True:
    if dem1_units == "Meter":
        area = "500 SquareMeters" #Area variable for meter
        unit = 1 #meter
    elif (dem1_units == "Foot_US") or (dem1_units == "Foot"):
        area = "5382 SquareFeet" #Area variable for feet
        unit = 3.28084 #feet in meter
    else:
        raise Exception("Units are not 'Meter', 'Foot_US', or 'Foot'.")
else:
    raise Exception("Linear units do not match.  Check spatial reference.")

# Local variables:
(workspace, filename) = os.path.split(outForest)
arcpy.env.workspace = workspace
arcpy.env.overwriteOutput = True
dem1 = Raster(dem1)
dem2 = Raster(dem2)
bare_earth = Raster(bare_earth)
nbr1 = NbrRectangle(3, 3, "CELL")
nbr2 = NbrRectangle(5, 5, "CELL")
nbr3 = NbrCircle(5, "CELL")

# Give units and multiplier
arcpy.AddMessage("Linear units are " + dem1_units + ". Using multiplier of " + str(unit) + "...")

arcpy.AddMessage("Processing DEMs...")
# Process: Raster Calculator (DEM1 - BE)
ndsm_dem1 = dem1 - bare_earth

# Process: Raster Calculator (DEM1 - DEM2)
d1_d2 = dem1 - dem2

# Process: Raster Calculator
threshold_d1d2 = (d1_d2 > (0.1 * unit))  &  (ndsm_dem1 >= (4.0 * unit))

# Process: Focal Statistics (max 3x3)
focal_max1 = FocalStatistics(threshold_d1d2, nbr1, "MAXIMUM", "DATA")

# Process: Focal Statistics (majority 5x5)
focal_majority = FocalStatistics(focal_max1, nbr2, "MAJORITY", "DATA")

# Process: Con
con_ndsm_dem1 = Con(ndsm_dem1 >= (4.0 * unit), focal_majority, focal_max1)
focal_majority = None
focal_max1 = None

# Process: Focal Statistics (min 3x3)
focal_min1 = FocalStatistics(con_ndsm_dem1, nbr1, "MINIMUM", "DATA")
con_ndsm_dem1 = None

# Process: Focal Statistics (min 3x3)
veg_mask = FocalStatistics(focal_min1, nbr1, "MINIMUM", "DATA")

# Process: Focal Statistics (max R5)
focal_max2 = FocalStatistics(ndsm_dem1, nbr3, "MAXIMUM", "DATA")

arcpy.AddMessage("Calculating tree points...")
# Process: Raster Calculator
tree_points = (veg_mask == 1) & (ndsm_dem1 == focal_max2) & (ndsm_dem1 >= (4.0 * unit))
ndsm_dem1 = None
focal_max2 = None

# Process: Raster Calculator
tree_pick = Pick(tree_points == 1, 1)
tree_points = None

# Process: Raster to Polygon
arcpy.RasterToPolygon_conversion(tree_pick, workspace + "\\tree_poly.shp", "SIMPLIFY", "Value")
tree_pick = None

# Process: Feature To Point
arcpy.AddMessage("Writing tree points...")
arcpy.env.workspace = workspace #reset workspace
arcpy.env.overwriteOutput = True #reset overwrite permission
arcpy.FeatureToPoint_management(workspace + "\\tree_poly.shp", outTree, "CENTROID")

arcpy.AddMessage("Calculating forest polygons...")
# Process: Focal Statistics (max 3x3)
forests = FocalStatistics(veg_mask, nbr1, "MAXIMUM", "DATA")
veg_mask = None

# Process: Raster Calculator
forest_pick = Pick(forests == 1, 1)

# Process: Raster to Polygon
arcpy.RasterToPolygon_conversion(forest_pick, workspace + "\\forest_poly.shp", "SIMPLIFY", "Value")

# Process: Eliminate Holes > 500 sq m (5382 sq ft)
arcpy.AddMessage("Writing forest polygons...")
arcpy.EliminatePolygonPart_management(workspace + "\\forest_poly.shp", outForest, "AREA", area, "0", "CONTAINED_ONLY")

# Clean up
arcpy.AddMessage("Cleaing up...")
arcpy.Delete_management(workspace + "\\tree_poly.shp")
arcpy.Delete_management(workspace + "\\forest_poly.shp")
Barbarroja
fuente
2

Estoy publicando esto como respuesta debido al límite de longitud en el comentario, sin esperanzas de créditos :). Cepillo muy amplio, siempre que tenga DEM.

  1. Extraer DEM para polígono individual a dem.
  2. Define los extremos de elevación de dem
  3. Establezca zCur + = - zStep. Paso que se debe encontrar por iteraciones de antemano, por ejemplo, caída razonable entre la elevación de la 'celda superior del árbol' y los vecinos
  4. Abajo = Con (dem => zCur, int (1))
  5. Agrupe las regiones de abajo. Cuenta lo suficientemente grande, que son 'árboles'. Definición requerida aquí por inspección visual, investigación preliminar?
  6. Vaya al paso 3 si zCur> zMin, paso 1 de lo contrario.

Número máximo de grupos en el proceso = recuento de árboles dentro de un polígono individual. Criterios adicionales, por ejemplo, la distancia entre 'árboles' dentro de los polígonos podría ayudar ... El suavizado DEM usando el núcleo también es una opción.

FelixIP
fuente
Creo que te estás refiriendo a un DSM y no a un DEM ... Por lo general, los árboles, las estructuras y otros elementos no deseados no se convierten en un DEM, sino que aparecen en el DSM (menos clases de ruido). DSM - DEM = DHM (modelo de altura). Todo esto se puede extraer razonablemente de los datos LiDAR, incluso si solo se clasifica tierra / no tierra, pero si tiene el DEM y no el LAS, está en el arroyo sin una paleta porque las características que busca no son ahi !
Michael Stimson
Yeap, DHM como lo describiste lo hará. Sé poco sobre Lidar.
FelixIP