Extraindo áreas da copa das árvores a partir de dados de sensoriamento remoto (imagens visuais e LiDAR)

12

Estou procurando um método para processar uma imagem de sensoriamento remoto e extrair as áreas de copa das árvores individuais da imagem.

Eu tenho imagens de área de comprimento de onda visual e dados da área. A localização em questão é uma área deserta, portanto a cobertura das árvores não é tão densa quanto uma área florestal. A resolução das imagens aéreas é de 0,5 pés por 0,5 pés. A resolução do lidar é de aproximadamente 1 x 1 pé. Os dados visuais e o lidar vêm de um conjunto de dados de Pima County, Arizona. Uma amostra do tipo de imagem aérea que tenho está no final deste post.

Esta pergunta Detecção de árvore única no ArcMap? parece ser o mesmo problema, mas não parece haver uma boa resposta lá.

Posso obter uma classificação razoável dos tipos de vegetação (e informações sobre a porcentagem geral de cobertura) na área usando a classificação Iso Cluster no Arcmap, mas isso fornece pouca informação sobre árvores individuais. O mais próximo que eu tenho do que eu quero é os resultados de passar a saída da classificação isocluster pelo recurso Raster to Polygon no Arcmap. O problema é que esse método se funde próximo a árvores em um único polígono.

Edit: Eu provavelmente deveria ter incluído mais alguns detalhes sobre o que tenho. Os conjuntos de dados brutos que tenho são:

  • Dados completos, e um raster tiff gerado a partir deles.
  • Imagens visuais (como a imagem de exemplo mostrada, mas cobrindo uma área muito maior)
  • Medições diretas manuais de um subconjunto de árvores na área.

A partir destes eu gerei:

  1. As classificações de solo / vegetação.
  2. Os rasters DEM / DSM.

amostra de imagens aéreas

Theodore Jones
fonte
Você tem mais dados que o link. Você tem os arquivos classificados ou apenas a varredura DEM / DSM (qual?)? É realmente não é fácil de fazer isso com comprimentos de onda apenas visuais com algum grau de precisão.
Michael Stimson
Eu provavelmente deveria ter incluído mais alguns detalhes sobre o que tenho. Os conjuntos de dados brutos que eu tenho são: 1.Dados completos e uma raster gerada a partir deles 2. Imagens visuais (como a imagem de exemplo mostrada, mas cobrindo uma área muito maior) 3. medições diretas manuais de um subconjunto de árvores em a área. A partir disso, eu gerei: 1. as classificações de solo / vegetação 2. os rasters DEM / DSM
Theodore Jones
Você tem acesso ao eCognition? Caso contrário, a que software de processamento de imagem ou linguagens de programação você tem acesso ou conhece?
Aaron
Não tenho uma cópia do eCognition, mas vou verificar se alguém que conheço no meu laboratório / universidade tem, porque parece popular para esse tipo de coisa. Eu tenho conhecimento em Python, C e Java. Eu tenho uma cópia do Matlab, mas sou praticamente um noob nisso. Eu tenho acesso a qualquer um dos softwares desta lista softwarelicense.arizona.edu/students , além, é claro, do ArcGIS.
Theodore Jones
Um pouco mais detalhadamente nas aplicações comerciais que tenho. Alguns dos itens nessa lista de softwares que vinculei são Matlab, Mathematica, JMP e outras ferramentas estatísticas, e ferramentas de desenvolvimento de software, como o Visual Studio.
Theodore Jones

Respostas:

10

Existe um corpo considerável de literatura sobre detecção individual de coroas em dados espectrais e de lidar. Métodos, talvez comece com:

Falkowski, MJ, AMS Smith, PE Gessler, AT Hudak, LA Vierling e JS Evans. (2008). A influência do dossel de floresta de coníferas cobre a precisão de dois algoritmos de medição de árvores individuais usando dados do lidar. Jornal Canadense de Sensoriamento Remoto 34 (2): 338-350.

Smith AMS, EK Strand, CM Steele, DB Hann, SR Garrity, MJ Falkowski, JS Evans (2008) Produção de mapas de estrutura espacial da vegetação por análise por objeto da invasão de zimbro em fotografias aéreas multitemporais. Canadian Journal Sensoriamento Remoto 34 (2): 268-285

Se você está interessado no método Wavelet (Smith et al., 2008), eu o codifiquei em Python, mas é muito lento. Se você tem experiência com o Matlab, é aqui que ele é implementado no modo de produção. Temos dois documentos em que identificamos ~ 6 milhões de acres de invasão de zimbro no leste do Oregon usando o método wavelet com imagens NAIP RGB-NIR, portanto, está bem comprovado.

