4.3 Main loop
The main loop runs through the tiles of tree crown polygons, thus executing the working steps outlined in the introduction. For this markup-demo this loop is not required as only a reduced project region is considered. Hence, the loop is de-activated and the working steps are executed only once.
<- 1 # the loop-index needs to be fixed as the loop is de-activated
i
#for (i in 1:length(list_tile_names)) {
# Read the attribute data (dbf) of the tree crown tile (shp-file)
<- read.dbf(list_tile_names[i])
df_crown_parameter str(df_crown_parameter) # display the data frame structure
## 'data.frame': 22437 obs. of 10 variables:
## $ ID : int 0 1 2 3 4 5 6 7 8 9 ...
## $ WOODTYPE : Factor w/ 1 level "BW": 1 1 1 1 1 1 1 1 1 1 ...
## $ EASTTOP : num 3444039 3444023 3444501 3444082 3444106 ...
## $ NORTHTOP : num 5374000 5374000 5374000 5374000 5374000 ...
## $ HEIGHTTOP: num 21.4 16.2 18.6 40.5 36.3 ...
## $ HEIGHTDTM: num 840 850 821 819 810 ...
## $ AREACROWN: num 25.4 40.9 20.8 69.9 66.9 ...
## $ RAD1ELLIP: num 3.45 3.02 2.45 4.58 5.11 ...
## $ RAD2ELLIP: num 2.56 4.41 2.78 4.82 4.57 ...
## $ DIRELLIP : num -0.948 0.883 1.422 0.84 0.124 ...
## - attr(*, "data_types")= chr [1:10] "N" "C" "N" "N" ...
<- df_crown_parameter[,c(1,3:7)] # eliminate unrequired fields
df_crown_parameter head(df_crown_parameter) # display the top lines of the data frame
## ID EASTTOP NORTHTOP HEIGHTTOP HEIGHTDTM AREACROWN
## 1 0 3444039 5374000 21.350010 840.24 25.41760
## 2 1 3444023 5374000 16.170004 849.59 40.89117
## 3 2 3444501 5374000 18.610003 820.95 20.80000
## 4 3 3444082 5374000 40.490001 819.28 69.92512
## 5 4 3444106 5374000 36.310002 809.81 66.88000
## 6 5 3444653 5374000 4.650006 868.15 3.68000
# only for demo purpose we plot the crown polygons
<- st_read(
sf_crown_parameter paste0(tools::file_path_sans_ext(list_tile_names[i]),".shp"))
## Reading layer `3444_5373_DSM_04m_br_G200_M50000000' from data source
## `Y:\Baumartenklassifikation\R_script\2021_publication\data_markup\data_input\lidar\tree_crown_polygon_shp\3444_5373_DSM_04m_br_G200_M50000000.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 22437 features and 10 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 3443993 ymin: 5372993 xmax: 3445006 ymax: 5374013
## CRS: NA
# as shp-file lacks a coordinate reference system (crs), it needs to be set
st_crs(sf_crown_parameter) = st_crs(sf_nlp_1000m_boundary_poly)
# subset polygons to markup_boundary
<- sf_crown_parameter[sf_markup_boundary_poly,]
sf_crown_parameter ggplot() +
geom_sf(data = st_geometry(sf_crown_parameter)) +
geom_sf(data = st_geometry(sf_markup_boundary_poly),
aes(colour = "markup_boundary", alpha=0.9))
4.3.0.1 Apply the crown property parameters to filter treetops
# Exclude crowns that are too small and too big or high or are outliers
<- subset(df_crown_parameter,
df_crown_parameter >=AREACROWN_min & AREACROWN<=AREACROWN_max &
AREACROWN>=HEIGHTTOP_min & HEIGHTTOP<=HEIGHTTOP_max )
HEIGHTTOP
# apply the same filter to the polygons and plot the subset
<- subset(sf_crown_parameter,
sf_crown_parameter >=AREACROWN_min & AREACROWN<=AREACROWN_max &
AREACROWN>=HEIGHTTOP_min & HEIGHTTOP<=HEIGHTTOP_max )
HEIGHTTOPggplot() +
geom_sf(data = st_geometry(sf_crown_parameter)) +
geom_sf(data = st_geometry(sf_markup_boundary_poly),
aes(colour = "markup_boundary", alpha=0.9))
#### Apply regions to filter treetops
# convert crown peak point coordinates to simple feature points
<- st_as_sf(x = df_crown_parameter,
sf_crown_points coords = c("EASTTOP", "NORTHTOP"),
crs = st_crs(sf_nlp_1000m_boundary_poly))
str(sf_crown_points) # display structure of simple feature points
## Classes 'sf' and 'data.frame': 14894 obs. of 5 variables:
## $ ID : int 0 1 2 3 4 6 7 8 9 11 ...
## $ HEIGHTTOP: num 21.4 16.2 18.6 40.5 36.3 ...
## $ HEIGHTDTM: num 840 850 821 819 810 ...
## $ AREACROWN: num 25.4 40.9 20.8 69.9 66.9 ...
## $ geometry :sfc_POINT of length 14894; first list element: 'XY' num 3444039 5374000
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
## ..- attr(*, "names")= chr [1:4] "ID" "HEIGHTTOP" "HEIGHTDTM" "AREACROWN"
# clip points by nlp 1000m buffer boundary polygon
<- sf_crown_points[sf_nlp_1000m_boundary_poly,]
sf_crown_points
# Generate field Feature_ID to join to raster-extract to lidar_metrics
# as Feature_ID starts at 0 we need to subtract 1 from nrow
$Feature_ID <-
sf_crown_pointsas.integer(seq(0, nrow(sf_crown_points)-1), 1)
# change the field order (only for cosmetic purpose)
<- sf_crown_points[,c(1,6,2,3,4,5)]
sf_crown_points
# limit the project area (for this markup only)
if(sw_test_run == TRUE) {
<- sf_crown_points[sf_markup_boundary_poly,]
sf_crown_points
ggplot() +
ggtitle("Treetop points inside the markdown_boundary")+
geom_sf(data = st_geometry(sf_crown_parameter)) +
geom_sf(data = st_geometry(sf_crown_points))+
geom_sf(data = st_geometry(sf_markup_boundary_poly),
aes(colour = "markup_boundary", alpha=0.9))
}
str(sf_crown_points)
## Classes 'sf' and 'data.frame': 361 obs. of 6 variables:
## $ ID : int 13300 13301 13302 13318 13319 13331 13340 13341 13349 13363 ...
## $ Feature_ID: int 8922 8923 8924 8936 8937 8948 8953 8954 8962 8971 ...
## $ HEIGHTTOP : num 21 17.2 18 20.7 17.7 ...
## $ HEIGHTDTM : num 947 950 950 949 951 ...
## $ AREACROWN : num 39.2 29.4 29.6 42.7 23 ...
## $ geometry :sfc_POINT of length 361; first list element: 'XY' num 3444782 5373450
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
## ..- attr(*, "names")= chr [1:5] "ID" "Feature_ID" "HEIGHTTOP" "HEIGHTDTM" ...
4.3.0.2 Generate point buffer to extract raster information
<- st_buffer(sf_crown_points,
sf_crown_point_buffers dist = buffer_radius)
ggplot() +
ggtitle("Treetop buffers") +
geom_sf(data = st_geometry(sf_crown_parameter)) +
geom_sf(data = st_geometry(sf_crown_point_buffers))+
geom_sf(data = st_geometry(sf_markup_boundary_poly),
aes(colour = "markup_boundary", alpha=0.9) )
# Exclude those buffers not falling completely within the markdown_boundary
# as the rasters are clipped to that boundary,
# hence part of the buffer would extract no information
<-
sf_crown_point_buffers st_within(sf_crown_point_buffers,
sf_crown_point_buffers[%>% lengths > 0,]
sf_markup_boundary_poly)
ggplot() +
ggtitle("Treetop buffers falling completely within the markdown_boundary") +
geom_sf(data = st_geometry(sf_crown_parameter)) +
geom_sf(data = st_geometry(sf_crown_point_buffers))+
geom_sf(data = st_geometry(sf_markup_boundary_poly),
aes(colour = "markup_boundary", alpha=0.9) )
str(sf_crown_point_buffers)
## Classes 'sf' and 'data.frame': 342 obs. of 6 variables:
## $ ID : int 13383 13384 13401 13415 13429 13436 13445 13446 13466 13467 ...
## $ Feature_ID: int 8988 8989 9002 9011 9022 9027 9034 9035 9051 9052 ...
## $ HEIGHTTOP : num 18.5 17.5 19.6 18.1 18.6 ...
## $ HEIGHTDTM : num 944 951 951 951 951 ...
## $ AREACROWN : num 17.3 29 33.1 35.2 29.9 ...
## $ geometry :sfc_POLYGON of length 342; first list element: List of 1
## ..$ : num [1:121, 1:2] 3444762 3444762 3444762 3444762 3444762 ...
## ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
## ..- attr(*, "names")= chr [1:5] "ID" "Feature_ID" "HEIGHTTOP" "HEIGHTDTM" ...
# as the number of buffers had been reduced,
# the sf_crown_points require adjustment
<- sf_crown_points[sf_crown_points$ID %in%
sf_crown_points $ID, ] sf_crown_point_buffers
4.3.0.3 Calculate ALS (Lidar)-metrics for treetop buffers
See chapter ‘Working Steps’ for details on the ALS data used here. Collect the ASCII-files holding the targets of the treetop buffers. File names indicate the Feature-ID of treetop buffer polygon.
# compile list of ASCII-files holding each the ASL-targets for one tree
# hence the list holds all the tree file names for one tile
<- list.files(
list_lidar_targets_file_names path = list_lidar_targets_dir_names[i],
pattern = "\\.asc$", full.names = TRUE)
head(list_lidar_targets_file_names)
## [1] "data_markup/data_input/lidar/points_altitude-normalized_asc/3444_5373_ACge10_ACse400_HTge15_HTle70/Rawdata_Polygon0.rwb.asc"
## [2] "data_markup/data_input/lidar/points_altitude-normalized_asc/3444_5373_ACge10_ACse400_HTge15_HTle70/Rawdata_Polygon1.rwb.asc"
## [3] "data_markup/data_input/lidar/points_altitude-normalized_asc/3444_5373_ACge10_ACse400_HTge15_HTle70/Rawdata_Polygon10.rwb.asc"
## [4] "data_markup/data_input/lidar/points_altitude-normalized_asc/3444_5373_ACge10_ACse400_HTge15_HTle70/Rawdata_Polygon100.rwb.asc"
## [5] "data_markup/data_input/lidar/points_altitude-normalized_asc/3444_5373_ACge10_ACse400_HTge15_HTle70/Rawdata_Polygon1000.rwb.asc"
## [6] "data_markup/data_input/lidar/points_altitude-normalized_asc/3444_5373_ACge10_ACse400_HTge15_HTle70/Rawdata_Polygon10000.rwb.asc"
# extract the tile_name from the data path
<- list_lidar_targets_dir_names[i] %>%
tile_name strsplit("/") %>%
unlist() %>%
tail(n=1)
# as calculation of lidar metrics is time consuming,
# during test runs read the results from a RDS file
if (sw_calculate_lidar_metrics == FALSE) {
<-
df_lidar_metrics_tile readRDS(paste0(path_output_data, "TileLidarMetrics_", tile_name, ".rds"))
<- 1 # set index value
j
else {
}
# loop through all crown polygon Lidar targets files of this tile
for (j in 1:length(list_lidar_targets_file_names)) {
# extract tree polygon Feature_ID from df_lidar_targets targets file name
<- list_lidar_targets_file_names[j] %>%
feature_id strsplit("/") %>% # split path string into list according to "/"
unlist() %>% # unlist list into vector
tail(n=1) %>% # file name is last element of vector
parse_number() # extract tree polygon number from file name
# only if a record for this feature_id
# exists in the sf_crown_points
if (feature_id %in% sf_crown_points$Feature_ID) {
<- read.table(list_lidar_targets_file_names[j],
df_lidar_targets header = TRUE, sep = ' ', dec = '.')
# add Feature_ID as additional column to table df_lidar_targets
$Feature_ID <- feature_id
df_lidar_targets
### dataframe to collect the aggregated targets parameters
# one record per tree
# calculate mean values for intensity, number of targets and
# target number for all points in df_lidar_targets
# Feature_ID is required for a JOIN to sf_crown_points
<- summarise_at(df_lidar_targets,
df_lidar_metrics_tree .vars = c("X.Intensity.", "X.TargetNumber.",
"X.NumTargets.", "Feature_ID"),
list(mean))
# sum number of targets
$Sum_NumTargets_all <-
df_lidar_metrics_treeas.numeric(count(df_lidar_targets))
# adjust column names
names(df_lidar_metrics_tree) <- c(
"Mean_Intensity_all",
"Mean_TargetNum_all",
"Mean_NumTargets_all",
"Feature_ID",
"Sum_Targets_all")
## subset targets of the top 1.5m of the tree
<- filter(df_lidar_targets,
df_lidar_targets_top150cm $X.PositionH. >
df_lidar_targetsmax(df_lidar_targets$X.PositionH.) - 1.5))
(
# for the targets of the top 1.5m of the tree
# calculate mean values for intensity, mean number of targets and
# the mean target number for all targets in df_lidar_targets
# Feature_ID is required for a JOIN to df_lidar_metrics_tree
<- summarise_at(df_lidar_targets_top150cm,
df_lidar_metrics_tree_top150cm .vars = c("X.Intensity.", "X.TargetNumber.",
"X.NumTargets.","Feature_ID"),
list(mean))
# count top 1.5m number of targets of the tree
$Sum_Targets_top150cm <-
df_lidar_metrics_tree_top150cmas.numeric(count(df_lidar_targets_top150cm))
# get ratio of number of top 1.5m targets to number of total targets
<-
df_lidar_metrics_tree_top150cm ::mutate(df_lidar_metrics_tree_top150cm,
dplyrRatio_Targets_top150cm_to_all =
/
Sum_Targets_top150cm $Sum_Targets_all)
df_lidar_metrics_tree
# adjust column names
names(df_lidar_metrics_tree_top150cm) <- c(
"Mean_Intensity_top150cm",
"Mean_TargetNum_top150cm",
"Mean_NumTargets_top150cm",
"Feature_ID",
"Sum_Targets_top150cm",
"Ratio_Targets_top150cm_to_all")
# JOIN statistics of top 1.5m to main data frame
<- dplyr::full_join(df_lidar_metrics_tree,
df_lidar_metrics_tree
df_lidar_metrics_tree_top150cm, by = "Feature_ID") %>%
::select(., c("Feature_ID",
dplyr"Ratio_Targets_top150cm_to_all"))
# collect results into a list that stores Lidar metrics results
# for all trees of one tile
<- df_lidar_metrics_tree
list_lidar_metrics_tree[[j]]
# if (feature_id %in% sf_crown_points$Feature_ID)
} # for (j in 1:length(list_lidar_targets_file_names))
}
# convert the collect list_lidar_metrics_tree to a data frame
<- as.data.frame(bind_rows(list_lidar_metrics_tree, .id = NULL))
df_lidar_metrics_tile str(df_lidar_metrics_tile)
if (sw_save_results_as_rds== TRUE) {
saveRDS(df_lidar_metrics_tile, file = paste0(path_output_data,
"TileLidarMetrics_", tile_name, ".rds"))
# if
}
# else i.e. (sw_calculate_lidar_metrics == TRUE)
}
# list to collect metrics of all tiles
<- df_lidar_metrics_tile list_lidar_metrics_all_tiles[[i]]
4.3.0.4 Extract raster values using the treetop buffers
# as raster extraction metrics is time consuming,
# during test runs read the results from a RDS file
if (sw_extract_rasters == FALSE) {
<-
sf_crown_points readRDS(paste0(path_output_data, "TileRasterExtracts_", tile_name, ".rds"))
else {
}
$Ratio_2015_Red.Mean <-
sf_crown_pointsas.numeric(raster::extract(r_ratio_2015_red,
fun = mean))
sf_crown_point_buffers, $Ratio_2015_Green.Mean <-
sf_crown_pointsas.numeric(raster::extract(r_ratio_2015_gre,
fun = mean))
sf_crown_point_buffers, $Ratio_2015_Blue.Mean <-
sf_crown_pointsas.numeric(raster::extract(r_ratio_2015_blu,
fun = mean))
sf_crown_point_buffers, $Ratio_2015_NIR.Mean <-
sf_crown_pointsas.numeric(raster::extract(r_ratio_2015_nir,
fun = mean))
sf_crown_point_buffers,
$Deri_1_NIR_2015.Mean <-
sf_crown_pointsas.numeric(raster::extract(r_deri1_nir_2015,
fun = mean))
sf_crown_point_buffers,
$PCA_1_2015.Mean <-
sf_crown_pointsas.numeric(raster::extract( r_pca_1_2015,
fun = mean))
sf_crown_point_buffers,
$NDVI_2014.Mean <-
sf_crown_pointsas.numeric(raster::extract(r_ndvi_2014,
fun = mean))
sf_crown_point_buffers, $NDVI_2015.Mean <-
sf_crown_pointsas.numeric(raster::extract(r_ndvi_2015,
fun = mean))
sf_crown_point_buffers, $NDVI_2014.Min <-
sf_crown_pointsas.numeric(raster::extract(r_ndvi_2014,
fun = min))
sf_crown_point_buffers,
<-
sf_crown_points ::mutate(sf_crown_points,
dplyrNDVI_Diff_2014_2015.Mean = NDVI_2014.Mean - NDVI_2015.Mean )
# collect the results of each tile in a list
<- sf_crown_points
list_raster_extract_all_tiles[[i]]
if (sw_save_results_as_rds== TRUE) {
saveRDS(sf_crown_points, file = paste0(path_output_data,
"TileRasterExtracts_", tile_name, ".rds"))
# if
} # else i.e. (sw_extract_rasters == TRUE)
}
str(sf_crown_points)
## Classes 'sf' and 'data.frame': 342 obs. of 17 variables:
## $ ID : int 13383 13384 13401 13415 13429 13436 13445 13446 13466 13467 ...
## $ Feature_ID : int 8988 8989 9002 9011 9022 9027 9034 9035 9051 9052 ...
## $ HEIGHTTOP : num 18.5 17.5 19.6 18.1 18.6 ...
## $ HEIGHTDTM : num 944 951 951 951 951 ...
## $ AREACROWN : num 17.3 29 33.1 35.2 29.9 ...
## $ geometry :sfc_POINT of length 342; first list element: 'XY' num 3444761 5373447
## $ Ratio_2015_Red.Mean : num 68.8 73.7 91.4 90.6 82.2 ...
## $ Ratio_2015_Green.Mean : num 57.9 65.8 85 80.3 68.9 ...
## $ Ratio_2015_Blue.Mean : num 29.8 21.2 54 35.2 32.3 ...
## $ Ratio_2015_NIR.Mean : num 251 247 202 216 234 ...
## $ Deri_1_NIR_2015.Mean : num 90 90 90 90 90 ...
## $ PCA_1_2015.Mean : num 127 123 125 124 126 ...
## $ PCA_2_2015.Mean : num 144 142 140 138 141 ...
## $ NDVI_2014.Mean : num 190 189 159 169 175 ...
## $ NDVI_2015.Mean : num 193 188 160 165 177 ...
## $ NDVI_2014.Min : num 172 169 148 153 148 132 145 132 169 169 ...
## $ NDVI_Diff_2014_2015.Mean: num -3.222 0.731 -0.972 4.185 -1.343 ...
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:16] "ID" "Feature_ID" "HEIGHTTOP" "HEIGHTDTM" ...
4.3.0.5 JOIN lidar-metrics to raster-extracted parameters
# the join fields Feature_ID must be of the same type, i.e. integer
$Feature_ID <-
df_lidar_metrics_tileas.integer(df_lidar_metrics_tile$Feature_ID)
# JOIN lidar-metrics to raster-extracted features
<- as.data.frame(inner_join(sf_crown_points,
df_tile_attributes
df_lidar_metrics_tile,by = "Feature_ID"))
if (sw_save_results_as_rds == TRUE) { # save results of one tile to rds-file
saveRDS(df_tile_attributes,
file = paste0(path_output_data, "TileJoinedParameters_", tile_name, ".rds"))
}
# collect results into a list that stores all results for all trees of one tile
<- df_tile_attributes
list_attributes_all_tiles[[i]]
# --- close the main for-loop loop
# } # for
# ----------------