主要用到的是R包：ForestGapR

1. 讀取lidar資料

lidR::readLAS

2. 生成CHM

lidR::grid_canopy

ps. 此處可有其他選擇，生成CHM講究也很多，生成raster

3. 識別林窗

ForestGapR::getForestGaps(chm_layer, threshold=10, size=c(10,10^3))

這裡根據地區調整林窗識別的條件：高度閾值和麵積檢測範圍

4.統計林窗

ForestGapR::GapStats(gap_layer, chm_layer)

5.計算GSFD

GapSizeFDist(gaps_stats,method,...)

6.出圖

程式碼如下：

 1 library(ForestGapR)
 2 library(raster)
 3 setwd('D:\\research\\......')
 4 #1.讀取lidar
 5 XXX_all <- lidR::readLAS('D:\\research\\ ......, threshold = 10, size = c(1,10000) )
 6 gaps_XXX %>% plot()
 7 dev.copy(png, 'GSFD XXX林冠間隙大小頻率分佈.png');dev.off()
 8 plot(gaps_XXX2 , col ="forestgreen
 
")
 9 
10 #2.生成chm
11 chm_XXX <- lidR::grid_canopy(las = XXX_all, res = 1, algorithm = pitfree(thresholds = c(0, 10, 20), max_edge = c(0, 1.5)))
12 plot(chm_XXX, col=col)
13 
14 #3.識別林窗，引數：10m，10~1000m2
15 gaps_XXX <- ForestGapR::getForestGaps(chm_layer = chm_XXX, threshold = 10, size = c(10,1000) )
16 gaps_XXX %>% plot()
 
17 dev.copy(png, 'GSFD XXX林冠間隙大小頻率分佈.png');dev.off()
18 
19 #4. 統計林窗
20 gaps_XXX_stats <- ForestGapR::GapStats(gap_layer = gaps_XXX, chm_layer = chm_XXX)
21 gaps_XXX_stats %>% DT::datatable()
22 
23 #5. GSFD
24 GapSizeFDist(
25   gaps_stats = gaps_XXX_stats,
26   method = 'Hanel_2017',
27   col = 'forestgreen',
28   pch = 16,
29   cex = 1,
30   axes = FALSE,
31   ylab = '間隙頻率',
32   xlab = as.expression(bquote('間隙大小' ~ (m^2)))
33 )
34 axis(1);axis(2)
35 grid(4,4)
36 dev.copy(png, 'GSFD間隙大小的頻率分佈.png');dev.off()
37 
38 #6.繪圖
39 par(mfrow=c(1,1))
40 #library(viridis)
41 plot(chm_XXX, main='XXX林冠間隙', col=viridis(10))
42 plot(gaps_XXX,add=TRUE, col='red',legend=FALSE)