Wyodrębnianie obszarów koron drzew z danych teledetekcyjnych (obrazy wizualne i LiDAR)

Szukam metody przetwarzania obrazu za pomocą teledetekcji i wydobywania z obrazu obszarów koron poszczególnych drzew.

Mam zarówno wizualne zdjęcia powierzchniowe fal, jak i dane lidarowe z tego obszaru. Lokalizacja, o której mowa, to obszar pustynny, więc pokrywa drzew nie jest tak gęsta jak obszar leśny. Rozdzielczość zdjęć lotniczych wynosi 0,5 stopy na 0,5 stopy. Rozdzielczość lidar wynosi około 1 x 1 stóp. Zarówno dane wizualne, jak i lidar pochodzą z zestawu danych hrabstwa Pima w Arizonie. Próbka rodzaju zdjęć lotniczych, które mam, znajduje się na końcu tego postu.

To pytanie Wykrywanie pojedynczego drzewa w ArcMap? wydaje się być tym samym problemem, ale wydaje się, że nie ma tam dobrej odpowiedzi.

Mogę uzyskać rozsądną klasyfikację rodzajów roślinności (i informacje o ogólnej procentowej powierzchni pokrycia) na danym obszarze, stosując klasyfikację klastrów Iso w Arcmap, ale to daje niewiele informacji na temat poszczególnych drzew. Muszę najbardziej zbliżyć się do wyników przekazywania wyników klasyfikacji isocluster przez funkcję Raster to Polygon w Arcmap. Problem polega na tym, że ta metoda łączy drzewa w jednym wielokącie.

Edycja: Prawdopodobnie powinienem podać więcej szczegółów na temat tego, co mam. Surowe zestawy danych, które mam, to:

• Pełne dane lasów i wygenerowany z nich tiff raster.
• Obrazy wizualne (jak pokazano przykładowy obraz, ale obejmujące znacznie szerszy obszar)
• Ręczne bezpośrednie pomiary podzbioru drzew w okolicy.

Z tych wygenerowałem:

1. Klasyfikacje gruntów / roślinności.
2. Rastry DEM / DSM.

Istnieje znaczna literatura na temat wykrywania poszczególnych koron w danych spektralnych i lidarowych. Metody mądre, być może zacznij od:

Falkowski, MJ, AMS Smith, PE Gessler, AT Hudak, LA Vierling i JS Evans. (2008). Wpływ pokrywy drzewostanu iglastego na dokładność dwóch indywidualnych algorytmów pomiarowych drzew z wykorzystaniem danych 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) Produkcja map struktury przestrzennej roślinności poprzez analizę poszczególnych obiektów wkroczenia jałowca na wielorasowe zdjęcia lotnicze. Canadian Journal Remote Sensing 34 (2): 268–285

Jeśli interesuje Cię metoda Wavelet (Smith i in., 2008), mam ją zakodowaną w Pythonie, ale jest bardzo wolna. Jeśli masz doświadczenie w Matlabie, jest to miejsce, w którym jest ono implementowane w trybie produkcyjnym. Mamy dwa artykuły, w których zidentyfikowaliśmy około 6 milionów akrów wkroczenia jałowca we wschodnim Oregonie przy użyciu metody falkowej ze zdjęciami NAIP RGB-NIR, więc jest to dobrze udowodnione.

Baruch-Mordo, S., JS Evans, J. Severson, JD Naugle, J. Kiesecker, J. Maestas i MJ Falkowski (2013) Ratowanie cietrzewia z drzew: Proaktywne rozwiązanie w celu zmniejszenia kluczowego zagrożenia dla kandydata gatunek Ochrona biologiczna 167: 233-241

Poznanovic, AJ, MJ Falkowski, AL Maclean i JS Evans (2014) Ocena dokładności algorytmów wykrywania drzew w lasach jałowca. Inżynieria fotogrametryczna i teledetekcja 80 (5): 627–637

Istnieje kilka interesujących podejść, w ogólnym rozkładzie obiektów, z literatury stosowanej w matematyce w przestrzeni stanów wykorzystującej wielorozdzielcze procesy gaussowskie do dekompozycji cech obiektów w skali. Używam tego typu modeli do opisywania procesu wieloskalowego w modelach ekologicznych, ale można go dostosować do dekompozycji cech obiektu obrazowego. Zabawne, ale nieco ezoteryczne.

Gramacy, RB i HKH Lee (2008) Bayesian potraktowali modele procesów Gaussa za pomocą aplikacji do modelowania komputerowego. Journal of American Statistics Association, 103 (483): 1119–1130

Kim, HM, BK Mallick i CC Holmes (2005) Analiza niestacjonarnych danych przestrzennych przy użyciu cząstkowych procesów Gaussa. Journal of American Statistics Association, 100 (470): 653–668

