Извлечение областей кроны деревьев из данных дистанционного зондирования (визуальные изображения и LiDAR)

13

Я ищу метод для обработки изображения дистанционного зондирования и извлечения областей кроны отдельных деревьев из изображения.

У меня есть как визуальные снимки с использованием длин волн, так и лидарные данные из этой области. Место, о котором идет речь, является пустынным районом, поэтому древесный покров не такой плотный, как лесной массив. Разрешение аэрофотоснимков составляет 0,5 фута на 0,5 фута. Разрешение лидара составляет приблизительно 1 x 1 фут. И визуальные данные, и лидар поступают из набора данных округа Пима, штат Аризона. Образец типа аэрофотоснимков, который у меня есть, находится в конце этого поста.

Этот вопрос Обнаружение единого дерева в ArcMap? кажется, та же самая проблема, но там, кажется, нет хорошего ответа.

Я могу получить разумную классификацию типов растительности (и информацию об общем процентном покрытии) в области, используя классификацию Iso Cluster в Arcmap, но это дает мало информации об отдельных деревьях. Ближе всего к тому, что я хочу, это результаты прохождения выходных данных классификации изокластера через растр в Polygon в Arcmap. Проблема в том, что этот метод объединяет деревья рядом в один многоугольник.

Изменить: я, вероятно, должен был бы включить некоторые более подробные сведения о том, что у меня есть. У меня есть наборы необработанных данных:

  • Полные данные LAS и сгенерированный из него TIFF растр.
  • Визуальные образы (например, показанный образец изображения, но охватывающий гораздо более широкую область)
  • Ручные прямые измерения подмножества деревьев в этом районе.

Из них я сгенерировал:

  1. Классификация почв / растительности.
  2. Растры DEM / DSM.

Образец аэрофотоснимков

Теодор Джонс
источник
У вас есть больше данных, чем ссылка. У вас есть секретные файлы las или только растр DEM / DSM (какой?)? Это действительно нелегко сделать с помощью только визуальных длин волн с любой степенью точности.
Майкл Стимсон
Я, наверное, должен был бы включить некоторые подробности о том, что у меня есть. У меня есть необработанные наборы данных: 1.Полные данные LAS и сгенерированный из них TIFF растр. 2. Визуальные образы (как показано на примере изображения, но охватывающего гораздо более широкую область). 3. Ручные прямые измерения подмножества деревьев в площадь. Из них я получил: 1. Классификации грунта / растительности 2. Растры DEM / DSM
Теодор Джонс
У вас есть доступ к eCognition? Если нет, то к какому программному обеспечению для обработки изображений или языкам программирования вы имеете доступ или знаете?
Аарон
У меня нет копии eCognition, но я проверю, есть ли у кого-нибудь из моих знакомых в моей лаборатории / университете, потому что она кажется популярной для такого рода вещей. Я хорошо разбираюсь в Python, C и Java. У меня есть копия Matlab, но я в значительной степени нуб. У меня есть доступ к любому программному обеспечению из этого списка softwarelicense.arizona.edu/students , плюс, конечно, ArcGIS.
Теодор Джонс
Чуть подробнее в коммерческих приложениях у меня есть. В этом списке программ, которые я связал, есть такие, как Matlab, Mathematica, JMP и другие инструменты статистики, а также инструменты разработки программного обеспечения, такие как Visual Studio.
Теодор Джонс

Ответы:

10

По спектральным и лидарным данным имеется значительное количество литературы по обнаружению индивидуальной короны. Мудрые методы, возможно, начнем с:

Фальковский, MJ, AMS Smith, PE Gessler, AT Hudak, LA Vierling и JS Evans. (2008). Влияние покрытия полога хвойного леса на точность двух отдельных алгоритмов измерения деревьев с использованием лидарных данных. Канадский журнал дистанционного зондирования 34 (2): 338-350.

Smith AMS, EK Strand, CM Steele, DB Hann, SR Garrity, MJ Falkowski, JS Evans (2008). Создание карт пространственной структуры растительности путем анализа каждого объекта вторжения можжевельника на многовременных аэрофотоснимках. Канадский журнал дистанционного зондирования 34 (2): 268-285

Если вас интересует метод Wavelet (Smith et al., 2008), я написал его на Python, но он очень медленный. Если у вас есть опыт работы с Matlab, то здесь он реализован в производственном режиме. У нас есть две работы, в которых мы идентифицировали ~ 6 миллионов акров вторжения можжевельника в восточном Орегоне, используя вейвлет-метод с изображениями NAIP RGB-NIR, поэтому это хорошо доказано.