Baruch-Mordo, S., JS Evans, J. Severson, JD Naugle, J. Kiesecker, J. Maestas e MJ Falkowski (2013). Conservação Biológica das espécies 167: 233-241

Poznanovic, AJ, MJ Falkowski, AL Maclean e JS Evans (2014) Uma Avaliação de Precisão de Algoritmos de Detecção de Árvores em Juniper Woodlands. Engenharia fotogramétrica e sensoriamento remoto 80 (5): 627–637

Existem algumas abordagens interessantes, na decomposição geral de objetos, a partir da literatura aplicada ao espaço de estados da matemática, usando processos gaussianos de várias soluções para decompor as características dos objetos através da escala. Eu uso esses tipos de modelos para descrever processos em várias escalas em modelos ecológicos, mas ele pode ser adaptado para decompor as características dos objetos de imagem. Divertido, mas um pouco esotérico.

Gramacy, RB e HKH Lee (2008) Bayesian utilizam modelos de processos gaussianos com uma aplicação para modelagem por computador. Jornal da Associação Estatística Americana, 103 (483): 1119-1130

Kim, HM, BK Mallick e CC Holmes (2005) Analisando dados espaciais não-estacionários usando processos Gaussianos fragmentados. Jornal da Associação Estatística Americana, 100 (470): 653–668

Jeffrey Evans
fonte
+1 Especialmente para a opção 4; Como o OP possui dados do Lidar, vale a pena executar o método wavelet em um modelo de superfície de cobertura. Embora, como você sabe, o método wavelet ainda não seja realmente popular (ou talvez nunca).
Aaron
Em uma ode ao ideal de tamanho único, vou começar a me referir a software comercial (por exemplo, ESRI, ERDAS) como software Big-box. Geralmente, a melhor solução, ou qualquer outra, não está disponível no "software Big-box". Freqüentemente, é preciso procurar nas comunidades acadêmicas ou de desenvolvimento respostas para problemas analíticos espaciais complexos. Isso leva você para fora do mainstream com muita pressa. Felizmente, essas comunidades gostam de compartilhar. É também por isso que é importante para um analista não confiar em soluções de botão de pressão.
Jeffrey Evans
2
Costumo concordar sobre o BBS para problemas espaciais complexos. No entanto, a extração de um único tipo de vegetação em um ambiente árido - especialmente se você tiver acesso aos dados do lidar - é bastante comum. Nesse caso, não há necessidade de reinventar a roda, desenvolvendo uma nova abordagem para a identificação simples de árvores. Meu pensamento é por que não usar uma abordagem pré-estabelecida por botões, especialmente em um pacote como o eCognition, que é altamente adequado para automação?
Aaron
1
Devo acrescentar que o eCognition tem capacidade para identificação individual da coroa. Como exemplo, você pode encontrar um conjunto de regras de amostra na comunidade eCog que usa uma abordagem de cultivo de sementes - procure "Conjunto de regras de amostra de delineação de palmeiras de óleo". A integração do novo algoritmo de correspondência de modelos do eCog e essa abordagem de cultivo de sementes pode ser um método muito poderoso.
Aaron
1
Estou interessado no código python que você mencionou para o método Wavelet de Smith (2008). Está disponível em algum lugar?
Alpheus 29/04
4

O eCognition é o melhor software para isso, eu fiz isso usando outro software, mas o eCognition é melhor. Aqui está a referência à literatura sobre o assunto:

Karlson, M., Reese, H. e Ostwald, M. (2014). Mapeamento de copas de árvores em florestas gerenciadas (áreas de parques) da África Ocidental semi-árida usando imagens do WorldView-2 e análise de imagens com base em objetos geográficos. Sensores, 14 (12), 22643-22669.

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

Além disso:

Zagalikis, G., Cameron, AD e Miller, DR (2005). Aplicação de fotogrametria digital e técnicas de análise de imagem para derivar características de árvores e suportes. Jornal canadense de pesquisa florestal, 35 (5), 1224-1237.

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

Giorgos Zagalikis
fonte
Você poderia expandir por que o eCognition é melhor? As respostas apenas ao link tendem a se tornar deficientes quando o link desaparece.
Aaron
1
O eCognition é um software de análise de imagem baseado em objetos que não o são desde agora. Eu usei uma abordagem semelhante. A aplicação de técnicas de fotogrametria digital e de análise de imagem para derivar características de árvores e suportes G Zagalikis, AD Cameron, DR Miller Jornal Canadense de Pesquisa Florestal, 2005, 35 (5): 1224-1237, 10.1139 / x05-030 nrcresearchpress.com/doi /abs/10.1139/x05-030#.VJmMb14gAA
Giorgos Zagalikis
1
Obrigado pela referência Giorgos. Acho que esses comentários funcionariam bem como uma edição da sua resposta.
Aaron
3

