Extrahieren von Baumkronenbereichen aus Fernerkundungsdaten (visuelle Bilder und LiDAR)

12

Ich suche nach einer Methode, um ein Fernerkundungsbild zu verarbeiten und die Kronenbereiche der einzelnen Bäume aus dem Bild zu extrahieren.

Ich habe sowohl visuelle Wellenlängen-Flächenbilder als auch Lidar-Daten aus der Fläche. Bei dem fraglichen Ort handelt es sich um ein Wüstengebiet, daher ist die Baumdecke nicht so dicht wie ein Waldgebiet. Die Auflösung der Luftbilder beträgt 0,5 mal 0,5 Fuß. Die Lidar-Auflösung beträgt ungefähr 1 x 1 Fuß. Sowohl die visuellen Daten als auch das Lidar stammen aus einem Datensatz von Pima County, Arizona. Ein Beispiel für die Art der Luftbilder, die ich habe, befindet sich am Ende dieses Beitrags.

Diese Frage Einzelbaumerkennung in ArcMap? scheint das gleiche Problem zu sein, aber es scheint keine gute Antwort zu geben.

Mit der Iso-Cluster-Klassifizierung in Arcmap kann ich eine vernünftige Klassifizierung der Vegetationstypen (und Informationen zur prozentualen Gesamtbedeckung) in dem Gebiet erhalten, die jedoch nur wenige Informationen zu einzelnen Bäumen enthält. Das Ergebnis, das dem gewünschten Ergebnis am nächsten kommt, ist das Übergeben der Ausgabe der Isocluster-Klassifizierung durch das Feature "Raster zu Polygon" in Arcmap. Das Problem ist, dass diese Methode in der Nähe von Bäumen zu einem einzigen Polygon verschmilzt.

Edit: Ich hätte wahrscheinlich etwas detaillierter darüber berichten sollen, was ich habe. Die Rohdatensätze, die ich habe, sind:

  • Volle Laserdaten und ein daraus generiertes TIFF-Raster.
  • Bildmaterial (wie das gezeigte Beispielbild, jedoch mit einem viel größeren Bereich)
  • Manuelle direkte Messungen einer Teilmenge der Bäume in der Region.

Daraus habe ich generiert:

  1. Die Boden- / Vegetationsklassifikationen.
  2. Die DEM / DSM-Raster.

Beispiel Luftbilder

Theodore Jones
quelle
Sie haben mehr Daten als der Link. Haben Sie die klassifizierten Las-Dateien oder nur das DEM / DSM-Raster (welches?)? Es ist wirklich nicht einfach, dies nur mit visuellen Wellenlängen mit beliebiger Genauigkeit zu tun.
Michael Stimson
Ich hätte wahrscheinlich etwas detaillierter darüber berichten sollen, was ich habe. Die Rohdatensätze, die ich habe, sind: 1. Vollständige Laserdaten und ein daraus erzeugtes TIFF-Raster 2. Visuelle Bilder (wie das gezeigte Beispielbild, das jedoch einen viel größeren Bereich abdeckt) 3. Manuelle direkte Messungen einer Teilmenge der Bäume in das Gebiet. Daraus habe ich generiert: 1. die Boden- / Vegetationsklassifikationen 2. die DEM- / DSM-Raster
Theodore Jones
Haben Sie Zugang zu eCognition? Wenn nicht, auf welche Bildverarbeitungssoftware oder Programmiersprachen haben Sie Zugriff oder kennen sie?
Aaron
Ich habe keine Kopie von eCognition, aber ich überprüfe, ob jemand, den ich in meinem Labor / an meiner Universität kenne, eine Kopie davon hat, weil sie für diese Art von Dingen beliebt zu sein scheint. Ich kenne mich in Python, C und Java aus. Ich habe eine Kopie von Matlab, aber ich bin so ziemlich ein Noob. Ich habe Zugriff auf die Software auf dieser Liste softwarelicense.arizona.edu/students und natürlich auch auf ArcGIS.
Theodore Jones
Ein bisschen detaillierter in den kommerziellen Anwendungen, die ich habe. Einige der auf dieser Liste der von mir verknüpften Software sind Matlab, Mathematica, JMP und andere Statistik-Tools sowie Software-Entwicklungstools wie Visual Studio.
Theodore Jones