Aby utworzyć DHM, odejmij DEM od DEM, można to zrobić w Esri Raster Calculator lub GDAL_CALC . To sprawi, że wszystkie twoje wzniesienia znajdą się na „równym polu gry”.

Składnia (Zamień pełne ścieżki dla DEM, DSM i DHM):

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

DHM będzie w większości wynosić 0 (lub wystarczająco blisko), co czynisz swoją wartością nodata. Za pomocą kalkulatora rastrowego lub GDAL_CALC można wyodrębnić wartości bardziej niż dowolną wartość na podstawie ilości hałasu zaobserwowanego w DHM. Ma to na celu zmniejszenie hałasu i podkreślenie tylko koron roślinności - w przypadku, gdy dwa „drzewa” sąsiadują ze sobą, powinno ono rozdzielić się na dwie wyraźne plamy.

Składnia (Zamień pełne ścieżki na binarne i DHM i zaobserwowaną wartość na wartość):

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

Teraz za pomocą GDAL_CALC lub Esri IsNull utwórz binarny raster, który można poligonizować za pomocą GDAL_Polygonize lub Esri Raster do Polygon .

Aby udoskonalić wielokąty, usuń zbyt małe wielokąty, a następnie porównaj je z pasmami RGB szukającymi podpisów, w Esri pomoże narzędzie Statystyka strefowa . Następnie możesz odrzucić wielokąty, które wyraźnie nie mają odpowiednich statystyk (na podstawie eksperymentów i twoich danych, nie mogę podać ci wartości).

To powinno doprowadzić cię do około 80% dokładności przy kreśleniu poszczególnych koron.

eCognition jest najlepszym oprogramowaniem do tego, zrobiłem to przy użyciu innego oprogramowania, ale eCognition jest lepszy. Oto odniesienie do literatury na ten temat:

Karlson, M., Reese, H., i Ostwald, M. (2014). Mapowanie korony drzew w zarządzanych lasach (Parklands) w półsuchej Afryce Zachodniej za pomocą obrazów WorldView-2 i analizy obrazów na podstawie obiektów geograficznych. Czujniki, 14 (12), 22643-22669.

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

Do tego:

Zagalikis, G., Cameron, AD, i Miller, DR (2005). Zastosowanie cyfrowej fotogrametrii i technik analizy obrazu w celu uzyskania charakterystyk drzew i drzewostanów. Kanadyjski dziennik badań leśnych, 35 (5), 1224-1237.

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

Miałem ten sam problem kilka lat temu. Mam rozwiązanie, które nie wymaga przefiltrowanych danych LAS ani innych danych pomocniczych. Jeśli masz dostęp do danych LiDAR i możesz generować DEMs / DSMs / DHMs (DEM poniżej, nie debatuję nad semantyką nazewnictwa modeli powierzchni) z różnych zwrotów, poniższy skrypt może być przydatny.

Skrypt arcpy pobiera 3 DEM i wyrzuca wielokąt leśny i pliki kształtów w kształcie drzewa. 3 DEM powinny mieć tę samą rozdzielczość przestrzenną (tj. 1 metr) i zakresy oraz reprezentować pierwsze powroty, ostatnie powroty i gołą ziemię. Miałem bardzo specyficzne parametry do ekstrakcji warzyw, ale parametry można zmienić w celu dostosowania do innych potrzeb. Jestem pewien, że proces ten można ulepszyć, ponieważ była to moja pierwsza poważna próba skryptowania w języku 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")

Podaję to jako odpowiedź ze względu na limit długości w komentarzu, brak nadziei na kredyty :). Bardzo szeroki pędzel, pod warunkiem, że masz DEM.

1. Wyodrębnij DEM dla pojedynczego wielokąta do dem.
2. Zdefiniuj ekstremum wysokości dema
3. Ustaw zCur + = - zStep. Krok do znalezienia przez iteracje wcześniej, np. Rozsądny spadek między wysokością „wierzchołka drzewa” a sąsiadami
4. Poniżej = Con (dem => zCur, int (1))
5. Zgrupuj regiony poniżej. Policz wystarczająco duże, czyli „drzewa”. Definicja wymagana tutaj przez kontrolę wizualną, wstępne badania?
6. Idź do kroku 3, jeśli zCur> zMin, krok 1 w przeciwnym razie.

Maksymalna liczba grup w procesie = liczba drzew wewnątrz pojedynczego wielokąta. Dodatkowe kryteria, np. Odległość między „drzewami” wewnątrz wielokątów, mogą pomóc ... Wygładzanie DEM przy użyciu jądra również jest opcją.