Baruch-Mordo, S., JS Evans, J. Severson, JD Naugle, J. Kiesecker, J. Maestas и MJ Falkowski (2013) Спасение шалфея от деревьев: упреждающее решение для снижения ключевой угрозы для кандидата Биологическая консервация видов 167: 233-241

Познанович, AJ, MJ Falkowski, AL Maclean и JS Evans (2014) Оценка точности алгоритмов обнаружения деревьев в лесах можжевельника. Фотограмметрическая инженерия и дистанционное зондирование 80 (5): 627–637

Есть несколько интересных подходов, в общем случае декомпозиции объектов, из прикладной математической литературы по пространству состояний, использующей многоразрешающие гауссовские процессы для декомпозиции характеристик объектов в масштабе. Я использую эти типы моделей для описания многомасштабного процесса в экологических моделях, но он может быть адаптирован для разложения характеристик объекта изображения. Весело, но немного эзотерично.

Gramacy, RB и HKH Lee (2008). Модели гауссовского процесса с байесовским анализом с применением к компьютерному моделированию. Журнал Американской статистической ассоциации, 103 (483): 1119–1130

Kim, HM, BK Mallick, и CC Holmes (2005). Анализ нестационарных пространственных данных с использованием кусочно-гауссовских процессов. Журнал Американской статистической ассоциации, 100 (470): 653–668

Джеффри Эванс
источник
+1 Особенно для варианта 4; Поскольку у ОП есть лидарные данные, стоило бы запустить вейвлет-метод на модели поверхности купола. Хотя, как вы знаете, вейвлет-метод еще не очень распространен (или, возможно, когда-либо).
Аарон
В оду идеалу «один размер подходит всем» я начну ссылаться на коммерческое программное обеспечение (например, ESRI, ERDAS) как программное обеспечение Big-box. Зачастую самое лучшее решение или его вообще нет в «Big-box software». Часто приходится искать ответы на сложные пространственно-аналитические задачи в отношении образовательных или академических сообществ. Это очень быстро выводит вас из мейнстрима. К счастью, эти сообщества любят делиться. Это также, почему для аналитика важно не полагаться на кнопочные решения.
Джеффри Эванс
2
Я склонен согласиться с BBS для сложных пространственных задач. Однако выделение одного типа растительности в засушливой среде - особенно если у вас есть доступ к лидарным данным - является довольно распространенным явлением. В этом случае нет необходимости заново изобретать колесо, разрабатывая новый подход к простой идентификации деревьев. Я думаю, почему бы не использовать заранее установленный, кнопочный подход, особенно в пакете, таком как eCognition, который очень подходит для автоматизации?
Аарон
1
Я должен добавить, что у eCognition есть возможность индивидуальной идентификации короны. Например, вы можете найти примерный набор правил в сообществе eCog, в котором используется подход к семеноводству - поиск «Пример набора правил разведения пальмового масла». Интеграция нового алгоритма сопоставления шаблонов eCog и этого семенного подхода потенциально может быть очень мощным методом.
Аарон
1
Меня интересует код Python, который вы упоминаете для метода вейвлета Смита (2008). Это доступно где-нибудь?
Алфеус
3

Чтобы создать DHM вычесть DEM из DEM, это можно сделать в Esri Raster Calculator или GDAL_CALC . Это поместит все ваши возвышения на «ровное игровое поле».

Синтаксис (замена полных путей для DEM, DSM и DHM):

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

Значение DHM будет в основном 0 (или достаточно близко), что вы задаете в качестве значения нодаты. С помощью Raster Calculator или GDAL_CALC вы можете извлекать значения, превышающие произвольное значение, на основе количества шума, которое вы наблюдаете в DHM. Цель этого состоит в том, чтобы уменьшить шум и выделить только кроны растительности - в случае, когда два «дерева» соседствуют, это должно быть разделено на две отдельные капли.

Синтаксис (замените полные пути для двоичного и DHM и наблюдаемое значение для значения):

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

Теперь с помощью GDAL_CALC или Esri IsNull создайте двоичный растр, который можно полигонизировать с помощью GDAL_Polygonize или Esri Raster to Polygon .

Чтобы уточнить полигоны, удалите чрезмерно маленькие полигоны и затем сравните их с полосами RGB, которые ищут подписи, в Esri поможет инструмент Zonal Statistics . Затем вы можете отказаться от полигонов, которые явно не имеют правильной статистики (основываясь на экспериментах и ​​ваших данных, я не могу дать вам значения).

Это должно дать вам около 80% точности при построении отдельных коронок.

