Extraction des zones de cime des arbres à partir des données de télédétection (images visuelles et LiDAR)

Je recherche une méthode pour traiter une image de télédétection et extraire les zones de cime des arbres individuels de l'image.

J'ai à la fois des images surfaciques de longueur d'onde visuelle et des données lidar de la région. L'endroit en question est une zone désertique, donc la couverture arborée n'est pas aussi dense qu'une zone forestière. La résolution de l'imagerie aérienne est de 0,5 pied par 0,5 pied. La résolution lidar est d'environ 1 x 1 pied. Les données visuelles et le lidar proviennent d'un ensemble de données du comté de Pima, en Arizona. Un échantillon du type d'imagerie aérienne que j'ai est à la fin de ce post.

Cette question Détection d'arbre unique dans ArcMap? semble être le même problème, mais il ne semble pas y avoir de bonne réponse.

Je peux obtenir une classification raisonnable des types de végétation (et des informations sur le pourcentage global de couverture) dans la zone en utilisant la classification Iso Cluster à Arcmap, mais cela fournit peu d'informations sur les arbres individuels. Le plus proche de ce que je veux est le résultat du passage de la sortie de la classification d'isocluster via la fonction Raster vers polygone dans Arcmap. Le problème est que cette méthode fusionne près des arbres en un seul polygone.

Edit: j'aurais probablement dû inclure plus de détails sur ce que j'ai. Les ensembles de données brutes que j'ai sont:

• Données las complètes et raster tiff généré à partir de celles-ci.
• Imagerie visuelle (comme l'exemple d'image illustré, mais couvrant une zone beaucoup plus large)
• Mesures directes manuelles d'un sous-ensemble d'arbres de la zone.

À partir de ceux-ci, j'ai généré:

1. Les classifications sol / végétation.
2. Les rasters DEM / DSM.

Il existe une littérature considérable sur la détection des couronnes individuelles dans les données spectrales et lidar. En ce qui concerne les méthodes, commencez peut-être par:

Falkowski, MJ, AMS Smith, PE Gessler, AT Hudak, LA Vierling et JS Evans. (2008). L'influence de la couverture du couvert forestier de conifères sur la précision de deux algorithmes de mesure d'arbres individuels utilisant des données lidar. Journal canadien de télédétection 34 (2): 338-350.

Smith AMS, EK Strand, CM Steele, DB Hann, SR Garrity, MJ Falkowski, JS Evans (2008) Production de cartes de la structure spatiale de la végétation par analyse par objet de l'empiètement du genévrier dans les photographies aériennes multi-temporelles. Canadian Journal Remote Sensing 34 (2): 268-285

Si vous êtes intéressé par la méthode des ondelettes (Smith et al., 2008), je l'ai codée en Python mais, elle est très lente. Si vous avez une expérience Matlab, c'est là qu'il est implémenté en mode production. Nous avons deux articles où nous avons identifié environ 6 millions d'acres d'empiètement de genévriers dans l'est de l'Oregon en utilisant la méthode des ondelettes avec l'imagerie NAIP RGB-NIR, donc c'est bien prouvé.

Baruch-Mordo, S., JS Evans, J. Severson, JD Naugle, J. Kiesecker, J. Maestas et MJ Falkowski (2013) Sauver le tétras des arbres: une solution proactive pour réduire une menace clé pour un candidat espèces Conservation biologique 167: 233-241

Poznanovic, AJ, MJ Falkowski, AL Maclean et JS Evans (2014) An Accuracy Assessment of Tree Detection Algorithms in Juniper Woodlands. Ingénierie photogrammétrique et télédétection 80 (5): 627–637

Il existe des approches intéressantes, dans la décomposition générale des objets, à partir de la littérature sur l'espace d'état des mathématiques appliquées utilisant des processus gaussiens multirésolutions pour décomposer les caractéristiques des objets à l'échelle. J'utilise ces types de modèles pour décrire des processus multi-échelles dans des modèles écologiques mais il pourrait être adapté pour décomposer les caractéristiques des objets image. Amusant, mais un peu ésotérique.

Gramacy, RB et HKH Lee (2008) Modèles de processus gaussiens arborés bayésiens avec une application à la modélisation informatique. Journal de l'American Statistical Association, 103 (483): 1119–1130

Kim, HM, BK Mallick et CC Holmes (2005) Analyse de données spatiales non stationnaires à l'aide de processus gaussiens par morceaux. Journal de l'American Statistical Association, 100 (470): 653–668