Antworten:

10

Es gibt eine beträchtliche Menge an Literatur zur Einzelkronendetektion in Spektral- und Lidardaten. Methoden weise beginnen vielleicht mit:

Falkowski, MJ, AMS Smith, PE Gessler, AT Hudak, LA Vierling und JS Evans. (2008). Der Einfluss der Nadelwaldüberdachung auf die Genauigkeit von zwei einzelnen Baummessalgorithmen unter Verwendung von Lidardaten. Canadian Journal of Remote Sensing 34 (2): 338-350.

Smith AMS, EK Strand, CM Steele, DB Hann, SR Garrity, MJ Falkowski, JS Evans (2008) Erstellung von Vegetationsraumstrukturkarten durch objektbezogene Analyse des Wacholderangriffs in mehrzeitigen Luftbildern. Canadian Journal Remote Sensing 34 (2): 268-285

Wenn Sie sich für die Wavelet-Methode interessieren (Smith et al., 2008), habe ich sie in Python codiert, aber sie ist sehr langsam. Wenn Sie Matlab-Erfahrung haben, wird es hier im Produktionsmodus implementiert. Wir haben zwei Veröffentlichungen, in denen wir nach der Wavelet-Methode mit NAIP-RGB-NIR-Bildern ~ 6 Millionen Morgen Wacholder-Übergriffe in Ost-Oregon identifiziert haben.

S. Baruch-Mordo, J. S. Evans, J. Severson, J. D. Naugle, J. Kiesecker, J. Maestas und MJ Falkowski (2013) Rettung von Salbei vor den Bäumen: Eine proaktive Lösung zur Reduzierung einer Schlüsselbedrohung für einen Kandidaten species Biological Conservation 167: 233 & ndash; 241

Poznanovic, AJ, MJ Falkowski, AL Maclean und JS Evans (2014) Eine Genauigkeitsbewertung von Baumerkennungsalgorithmen in Juniper Woodlands. Photogrammetric Engineering & Remote Sensing 80 (5): 627–637

Es gibt einige interessante Ansätze zur allgemeinen Objektzerlegung aus der Literatur zum angewandten mathematischen Zustandsraum unter Verwendung von Gaußschen Prozessen mit mehreren Auflösungen, um Objekteigenschaften über den Maßstab hinweg zu zerlegen. Ich benutze diese Arten von Modellen, um Prozesse auf mehreren Skalen in ökologischen Modellen zu beschreiben, aber es könnte angepasst werden, um die Eigenschaften von Bildobjekten zu zerlegen. Spaß, aber ein bisschen esoterisch.

Gramacy, RB und HKH Lee (2008) Bayesian treed Gaussian Prozessmodelle mit einer Anwendung zur Computermodellierung. Journal of the American Statistical Association, 103 (483): 1119–1130

Kim, HM, BK Mallick und CC Holmes (2005) analysieren nichtstationäre räumliche Daten mit stückweisen Gaußschen Prozessen. Journal of the American Statistical Association, 100 (470): 653–668

