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

13

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.

przykładowe zdjęcia lotnicze

Theodore Jones
źródło
Masz więcej danych niż link. Czy masz sklasyfikowane pliki Las, czy tylko raster DEM / DSM (który?)? To naprawdę nie jest łatwo to zrobić z fal tylko wizualnych z jakimkolwiek stopniem dokładności.
Michael Stimson,
Prawdopodobnie powinienem był podać więcej szczegółów na temat tego, co mam. Surowe zestawy danych, które posiadam to: 1.Pełne dane lasów i wygenerowany z nich riff tiff 2. Obrazy wizualne (jak pokazano przykładowy obraz, ale obejmujący znacznie szerszy obszar) 3. Ręczne bezpośrednie pomiary podzbioru drzew w strefa. Z nich wygenerowałem: 1. klasyfikacje gruntów / roślinności 2. rastry DEM / DSM
Theodore Jones
Czy masz dostęp do eCognition? Jeśli nie, do jakiego oprogramowania do przetwarzania obrazu lub języków programowania masz dostęp lub znasz?
Aaron
Nie mam kopii eCognition, ale sprawdzę, czy ma ją ktoś, kogo znam w moim laboratorium / na uniwersytecie, ponieważ wydaje się ona popularna w tego typu sprawach. Mam wiedzę w języku Python, C i Javie. Mam kopię Matlaba, ale jestem w tym prawie nieobcy. Mam dostęp do dowolnego oprogramowania z tej listy softwarelicense.arizona.edu/students oraz, oczywiście, ArcGIS.
Theodore Jones,
Trochę więcej szczegółów w moich komercyjnych aplikacjach. Niektóre z tej listy oprogramowania, które połączyłem, to Matlab, Mathematica, JMP i inne narzędzia statystyczne oraz narzędzia programistyczne, takie jak Visual Studio.
Theodore Jones,

Odpowiedzi:

10

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

Jeffrey Evans
źródło
+1 Szczególnie dla opcji 4; Ponieważ OP ma dane lidarowe, warto byłoby zastosować metodę falkową na modelu powierzchni czaszy. Chociaż, jak wiadomo, metoda falkowa nie jest tak naprawdę głównym nurtem (a może nigdy).
Aaron
W odejściu od jednego uniwersalnego ideału zacznę odnosić się do oprogramowania komercyjnego (np. ESRI, ERDAS) jako oprogramowania Big-box. Często najlepsze rozwiązanie lub jakiekolwiek inne nie jest dostępne w „oprogramowaniu Big-box”. Często trzeba szukać społeczności rozwojowych lub akademickich, aby znaleźć odpowiedzi na złożone przestrzenne problemy analityczne. To bardzo szybko sprowadza cię z głównego nurtu. Na szczęście społeczności te lubią się dzielić. Dlatego ważne jest, aby analityk nie polegał na rozwiązaniach obsługiwanych przyciskiem.
Jeffrey Evans,
2
Zgadzam się na temat BBS w przypadku złożonych problemów przestrzennych. Jednak wyodrębnianie jednego rodzaju roślinności w suchym środowisku - szczególnie jeśli masz dostęp do danych lidar - jest dość powszechne. W takim przypadku nie trzeba ponownie wymyślać koła, opracowując nowe podejście do prostej identyfikacji drzewa. Zastanawiam się, dlaczego nie zastosować wcześniej ustalonego podejścia opartego na naciśnięciu przycisku, zwłaszcza w pakiecie takim jak eCognition, który doskonale nadaje się do automatyzacji?
Aaron
1
Powinienem dodać, że eCognition ma zdolność do indywidualnego identyfikowania korony. Jako przykład możesz znaleźć przykładowy zestaw reguł w społeczności eCog, który stosuje podejście do uprawy nasion - wyszukaj „Próbka zestawu reguł oliwy z drzewa palmowego”. Integracja nowego algorytmu dopasowywania szablonów eCog i tego podejścia do uprawy nasion może potencjalnie być bardzo skuteczną metodą.
Aaron
1
Interesuje mnie kod python, o którym wspomniałeś dla metody Wavelet Smitha (2008). Czy jest dostępny gdziekolwiek?
Alpheus
3

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.

Michael Stimson
źródło
Dzięki. Zobaczę, czy uzyskam dobre wyniki dzięki tej metodzie.
Theodore Jones,
Będziesz musiał przeprowadzić pewne eksperymenty, aby uzyskać odpowiednie wartości, sugeruję obcinanie małych obszarów jako próbek, które wskazują (podobne do) najlepszych / najgorszych obszarów w twoich danych. Uzyskanie parametrów może zająć pół tuzina próbek, ale to wciąż musi być lepsze niż wykreślanie ich ręcznie.
Michael Stimson,
3

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

Giorgos Zagalikis
źródło
Czy mógłby Pan wyjaśnić, dlaczego eCognition jest lepszy? Tylko odpowiedzi linku stają się nieczynne, gdy link znika.
Aaron
1
eCognition to oprogramowanie do analizy obrazów oparte na obiektach, których nie ma od chwili obecnej. Zastosowałem podobne podejście. Zastosowanie cyfrowej fotogrametrii i technik analizy obrazu w celu uzyskania charakterystyk drzew i drzewostanów 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
Dzięki za referencje Giorgo. Sądzę, że te komentarze sprawdzą się jako odpowiedź na twoją odpowiedź.
Aaron
3

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")
Barbarossa
źródło
2

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ą.

FelixIP
źródło
Uważam, że masz na myśli DSM, a nie DEM ... Zazwyczaj drzewa, struktury i inne śmieci nie przekształcają się w DEM, ale są funkcją DSM (bez klas szumów). DSM - DEM = DHM (model wysokości). Wszystkie te dane można w rozsądny sposób wyodrębnić z danych LiDAR, nawet jeśli są one tylko klasyfikowane jako naziemne / niepoprawne, ale jeśli masz DEM, a nie LAS, jesteś w potoku bez wiosła, ponieważ funkcje, których szukasz, nie są tam !
Michael Stimson,
Tak, DHM, jak opisano, zrobi. Niewiele wiem o Lidar.
FelixIP,