Майкл Стимсон
источник
Благодарю. Я посмотрю, получу ли я хорошие результаты от этого метода.
Теодор Джонс
Вам нужно будет поэкспериментировать, чтобы получить подходящие значения, я предлагаю вырезать небольшие области в качестве образцов, которые указывают (аналогично) лучшие / худшие области в ваших данных. Для получения ваших параметров может потребоваться полдюжины пробных прогонов, но это все же лучше, чем составлять их вручную.
Майкл Стимсон
3

eCognition - лучшее программное обеспечение для этого, я сделал это, используя другое программное обеспечение, но eCognition лучше. Вот ссылка на литературу на эту тему:

Karlson, M., Reese, H. & Ostwald, M. (2014). Картографирование крон деревьев в управляемых лесных массивах (парковых насаждениях) полузасушливой Западной Африки с использованием изображений WorldView-2 и анализа изображений на основе географических объектов. Датчики, 14 (12), 22643-22669.

например, http://www.mdpi.com/1424-8220/14/12/22643

Дополнительно:

Загаликис Г., Камерон А.Д. и Миллер Д.Р. (2005). Применение методов цифровой фотограмметрии и анализа изображений для получения характеристик деревьев и древостоев. Канадский журнал исследований леса, 35 (5), 1224-1237.

например, http://www.nrcresearchpress.com/doi/abs/10.1139/x05-030#.VJmMb14gAA

Гиоргос Загаликис
источник
Не могли бы вы рассказать, почему eCognition лучше? Ответы только на ссылки, как правило, перестают существовать, когда ссылка исчезает.
Аарон
1
eCognition - это объектно-ориентированное программное обеспечение для анализа изображений. Я использовал подобный подход. Применение методов цифровой фотограмметрии и анализа изображений для определения характеристик деревьев и древостоев Г. Загликис, А. Д. Кэмерон, Д. Р. Миллер, Канадский журнал исследований леса, 2005, 35 (5): 1224-1237, 10.1139 / x05-030 nrcresearchpress.com/doi /abs/10.1139/x05-030#.VJmMb14gAA
Гиоргос Загаликис
1
Спасибо за ссылку Гиоргос. Я думаю, что эти комментарии будут работать как редактирование вашего ответа.
Аарон
3

У меня была такая же проблема пару лет назад. У меня есть решение, которое не требует отфильтрованных данных LAS или других вспомогательных данных. Если у вас есть доступ к данным LiDAR, и вы можете генерировать DEM / DSM / DHM (далее DEM, я не обсуждаю семантику номенклатуры модели поверхности) из разных возвратов, может пригодиться следующий скрипт.

Скрипт arcpy принимает 3 матрицы высот и выплевывает лесной полигон и шейп-файлы деревьев. 3 матрицы высот должны иметь одинаковое пространственное разрешение (т.е. 1 метр) и экстенты, и представлять первые возвращения, последние возвращения и голую землю. У меня были очень специфические параметры для извлечения овощей, но параметры могут быть изменены в соответствии с другими потребностями. Я уверен, что этот процесс можно улучшить, так как это была моя первая серьезная попытка написания скриптов на 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
источник
2

Я отправляю это как ответ из-за ограничения длины в комментарии, никаких надежд на кредиты :). Очень широкая кисть, при условии, что у вас есть DEM.

  1. Извлечь DEM для отдельного полигона в дем.
  2. Определить крайние значения высоты демона
  3. Установите zCur + = - zStep. Шаг должен быть найден заранее итерациями, например, разумное падение между высотой 'верхней ячейки дерева' и соседями
  4. Ниже = Con (dem => zCur, int (1))
  5. Группа регионов ниже. Считай достаточно большой, это «деревья». Определение требуется здесь при визуальном осмотре, предварительном исследовании?
  6. Перейдите к шагу 3, если zCur> zMin, в противном случае к шагу 1.

Максимальное количество групп в процессе = количество деревьев внутри отдельного многоугольника. Дополнительные критерии, например, расстояние между «деревьями» внутри многоугольников, могут помочь ... Сглаживание матрицы высот с использованием ядра также вариант.

FelixIP
источник
Я полагаю, что вы имеете в виду DSM, а не DEM ... Обычно деревья, структуры и другие нежелательные объекты не превращаются в DEM, а используются в DSM (за исключением классов шума). DSM - DEM = DHM (модель высоты). Все они могут быть разумно извлечены из данных LiDAR, даже если они классифицируются только как наземные / наземные, но если у вас есть DEM, а не LAS, вы идете вверх по течению без весла, потому что функции, которые вы ищете, не являются там !
Майкл Стимсон
Да, DHM, как вы описали, будет делать. Я мало знаю о Лидаре.
FelixIP