В этом примере показано, как извлечь лесные метрики, и отдельное дерево приписывает из воздушных данных о лидаре.
Лесное исследование и приложения все больше используют данные о лидаре, полученные от бортовых лазерных систем сканирования. Данные об облаке точек из лидара высокой плотности включают измерение не только лесные метрики, но также и атрибуты отдельных деревьев.
Этот пример использует данные об облаке точек из файла LAZ, полученного бортовой системой лидара, как введено. В этом примере вы сначала извлекаете лесные метрики путем классификации данных об облаке точек в землю и точки растительности, и затем извлекаете отдельные древовидные атрибуты путем сегментации точек растительности на отдельные деревья. Этот рисунок предоставляет обзор процесса.
Разархивируйте forestData.zip к временной директории и загрузке данные об облаке точек из файла LAZ, forestData.laz, в рабочую область MATLAB®. Данные получены из Открытого Набора данных Топографии [1]. Облако точек, в основном, содержит точки растительности и земля. Загрузите данные об облаке точек в рабочую область с помощью readPointCloud функция lasFileReader объект. Визуализируйте облако точек с помощью pcshow функция.
dataFolder = fullfile(tempdir,"forestData",filesep); dataFile = dataFolder + "forestData.laz"; % Check whether the folder and data file already exist or not folderExists = exist(dataFolder,'dir'); fileExists = exist(dataFile,'file'); % Create a new folder if it doesn't exist if ~folderExists mkdir(dataFolder); end % Extract aerial data file if it doesn't exist if ~fileExists unzip('forestData.zip',dataFolder); end % Read LAZ data from file lasReader = lasFileReader(dataFile); % Read point cloud along with corresponding scan angle information [ptCloud,pointAttributes] = readPointCloud(lasReader,"Attributes","ScanAngle"); % Visualize the input point cloud figure pcshow(ptCloud.Location) title("Input Point Cloud")
Основывайтесь сегментация является шагом предварительной обработки, чтобы изолировать данные о растительности для извлечения лесных метрик. Сегментируйте данные, загруженные из файла LAZ в наземные и неназемные точки с помощью segmentGroundSMRF функция.
% Segment Ground and extract non-ground and ground points groundPtsIdx = segmentGroundSMRF(ptCloud); nonGroundPtCloud = select(ptCloud,~groundPtsIdx); groundPtCloud = select(ptCloud,groundPtsIdx); % Visualize non-ground and ground points in magenta and green, respectively figure pcshowpair(nonGroundPtCloud,groundPtCloud) title("Segmented Non-Ground and Ground Points")
Используйте нормализацию вертикального изменения, чтобы устранить эффект ландшафта на ваших данных о растительности. Используйте точки с нормированным вертикальным изменением, как введено для вычислительных лесных метрик и древовидных атрибутов. Это шаги для нормализации вертикального изменения.
Устраните дублирующиеся точки вдоль x-и осей Y, если таковые имеются, при помощи groupsummary функция.
Создайте interpolant использование scatteredInterpolant объект, чтобы оценить землю в каждой точке в данных об облаке точек.
Нормируйте вертикальное изменение каждой точки путем вычитания интерполированного наземного вертикального изменения из исходного вертикального изменения.
groundPoints = groundPtCloud.Location; % Eliminate duplicate points along x- and y-axes [uniqueZ,uniqueXY] = groupsummary(groundPoints(:,3),groundPoints(:,1:2),@mean); uniqueXY = [uniqueXY{:}]; % Create interpolant and use it to estimate ground elevation at each point F = scatteredInterpolant(double(uniqueXY),double(uniqueZ),"natural"); estElevation = F(double(ptCloud.Location(:,1)),double(ptCloud.Location(:,2))); % Normalize elevation by ground normalizedPoints = ptCloud.Location; normalizedPoints(:,3) = normalizedPoints(:,3) - estElevation; % Visualize normalized points figure pcshow(normalizedPoints) title("Point Cloud with Normalized Elevation")
Извлеките лесные метрики из нормированных точек с помощью helperExtractForestMetrics функция помощника, присоединенная к этому примеру как вспомогательный файл. Функция помощника сначала делит облако точек на сетки на основе обеспеченного gridSize, и затем вычисляет лесные метрики. Функция помощника принимает что все точки с нормированной высотой ниже, чем cutoffHeight земля, и остающиеся точки являются растительностью. Вычислите эти лесные метрики.
Навес (CC) — Навес [2] является пропорцией леса, покрытого вертикальной проекцией древовидных корон. Вычислите его, когда отношение растительности возвращается относительно общего количества возвратов.
Часть разрыва (GF) — часть Разрыва [3] является вероятностью луча света, проходящего через навес, не сталкиваясь с листвой или другими элементами объекта. Вычислите его, когда отношение земли возвращается относительно общего количества возвратов.
Индекс [3] области Leaf area index (LAI) — Leaf является суммой односторонней листовой области на модуль земельного участка. Значение LAI вычисляется с помощью уравнения , где ang средний угол сканирования, GF часть разрыва и k коэффициент исчезновения, который тесно связан с распределением листового угла.
% Set grid size to 10 meters per pixel and cutOffHeight to 2 meters gridSize = 10; cutOffHeight = 2; leafAngDistribution = 0.5; % Extract forest metrics [canopyCover,gapFraction,leafAreaIndex] = helperExtractForestMetrics(normalizedPoints, ... pointAttributes.ScanAngle,gridSize,cutOffHeight,leafAngDistribution); % Visualize forest metrics hForestMetrics = figure; axCC = subplot(2,2,1,Parent=hForestMetrics); axCC.Position = [0.05 0.51 0.4 0.4]; imagesc(canopyCover,Parent=axCC) title(axCC,"Canopy Cover") axis off colormap(gray) axGF = subplot(2,2,2,Parent=hForestMetrics); axGF.Position = [0.55 0.51 0.4 0.4]; imagesc(gapFraction,'Parent',axGF) title(axGF,"Gap Fraction") axis off colormap(gray) axLAI = subplot(2,2,[3 4],Parent=hForestMetrics); axLAI.Position = [0.3 0.02 0.4 0.4]; imagesc(leafAreaIndex,Parent=axLAI) title(axLAI,"Leaf Area Index") axis off colormap(gray)
Модели высоты навеса (CHMs) являются растровыми представлениями высоты деревьев, созданий и других структур над землей топография. Используйте CHM в качестве входа для древовидного обнаружения и сегментации. Сгенерируйте CHM от своих нормированных значений вертикального изменения с помощью pc2dem функция.
% Set grid size to 0.5 meters per pixel gridRes = 0.5; % Generate CHM canopyModel = pc2dem(pointCloud(normalizedPoints),gridRes,CornerFillMethod="max"); % Clip invalid and negative CHM values to zero canopyModel(isnan(canopyModel) | canopyModel<0) = 0; % Perform gaussian smoothing to remove noise effects H = fspecial("gaussian",[5 5],1); canopyModel = imfilter(canopyModel,H,'replicate','same'); % Visualize CHM figure imagesc(canopyModel) title('Canopy Height Model') axis off colormap(gray)
Обнаружьте вершины деревьев с помощью helperDetectTreeTops функция помощника, присоединенная к этому примеру как вспомогательный файл. Функция помощника обнаруживает вершины деревьев путем нахождения локальных максимумов в переменных размерах окна [4] в CHM. Для обнаружения вершины дерева функция помощника рассматривает только вопросы с нормированной высотой, больше, чем minTreeHeight.
% Set minTreeHeight to 5 m minTreeHeight = 5; % Detect tree tops [treeTopRowId,treeTopColId] = helperDetectTreeTops(canopyModel,gridRes,minTreeHeight); % Visualize treetops figure imagesc(canopyModel) hold on plot(treeTopColId,treeTopRowId,"rx",MarkerSize=3) title("CHM with Detected Tree Tops") axis off colormap("gray")
Деревья индивидуума сегмента с помощью helperSegmentTrees функция помощника, присоединенная к этому примеру как вспомогательный файл. Функция помощника использует управляемую маркером сегментацию водораздела [5], чтобы сегментировать отдельные деревья. Во-первых, функция создает бинарное изображение маркера с местоположениями вершины дерева, обозначенными значением 1 . Затем функция фильтрует дополнение CHM наложением минимумов, чтобы удалить минимумы, которые не являются вершинами деревьев. Функция затем выполняет сегментацию водораздела на отфильтрованном дополнении CHM, чтобы сегментировать отдельные деревья. После сегментации визуализируйте отдельные древовидные сегменты.
% Segment individual trees label2D = helperSegmentTrees(canopyModel,treeTopRowId,treeTopColId,minTreeHeight); % Identify row and column id of each point in label2D and transfer labels % to each point rowId = ceil((ptCloud.Location(:,2) - ptCloud.YLimits(1))/gridRes) + 1; colId = ceil((ptCloud.Location(:,1) - ptCloud.XLimits(1))/gridRes) + 1; ind = sub2ind(size(label2D),rowId,colId); label3D = label2D(ind); % Extract valid labels and corresponding points validSegIds = label3D ~= 0; ptVeg = select(ptCloud,validSegIds); veglabel3D = label3D(validSegIds); % Assign color to each label numColors = max(veglabel3D); colorMap = randi([0 255],numColors,3)/255; labelColors = label2rgb(veglabel3D,colorMap,OutputFormat="triplets"); % Visualize tree segments figure pcshow(ptVeg.Location,labelColors) title("Individual Tree Segments") view(2)
Извлеките отдельные древовидные атрибуты с помощью helperExtractTreeMetrics функция помощника, присоединенная к этому примеру как вспомогательный файл. Во-первых, функция идентифицирует точки, принадлежащие отдельным деревьям от меток. Затем функция извлекает древовидные атрибуты, такие как древовидное местоположение вершины вдоль x-и осей Y, аппроксимированной древовидной высоты, древовидного диаметра короны и области. Функция помощника возвращает атрибуты как таблицу, где каждая строка представляет атрибуты отдельного дерева.
% Extract tree attributes treeMetrics = helperExtractTreeMetrics(normalizedPoints,label3D); % Display first 5 tree segments metrics disp(head(treeMetrics,5));
TreeId NumPoints TreeApexLocX TreeApexLocY TreeHeight CrownDiameter CrownArea
______ _________ ____________ ____________ __________ _____________ _________
1 388 2.6867e+05 4.1719e+06 29.509 7.5325 44.562
2 22 2.6867e+05 4.1719e+06 21.464 0.99236 0.77344
3 243 2.6867e+05 4.1719e+06 24.201 5.7424 25.898
4 101 2.6867e+05 4.1719e+06 21.927 3.4571 9.3867
5 54 2.6867e+05 4.1719e+06 19.515 3.0407 7.2617
[1] Томпсон, Обзор Лидара Области с Ручьем. Иллилуетт, Йосемитская долина, CA 2018. Национальный Центр Бортового Лазера, Сопоставляющего (NCALM). Распределенный OpenTopography. https://doi.org/10.5069/G96M351N. Полученный доступ: 2021-05-14
[2] Мама, Цинь, Янджун Су и Цинхуа Го. “Сравнение Оценок Навеса От Бортового LiDAR, Воздушного Формирования изображений и Спутниковых снимков”. Журнал IEEE Выбранных Тем в Прикладных наблюдениях Земли и Дистанционном зондировании 10, № 9 (сентябрь 2017): 4225–36. https://doi.org/10.1109/JSTARS.2017.2711482.
[3] Ричардсон, Джеффри Дж., Л. Моника Москаль и Су-Хюнг Ким. "Моделирование Подходов, чтобы Оценить Эффективный Листовой индекс области от Дискретной Антенны - Возвращает ЛИДАР". Сельскохозяйственный и Лесная Метеорология 149, № 6-7 (июнь 2009): 1152–60. https://doi.org/10.1016/j.agrformet.2009.02.007.
[4] Pitkänen, J., М. Мэлтэмо, Дж. Хииппэ и Кс. Ю. "Методы адаптации для Отдельного Древовидного Обнаружения на Бортовой основанной на лазере Модели Высоты Навеса". Международные Архивы Фотограмметрии, Дистанционного зондирования и Пространственной Информатики 36, № 8 (январь 2004): 187–91.
[5] Чен, Ци, Деннис Болдоччи, Пэн Гун и Маджи Келли. “Изолируя Отдельные Деревья в Лесистой местности Саванны Используя Маленькие Данные о Лидаре Места”. Photogrammetric Engineering & Remote Sensing 72, № 8 (1 августа 2006): 923–32. https://doi.org/10.14358/PERS.72.8.923.