Jeffrey Evans
quelle
+1 Speziell für Option 4; Da das OP LIDAR-Daten enthält, ist es empfehlenswert, die Wavelet-Methode auf einem Baldachin-Oberflächenmodell auszuführen. Obwohl, wie Sie wissen, die Wavelet-Methode noch nicht wirklich Mainstream ist (oder vielleicht jemals).
Aaron
In einer Ode an das Ideal der Einheitsgröße werde ich anfangen, kommerzielle Software (z. B. ESRI, ERDAS) als Big-Box-Software zu bezeichnen. Oft ist die beste oder überhaupt keine Lösung in der "Big-Box-Software" verfügbar. Oft muss man sich an die Entwicklungs- oder Akademikergemeinschaften wenden, um Antworten auf komplexe räumliche Analyseprobleme zu finden. Dies bringt Sie in großer Eile aus dem Mainstream. Zum Glück teilen diese Gemeinschaften gerne. Aus diesem Grund ist es für einen Analysten wichtig, sich nicht auf Druckknopflösungen zu verlassen.
Jeffrey Evans
2
Über BBS stimme ich bei komplexen räumlichen Problemen eher zu. Das Extrahieren eines einzelnen Vegetationstyps in einer trockenen Umgebung - insbesondere wenn Sie Zugriff auf Lidar-Daten haben - ist jedoch ein ziemlich normaler Vorgang. In diesem Fall ist es nicht erforderlich, das Rad neu zu erfinden, indem ein neuer Ansatz zur einfachen Baumidentifizierung entwickelt wird. Meine Gedanken sind, warum nicht einen vorgefertigten Druckknopf-Ansatz verwenden, insbesondere in einem Paket wie eCognition, das sich hervorragend für die Automatisierung eignet?
Aaron
1
Ich sollte hinzufügen, dass eCognition die Fähigkeit zur individuellen Kronen-ID besitzt. Als Beispiel finden Sie in der eCog-Community einen Beispiel-Regelsatz, bei dem ein Saatgutwachstumsansatz verwendet wird. Suchen Sie nach "Beispiel-Regelsatz zur Abgrenzung von Ölpalmen". Die Integration des neuen Template-Matching-Algorithmus von eCog und dieses Ansatzes zur Saatgutvermehrung könnte möglicherweise eine sehr leistungsfähige Methode sein.
Aaron
1
Ich interessiere mich für den Python-Code, den Sie für die Smith (2008) Wavelet-Methode erwähnen. Ist es überall verfügbar?
Alpheus
4

eCognition ist die beste Software dafür, ich habe das mit einer anderen Software gemacht, aber eCognition ist besser. Hier ist der Verweis auf Literatur zu diesem Thema:

Karlson, M., Reese, H. & Ostwald, M. (2014). Baumkronenkartierung in bewirtschafteten Wäldern (Parklands) in halbtrockenem Westafrika mit WorldView-2-Bildern und geografischer objektbasierter Bildanalyse. Sensors, 14 (12), 22643 & ndash; 22669.

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

Zusätzlich:

Zagalikis, G., Cameron, AD & Miller, DR (2005). Die Anwendung digitaler Photogrammetrie- und Bildanalysetechniken zur Ableitung von Baum- und Bestandsmerkmalen. Canadian Journal of Forest Research, 35 (5), 1224-1237.

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

Giorgos Zagalikis
quelle
Könnten Sie bitte erläutern, warum eCognition besser ist? Antworten, die nur mit Links beantwortet werden, sind in der Regel nicht mehr funktionsfähig, wenn der Link nicht mehr funktioniert.
Aaron
1
eCognition ist eine objektbasierte Bildanalyse-Software, andere gibt es seitdem nicht mehr. Ich habe einen ähnlichen Ansatz gewählt. Anwendung digitaler Photogrammetrie- und Bildanalysetechniken zur Ableitung von Baum- und Bestandsmerkmalen 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
Danke für den Hinweis Giorgos. Ich denke, dass diese Kommentare als Bearbeitung Ihrer Antwort gut funktionieren würden.
Aaron
3

Um ein DHM zu erstellen, subtrahieren Sie das DEM vom DEM. Dies kann in Esri Raster Calculator oder GDAL_CALC erfolgen . Dadurch werden alle Ihre Erhöhungen auf ein „Level Playing Field“ gestellt.

Syntax (Vollständige Pfade für DEM, DSM und DHM ersetzen):

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

Das DHM ist meistens 0 (oder nahe genug), was Sie zu Ihrem Nodata-Wert machen. Mit Raster Calculator oder GDAL_CALC können Sie Werte extrahieren, die mehr als ein willkürlicher Wert sind, basierend auf der Menge an Rauschen, die Sie im DHM beobachten. Ziel ist es, das Rauschen zu reduzieren und nur die Vegetationskronen hervorzuheben - in dem Fall, in dem zwei „Bäume“ benachbart sind, sollte dies in zwei verschiedene Flecken aufgeteilt werden.

Syntax (Ersetzen Sie Binary & DHM durch vollständige Pfade und Value durch den beobachteten Wert):

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

Erstellen Sie jetzt mit GDAL_CALC oder Esri IsNull ein binäres Raster, das mit GDAL_Polygonize oder Esri Raster to Polygon polygonisiert werden kann .