Pour créer un DHM, soustrayez le DEM du DEM, cela peut être fait dans Esri Raster Calculator ou GDAL_CALC . Cela mettra toutes vos élévations sur un «terrain de jeu égal».

Syntaxe (Remplacer les chemins d'accès complets pour DEM, DSM et DHM):

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

Le DHM sera principalement 0 (ou assez proche), ce que vous faites votre valeur nodata. Avec Raster Calculator ou GDAL_CALC, vous pouvez extraire des valeurs plus qu'une valeur arbitraire en fonction de la quantité de bruit que vous observez dans le DHM. Le but de ceci est de réduire le bruit et de mettre en évidence uniquement les couronnes de végétation - dans le cas où deux «arbres» sont adjacents, cela devrait se diviser en deux taches distinctes.

Syntaxe (Substituez les chemins complets pour Binary & DHM et la valeur observée pour Value):

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

Maintenant, avec GDAL_CALC ou Esri IsNull, créez un raster binaire, qui peut être polygonisé avec GDAL_Polygonize ou Esri Raster to Polygon .

Pour affiner les polygones, supprimez les polygones excessivement petits, puis comparez-les aux bandes RVB à la recherche de signatures. Dans Esri, l' outil Statistiques zonales vous aidera. Ensuite, vous pouvez éliminer les polygones qui n'ont clairement pas les bonnes statistiques (sur la base de l'expérimentation et de vos données, je ne peux pas vous donner les valeurs).

Cela devrait vous permettre d'atteindre une précision d'environ 80% lors du traçage des couronnes individuelles.

eCognition est le meilleur logiciel pour ça, je l'ai fait en utilisant d'autres logiciels mais eCognition c'est mieux. Voici la référence à la littérature sur le sujet:

Karlson, M., Reese, H., et Ostwald, M. (2014). Cartographie des cimes des arbres dans les zones boisées gérées (parcs) de l'Afrique de l'Ouest semi-aride à l'aide de l'imagerie WorldView-2 et de l'analyse d'images basée sur des objets géographiques. Capteurs, 14 (12), 22643-22669.

par exemple http://www.mdpi.com/1424-8220/14/12/22643

Aditionellement:

Zagalikis, G., Cameron, AD et Miller, DR (2005). L'application de techniques de photogrammétrie numérique et d'analyse d'images pour dériver les caractéristiques des arbres et des peuplements. Revue canadienne de recherche forestière, 35 (5), 1224-1237.

par exemple http://www.nrcresearchpress.com/doi/abs/10.1139/x05-030#.VJmMb14gAA

J'ai eu le même problème il y a quelques années. J'ai une solution qui ne nécessite pas de données LAS filtrées ou d'autres données auxiliaires. Si vous avez accès aux données LiDAR et pouvez générer des DEM / DSM / DHM (DEM ci-après, je ne débat pas de la sémantique de la nomenclature des modèles de surface) à partir de différents retours, le script suivant peut être utile.

Le script arcpy ingère 3 DEM et crache un polygone de forêt et des fichiers de formes de point d'arbre. Les 3 DEM devraient avoir la même résolution spatiale (c.-à-d. 1 mètre) et les mêmes étendues, et représenter les premiers retours, les derniers retours et la terre nue. J'avais des paramètres très spécifiques pour l'extraction des légumes, mais les paramètres peuvent être modifiés pour répondre à d'autres besoins. Je suis sûr que le processus peut être amélioré, car c'était ma première tentative sérieuse de scripting 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")

Je poste ceci comme réponse en raison de la limite de longueur en commentaire, aucun espoir de crédits :). Pinceau très large, à condition d'avoir DEM.

1. Extraire DEM pour un polygone individuel à dém.
2. Définir les extrêmes d'élévation de dem
3. Définissez zCur + = - zStep. Étape à trouver par les itérations au préalable, par exemple, chute raisonnable entre l'élévation de la «cellule supérieure de l'arbre» ​​et les voisins
4. En dessous = Con (dem => zCur, int (1))
5. Regroupez les régions ci-dessous. Comptez assez gros, ce sont des «arbres». Définition requise ici par inspection visuelle, recherche préliminaire?
6. Passez à l'étape 3 si zCur> zMin, étape 1 sinon.

Nombre maximum de groupes dans le processus = nombre d'arbres dans un polygone individuel. Des critères supplémentaires, par exemple la distance entre les «arbres» à l'intérieur des polygones pourraient aider ... Le lissage DEM à l'aide du noyau est également une option.