Para criar um DHM subtrair o DEM do DEM, isso pode ser feito na Esri Raster Calculator ou GDAL_CALC . Isso colocará todas as suas elevações em um 'campo de jogo nivelado'.

Sintaxe (substituir caminhos completos para DEM, DSM e DHM):

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

O DHM será principalmente 0 (ou próximo o suficiente), o que você faz com o valor do seu nodata. Com a Raster Calculator ou GDAL_CALC, você pode extrair valores mais que um valor arbitrário com base na quantidade de ruído observada no DHM. O objetivo disso é reduzir o ruído e destacar apenas as copas da vegetação - no caso em que duas 'árvores' estão adjacentes, isso deve se dividir em duas bolhas distintas.

Sintaxe (substitua os caminhos completos por Binário e DHM e o valor observado por Valor):

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

Agora, com GDAL_CALC ou Esri IsNull, crie uma varredura binária, que pode ser poligonalizada com GDAL_Polygonize ou Esri Raster to Polygon .

Para refinar os polígonos, remover polígonos excessivamente pequenos e compará-los às bandas RGB que procuram assinaturas. No Esri, a ferramenta Estatísticas Zonais ajudará. Então você pode descartar os polígonos que claramente não possuem as estatísticas corretas (com base na experimentação e nos seus dados, não posso fornecer os valores).

Isso deve levar você a cerca de 80% de precisão na plotagem de coroas individuais.

Michael Stimson
fonte
Obrigado. Vou ver se consigo bons resultados com esse método.
Theodore Jones
Você precisará fazer algumas experiências para obter os valores apropriados. Sugiro recortar pequenas áreas como amostras indicativas (semelhantes às) das melhores / piores áreas dos seus dados. Pode levar meia dúzia de execuções de amostra para obter seus parâmetros, mas isso ainda precisa ser melhor do que plotá-los manualmente.
Michael Stimson
3

Eu tive o mesmo problema alguns anos atrás. Eu tenho uma solução que não requer dados LAS filtrados ou outros dados auxiliares. Se você tiver acesso aos dados LiDAR e puder gerar DEMs / DSMs / DHMs (DEM a seguir, não discutirei a semântica da nomenclatura do modelo de superfície) de retornos diferentes, o script a seguir pode ser útil.

O script arcpy ingere 3 DEMs e cospe um polígono de floresta e arquivos de forma de ponto de árvore. Os 3 DEMs devem ter a mesma resolução espacial (ou seja, 1 metro) e extensões, e representam primeiros retornos, últimos retornos e terra nua. Eu tinha parâmetros muito específicos para a extração de vegetais, mas os parâmetros podem ser alterados para atender a outras necessidades. Estou certo de que o processo pode ser aprimorado, pois essa foi minha primeira tentativa séria de script em 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
fonte
2

Estou postando isso como resposta devido ao limite de tamanho no comentário, sem esperanças de créditos :). Pincel muito amplo, desde que você tenha DEM.

  1. Extraia DEM para polígono individual para dem.
  2. Definir extremos de elevação do dem
  3. Defina zCur + = - zStep. Etapa a ser encontrada pelas iterações de antemão, por exemplo, queda razoável entre a elevação da célula da árvore e os vizinhos
  4. Abaixo = Con (dem => zCur, int (1))
  5. Regiões do grupo Abaixo. Conte grande o suficiente, que são 'árvores'. Definição exigida aqui por inspeção visual, pesquisa preliminar?
  6. Vá para a etapa 3 se zCur> zMin, etapa 1 caso contrário.

Número máximo de grupos no processo = contagem de árvores no polígono individual. Critérios adicionais, por exemplo, a distância entre 'árvores' dentro dos polígonos podem ajudar ... A suavização do DEM usando o kernel também é uma opção.

FelixIP
fonte
Acredito que você esteja se referindo a um DSM e não a um DEM ... Normalmente, árvores, estruturas e outros itens indesejados não fazem dele um DEM, mas são apresentados no DSM (menos classes de ruído). DSM - DEM = DHM (modelo de altura). Tudo isso pode ser extraído razoavelmente a partir dos dados do LiDAR, mesmo que sejam apenas classificados como terra / não terra, mas se você tiver o DEM e não o LAS, estará no riacho sem remo, porque os recursos que você procura não são ai !
Michael Stimson
Sim, DHM como você descreveu fará. Eu sei pouco sobre Lidar.
precisa saber é o seguinte