Um die Polygone zu verfeinern, entfernen Sie zu kleine Polygone und vergleichen Sie sie mit den RGB-Bändern, um nach Signaturen zu suchen. In Esri hilft das Zonal Statistics-Tool . Dann können Sie die Polygone verwerfen, die eindeutig nicht die richtigen Statistiken haben (basierend auf Experimenten und Ihren Daten kann ich Ihnen die Werte nicht geben).

Dadurch sollten Sie beim Plotten einzelner Kronen eine Genauigkeit von ca. 80% erreichen.

Michael Stimson
quelle
Vielen Dank. Ich werde sehen, ob ich mit dieser Methode gute Ergebnisse erziele.
Theodore Jones
Sie müssen einige Experimente durchführen, um die geeigneten Werte zu erhalten. Ich empfehle, kleine Bereiche als Beispiele zu beschneiden, die die besten / schlechtesten Bereiche in Ihren Daten anzeigen (ähnlich). Es kann ein halbes Dutzend Probeläufe dauern, um Ihre Parameter zu erhalten, aber das muss immer noch besser sein, als sie manuell zu zeichnen.
Michael Stimson
3

Ich hatte vor ein paar Jahren das gleiche Problem. Ich habe eine Lösung, die keine gefilterten LAS-Daten oder andere Zusatzdaten erfordert. Wenn Sie Zugriff auf LiDAR-Daten haben und DEMs / DSMs / DHMs (im Folgenden: DEM) aus verschiedenen Rückgaben generieren können, ist das folgende Skript möglicherweise hilfreich.

Das arcpy-Skript nimmt 3 DEMs auf und spuckt ein Waldpolygon und Baumpunkt-Shapefiles aus. Die 3 DEMs sollten die gleiche räumliche Auflösung (dh 1 Meter) und Ausdehnung haben und die ersten Renditen, die letzten Renditen und die nackte Erde darstellen. Ich hatte sehr spezifische Parameter für die Gemüseextraktion, aber die Parameter können an andere Bedürfnisse angepasst werden. Ich bin sicher, dass der Prozess verbessert werden kann, da dies mein erster ernsthafter Versuch war, Python-Skripte zu erstellen.

# 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")
Barbarossa
quelle
2

Ich poste dies als Antwort wegen Längenbeschränkung im Kommentar, keine Hoffnung auf Credits :). Sehr breiter Pinsel, vorausgesetzt du hast DEM.

  1. Extrahiere DEM für ein einzelnes Polygon, um es zu demen.
  2. Definieren Sie die Elevationsextreme von dem
  3. Setze zCur + = - zStep. Schritt, der vorher durch Iterationen gefunden werden muss, z. B. angemessener Abstand zwischen der Höhe der obersten Baumzelle und den Nachbarn
  4. Below = Con (dem => zCur, int (1))
  5. Gruppenregionen von Below. Zählen Sie groß genug, das sind "Bäume". Definition hier erforderlich durch Sichtprüfung, Vorrecherche?
  6. Fahren Sie mit Schritt 3 fort, wenn zCur> zMin, andernfalls mit Schritt 1.

Maximale Anzahl der Gruppen im Prozess = Anzahl der Bäume innerhalb eines einzelnen Polygons. Zusätzliche Kriterien, z. B. der Abstand zwischen 'Bäumen' in Polygonen, könnten helfen ... Eine DEM-Glättung mit dem Kernel ist ebenfalls eine Option.

FelixIP
quelle
Ich glaube, Sie beziehen sich auf ein DSM und nicht auf ein DEM ... Normalerweise werden Bäume, Strukturen und anderer Müll nicht zu einem DEM, sondern zu einer Funktion in DSM (minus Rauschklassen). DSM - DEM = DHM (Höhenmodell). Alle diese Daten können vernünftigerweise aus LiDAR-Daten extrahiert werden, auch wenn sie nur als Grund / Nichtgrund klassifiziert sind. Wenn Sie jedoch das DEM und nicht das LAS haben, sind Sie ohne Paddel auf dem Laufenden, da die gewünschten Funktionen nicht vorhanden sind da !
Michael Stimson
Ja, DHM, wie du beschrieben hast. Ich weiß wenig über Lidar.
FelixIP