Environmental Data Sherwood

Author

Josh Shaw

Load Packages

library(sf)
Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(dssd)
library(mapview)
library(ggplot2)
library(lidR)

Attaching package: 'lidR'
The following object is masked from 'package:sf':

    st_concave_hull
library(terra)
terra 1.8.15

Attaching package: 'terra'
The following objects are masked from 'package:lidR':

    area, crs, crs<-, is.empty, watershed
library(gstat)
library(ForestGapR)
library(viridis)
Loading required package: viridisLite
library(raster)
Loading required package: sp

Attaching package: 'sp'
The following object is masked from 'package:lidR':

    wkt

Attaching package: 'raster'
The following objects are masked from 'package:lidR':

    projection, projection<-
library(tidyverse)
Warning: package 'tidyverse' was built under R version 4.4.3
Warning: package 'tidyr' was built under R version 4.4.3
Warning: package 'readr' was built under R version 4.4.3
Warning: package 'dplyr' was built under R version 4.4.3
Warning: package 'forcats' was built under R version 4.4.3
Warning: package 'lubridate' was built under R version 4.4.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.4     ✔ tibble    3.2.1
✔ purrr     1.1.0     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract(), terra::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ dplyr::select()  masks raster::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(entropy)
Warning: package 'entropy' was built under R version 4.4.3

Attaching package: 'entropy'

The following object is masked from 'package:lidR':

    entropy
library(e1071)

Attaching package: 'e1071'

The following object is masked from 'package:raster':

    interpolate

The following object is masked from 'package:terra':

    interpolate
library(purrr)

Importing Data

Shapefiles

# Import Site Boundary
sssi <- read_sf(dsn = "./shapefiles", layer = "SSSI_boundary_revised")

# Import site boundary +200m for landscape scale measurements
boundary <- read_sf(dsn = "./shapefiles", layer = "boundary")

# Import Survey Plots
points <- read_sf(dsn = "./shapefiles", layer = "random_points")
transects <- read_sf(dsn = "./shapefiles", layer = "survey_plots")

# Assign BNG CRS
transects <- st_transform(transects, 27700)

# Remove Inaccessible points
points <- points %>%
  filter(!Number %in% c(18, 21, 25, 26, 27, 28))
transects <- transects %>%
  filter(!Number %in% c(18, 21, 25, 26, 27, 28))

# Check Alignment
ggplot() +
  geom_sf(data = transects[transects$Number == 1, ])+
  geom_sf(data = points[points$Number == 1, ])

# Check CRS match of all
crs(transects)
[1] "PROJCRS[\"OSGB36 / British National Grid\",\n    BASEGEOGCRS[\"OSGB36\",\n        DATUM[\"Ordnance Survey of Great Britain 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.\"],\n        BBOX[49.75,-9.01,61.01,2.01]],\n    ID[\"EPSG\",27700]]"
crs(points)
[1] "PROJCRS[\"OSGB36 / British National Grid\",\n    BASEGEOGCRS[\"OSGB36\",\n        DATUM[\"Ordnance Survey of Great Britain 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.\"],\n        BBOX[49.75,-9.01,61.01,2.01]],\n    ID[\"EPSG\",27700]]"
crs(boundary)
[1] "PROJCRS[\"OSGB36 / British National Grid\",\n    BASEGEOGCRS[\"OSGB36\",\n        DATUM[\"Ordnance Survey of Great Britain 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.\"],\n        BBOX[49.75,-9.01,61.01,2.01]],\n    ID[\"EPSG\",27700]]"
crs(sssi)
[1] "PROJCRS[\"OSGB36 / British National Grid\",\n    BASEGEOGCRS[\"OSGB36\",\n        DATUM[\"Ordnance Survey of Great Britain 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.\"],\n        BBOX[49.75,-9.01,61.01,2.01]],\n    ID[\"EPSG\",27700]]"
# All CRS match

Lidar

# Import LAS file
las <- readLAS("./lidar files/site.las")

# Normalise height of Las
las <- normalize_height(las, knnidw())
Inverse distance weighting: [==------------------------------------------------] 5% (4 threads)
Inverse distance weighting: [===-----------------------------------------------] 6% (4 threads)
Inverse distance weighting: [===-----------------------------------------------] 7% (4 threads)
Inverse distance weighting: [====----------------------------------------------] 8% (4 threads)
Inverse distance weighting: [====----------------------------------------------] 9% (4 threads)
Inverse distance weighting: [=====---------------------------------------------] 10% (4 threads)
Inverse distance weighting: [=====---------------------------------------------] 11% (4 threads)
Inverse distance weighting: [======--------------------------------------------] 12% (4 threads)
Inverse distance weighting: [======--------------------------------------------] 13% (4 threads)
Inverse distance weighting: [=======-------------------------------------------] 14% (4 threads)
Inverse distance weighting: [=======-------------------------------------------] 15% (4 threads)
Inverse distance weighting: [========------------------------------------------] 16% (4 threads)
Inverse distance weighting: [========------------------------------------------] 17% (4 threads)
Inverse distance weighting: [=========-----------------------------------------] 18% (4 threads)
Inverse distance weighting: [=========-----------------------------------------] 19% (4 threads)
Inverse distance weighting: [==========----------------------------------------] 20% (4 threads)
Inverse distance weighting: [==========----------------------------------------] 21% (4 threads)
Inverse distance weighting: [===========---------------------------------------] 22% (4 threads)
Inverse distance weighting: [===========---------------------------------------] 23% (4 threads)
Inverse distance weighting: [============--------------------------------------] 24% (4 threads)
Inverse distance weighting: [============--------------------------------------] 25% (4 threads)
Inverse distance weighting: [=============-------------------------------------] 26% (4 threads)
Inverse distance weighting: [=============-------------------------------------] 27% (4 threads)
Inverse distance weighting: [==============------------------------------------] 28% (4 threads)
Inverse distance weighting: [==============------------------------------------] 29% (4 threads)
Inverse distance weighting: [===============-----------------------------------] 30% (4 threads)
Inverse distance weighting: [===============-----------------------------------] 31% (4 threads)
Inverse distance weighting: [================----------------------------------] 32% (4 threads)
Inverse distance weighting: [================----------------------------------] 33% (4 threads)
Inverse distance weighting: [=================---------------------------------] 34% (4 threads)
Inverse distance weighting: [=================---------------------------------] 35% (4 threads)
Inverse distance weighting: [==================--------------------------------] 36% (4 threads)
Inverse distance weighting: [==================--------------------------------] 37% (4 threads)
Inverse distance weighting: [===================-------------------------------] 38% (4 threads)
Inverse distance weighting: [===================-------------------------------] 39% (4 threads)
Inverse distance weighting: [====================------------------------------] 40% (4 threads)
Inverse distance weighting: [====================------------------------------] 41% (4 threads)
Inverse distance weighting: [=====================-----------------------------] 42% (4 threads)
Inverse distance weighting: [=====================-----------------------------] 43% (4 threads)
Inverse distance weighting: [======================----------------------------] 44% (4 threads)
Inverse distance weighting: [======================----------------------------] 45% (4 threads)
Inverse distance weighting: [=======================---------------------------] 46% (4 threads)
Inverse distance weighting: [=======================---------------------------] 47% (4 threads)
Inverse distance weighting: [========================--------------------------] 48% (4 threads)
Inverse distance weighting: [========================--------------------------] 49% (4 threads)
Inverse distance weighting: [=========================-------------------------] 50% (4 threads)
Inverse distance weighting: [=========================-------------------------] 51% (4 threads)
Inverse distance weighting: [==========================------------------------] 52% (4 threads)
Inverse distance weighting: [==========================------------------------] 53% (4 threads)
Inverse distance weighting: [===========================-----------------------] 54% (4 threads)
Inverse distance weighting: [===========================-----------------------] 55% (4 threads)
Inverse distance weighting: [============================----------------------] 56% (4 threads)
Inverse distance weighting: [============================----------------------] 57% (4 threads)
Inverse distance weighting: [=============================---------------------] 58% (4 threads)
Inverse distance weighting: [=============================---------------------] 59% (4 threads)
Inverse distance weighting: [==============================--------------------] 60% (4 threads)
Inverse distance weighting: [==============================--------------------] 61% (4 threads)
Inverse distance weighting: [===============================-------------------] 62% (4 threads)
Inverse distance weighting: [===============================-------------------] 63% (4 threads)
Inverse distance weighting: [================================------------------] 64% (4 threads)
Inverse distance weighting: [================================------------------] 65% (4 threads)
Inverse distance weighting: [=================================-----------------] 66% (4 threads)
Inverse distance weighting: [=================================-----------------] 67% (4 threads)
Inverse distance weighting: [==================================----------------] 68% (4 threads)
Inverse distance weighting: [==================================----------------] 69% (4 threads)
Inverse distance weighting: [===================================---------------] 70% (4 threads)
Inverse distance weighting: [===================================---------------] 71% (4 threads)
Inverse distance weighting: [====================================--------------] 72% (4 threads)
Inverse distance weighting: [====================================--------------] 73% (4 threads)
Inverse distance weighting: [=====================================-------------] 74% (4 threads)
Inverse distance weighting: [=====================================-------------] 75% (4 threads)
Inverse distance weighting: [======================================------------] 76% (4 threads)
Inverse distance weighting: [======================================------------] 77% (4 threads)
Inverse distance weighting: [=======================================-----------] 78% (4 threads)
Inverse distance weighting: [=======================================-----------] 79% (4 threads)
Inverse distance weighting: [========================================----------] 80% (4 threads)
Inverse distance weighting: [========================================----------] 81% (4 threads)
Inverse distance weighting: [=========================================---------] 82% (4 threads)
Inverse distance weighting: [=========================================---------] 83% (4 threads)
Inverse distance weighting: [==========================================--------] 84% (4 threads)
Inverse distance weighting: [==========================================--------] 85% (4 threads)
Inverse distance weighting: [===========================================-------] 86% (4 threads)
Inverse distance weighting: [===========================================-------] 87% (4 threads)
Inverse distance weighting: [============================================------] 88% (4 threads)
Inverse distance weighting: [============================================------] 89% (4 threads)
Inverse distance weighting: [=============================================-----] 90% (4 threads)
Inverse distance weighting: [=============================================-----] 91% (4 threads)
Inverse distance weighting: [==============================================----] 92% (4 threads)
Inverse distance weighting: [==============================================----] 93% (4 threads)
Inverse distance weighting: [===============================================---] 94% (4 threads)
Inverse distance weighting: [===============================================---] 95% (4 threads)
Inverse distance weighting: [================================================--] 96% (4 threads)
Inverse distance weighting: [================================================--] 97% (4 threads)
Inverse distance weighting: [=================================================-] 98% (4 threads)
Inverse distance weighting: [=================================================-] 99% (4 threads)
Inverse distance weighting: [==================================================] 100% (4 threads)
# Check ground using hist, all should be 0
hist(filter_ground(las)$Z, breaks = seq(-0.6, 0.6, 0.01), main = "", xlab = "Elevation")

# Check las file
las_check(las)

 Checking the data
  - Checking coordinates... ✓
  - Checking coordinates type... ✓
  - Checking coordinates range... ✓
  - Checking coordinates quantization... ✓
  - Checking attributes type... ✓
  - Checking ReturnNumber validity... ✓
  - Checking NumberOfReturns validity... ✓
  - Checking ReturnNumber vs. NumberOfReturns... ✓
  - Checking RGB validity... ✓
  - Checking absence of NAs... ✓
  - Checking duplicated points... ✓
  - Checking degenerated ground points...
    ⚠ There were 7 degenerated ground points. Some X Y coordinates were repeated but with different Z coordinates
  - Checking attribute population...
    🛈 'ScanDirectionFlag' attribute is not populated
    🛈 'EdgeOfFlightline' attribute is not populated
  - Checking gpstime incoherances
    ✗ 38012 pulses (points with the same gpstime) have points with identical ReturnNumber
  - Checking flag attributes... ✓
  - Checking user data attribute...
    🛈 15227240 points have a non 0 UserData attribute. This probably has a meaning
 Checking the header
  - Checking header completeness... ✓
  - Checking scale factor validity... ✓
  - Checking point data format ID validity... ✓
  - Checking extra bytes attributes validity... ✓
  - Checking the bounding box validity... ✓
  - Checking coordinate reference system... ✓
 Checking header vs data adequacy
  - Checking attributes vs. point format... ✓
  - Checking header bbox vs. actual content... ✓
  - Checking header number of points vs. actual content... ✓
  - Checking header return number vs. actual content... ✓
 Checking coordinate reference system...
  - Checking if the CRS was understood by R... ✓
 Checking preprocessing already done 
  - Checking ground classification... yes
  - Checking normalization... yes
  - Checking negative outliers...
    ⚠ 21687 points below 0
  - Checking flightline classification... yes
 Checking compression
  - Checking attribute compression...
   -  ScanDirectionFlag is compressed
   -  EdgeOfFlightline is compressed
   -  Synthetic_flag is compressed
   -  Keypoint_flag is compressed
   -  Withheld_flag is compressed
   -  Overlap_flag is compressed
# Some points returning below 0 - check...
sub_zero <- filter_poi(las, Z < 0)
table(sub_zero@data$Classification) # Unclassified, low vegetation, medium vegetation and noise

    1     3     4     7 
 7682 12755  1205    45 
rm(sub_zero)

# Remove points below 0
las <- filter_poi(las, Z >= 0)

las_check(las)

 Checking the data
  - Checking coordinates... ✓
  - Checking coordinates type... ✓
  - Checking coordinates range... ✓
  - Checking coordinates quantization... ✓
  - Checking attributes type... ✓
  - Checking ReturnNumber validity... ✓
  - Checking NumberOfReturns validity... ✓
  - Checking ReturnNumber vs. NumberOfReturns... ✓
  - Checking RGB validity... ✓
  - Checking absence of NAs... ✓
  - Checking duplicated points... ✓
  - Checking degenerated ground points...
    ⚠ There were 7 degenerated ground points. Some X Y coordinates were repeated but with different Z coordinates
  - Checking attribute population...
    🛈 'ScanDirectionFlag' attribute is not populated
    🛈 'EdgeOfFlightline' attribute is not populated
  - Checking gpstime incoherances
    ✗ 38011 pulses (points with the same gpstime) have points with identical ReturnNumber
  - Checking flag attributes... ✓
  - Checking user data attribute...
    🛈 15213280 points have a non 0 UserData attribute. This probably has a meaning
 Checking the header
  - Checking header completeness... ✓
  - Checking scale factor validity... ✓
  - Checking point data format ID validity... ✓
  - Checking extra bytes attributes validity... ✓
  - Checking the bounding box validity... ✓
  - Checking coordinate reference system... ✓
 Checking header vs data adequacy
  - Checking attributes vs. point format... ✓
  - Checking header bbox vs. actual content... ✓
  - Checking header number of points vs. actual content... ✓
  - Checking header return number vs. actual content... ✓
 Checking coordinate reference system...
  - Checking if the CRS was understood by R... ✓
 Checking preprocessing already done 
  - Checking ground classification... yes
  - Checking normalization... yes
  - Checking negative outliers... ✓
  - Checking flightline classification... yes
 Checking compression
  - Checking attribute compression...
   -  ScanDirectionFlag is compressed
   -  EdgeOfFlightline is compressed
   -  Synthetic_flag is compressed
   -  Keypoint_flag is compressed
   -  Withheld_flag is compressed
   -  Overlap_flag is compressed
print(las)
class        : LAS (v1.4 format 6)
memory       : 1.7 Gb 
extent       : 458886.5, 463010.6, 366086.5, 368770.7 (xmin, xmax, ymin, ymax)
coord. ref.  : OSGB36 / British National Grid 
area         : 6.54 km²
points       : 23.99 million points
density      : 3.67 points/m²
density      : 2.11 pulses/m²
st_crs(las)
Coordinate Reference System:
  User input: PROJCS["OSGB36 / British National Grid",GEOGCS["OSGB36",DATUM["osgb_1936",SPHEROID["Airy 1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],TOWGS84[375,-111,431,0,0,0,0],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9102"]],AXIS["Geodetic longitude",EAST],AXIS["Geodetic latitude",NORTH],AUTHORITY["EPSG","4277"]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["latitude_of_origin",49],UNIT["Meter",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","27700"]] 
  wkt:
BOUNDCRS[
    SOURCECRS[
        PROJCRS["OSGB36 / British National Grid",
            BASEGEOGCRS["OSGB36",
                DATUM["Ordnance Survey of Great Britain 1936",
                    ELLIPSOID["Airy 1830",6377563.396,299.3249646,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]]],
            CONVERSION["unnamed",
                METHOD["Transverse Mercator",
                    ID["EPSG",9807]],
                PARAMETER["False easting",400000,
                    LENGTHUNIT["Meter",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",-100000,
                    LENGTHUNIT["Meter",1],
                    ID["EPSG",8807]],
                PARAMETER["Longitude of natural origin",-2,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["Scale factor at natural origin",0.9996012717,
                    SCALEUNIT["unity",1],
                    ID["EPSG",8805]],
                PARAMETER["Latitude of natural origin",49,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8801]]],
            CS[Cartesian,2],
                AXIS["easting",east,
                    ORDER[1],
                    LENGTHUNIT["Meter",1]],
                AXIS["northing",north,
                    ORDER[2],
                    LENGTHUNIT["Meter",1]],
            ID["EPSG",27700]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from OSGB36 to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",375,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",-111,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",431,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]

Creating basefiles

CHM

# create raster for canopy height model
chm <- rasterize_canopy(las, res = 0.5, algorithm = dsmtin())
Delaunay rasterization[=-------------------------------------------------] 2% (4 threads)
Delaunay rasterization[=-------------------------------------------------] 3% (4 threads)
Delaunay rasterization[==------------------------------------------------] 4% (4 threads)
Delaunay rasterization[==------------------------------------------------] 5% (4 threads)
Delaunay rasterization[===-----------------------------------------------] 6% (4 threads)
Delaunay rasterization[===-----------------------------------------------] 7% (4 threads)
Delaunay rasterization[====----------------------------------------------] 8% (4 threads)
Delaunay rasterization[====----------------------------------------------] 9% (4 threads)
Delaunay rasterization[=====---------------------------------------------] 10% (4 threads)
Delaunay rasterization[=====---------------------------------------------] 11% (4 threads)
Delaunay rasterization[======--------------------------------------------] 12% (4 threads)
Delaunay rasterization[======--------------------------------------------] 13% (4 threads)
Delaunay rasterization[=======-------------------------------------------] 14% (4 threads)
Delaunay rasterization[=======-------------------------------------------] 15% (4 threads)
Delaunay rasterization[========------------------------------------------] 16% (4 threads)
Delaunay rasterization[========------------------------------------------] 17% (4 threads)
Delaunay rasterization[=========-----------------------------------------] 18% (4 threads)
Delaunay rasterization[=========-----------------------------------------] 19% (4 threads)
Delaunay rasterization[==========----------------------------------------] 20% (4 threads)
Delaunay rasterization[==========----------------------------------------] 21% (4 threads)
Delaunay rasterization[===========---------------------------------------] 22% (4 threads)
Delaunay rasterization[===========---------------------------------------] 23% (4 threads)
Delaunay rasterization[============--------------------------------------] 24% (4 threads)
Delaunay rasterization[============--------------------------------------] 25% (4 threads)
Delaunay rasterization[=============-------------------------------------] 26% (4 threads)
Delaunay rasterization[=============-------------------------------------] 27% (4 threads)
Delaunay rasterization[==============------------------------------------] 28% (4 threads)
Delaunay rasterization[==============------------------------------------] 29% (4 threads)
Delaunay rasterization[===============-----------------------------------] 30% (4 threads)
Delaunay rasterization[===============-----------------------------------] 31% (4 threads)
Delaunay rasterization[================----------------------------------] 32% (4 threads)
Delaunay rasterization[================----------------------------------] 33% (4 threads)
Delaunay rasterization[=================---------------------------------] 34% (4 threads)
Delaunay rasterization[=================---------------------------------] 35% (4 threads)
Delaunay rasterization[==================--------------------------------] 36% (4 threads)
Delaunay rasterization[==================--------------------------------] 37% (4 threads)
Delaunay rasterization[===================-------------------------------] 38% (4 threads)
Delaunay rasterization[===================-------------------------------] 39% (4 threads)
Delaunay rasterization[====================------------------------------] 40% (4 threads)
Delaunay rasterization[====================------------------------------] 41% (4 threads)
Delaunay rasterization[=====================-----------------------------] 42% (4 threads)
Delaunay rasterization[=====================-----------------------------] 43% (4 threads)
Delaunay rasterization[======================----------------------------] 44% (4 threads)
Delaunay rasterization[======================----------------------------] 45% (4 threads)
Delaunay rasterization[=======================---------------------------] 46% (4 threads)
Delaunay rasterization[=======================---------------------------] 47% (4 threads)
Delaunay rasterization[========================--------------------------] 48% (4 threads)
Delaunay rasterization[========================--------------------------] 49% (4 threads)
Delaunay rasterization[=========================-------------------------] 50% (4 threads)
Delaunay rasterization[=========================-------------------------] 51% (4 threads)
Delaunay rasterization[==========================------------------------] 52% (4 threads)
Delaunay rasterization[==========================------------------------] 53% (4 threads)
Delaunay rasterization[===========================-----------------------] 54% (4 threads)
Delaunay rasterization[===========================-----------------------] 55% (4 threads)
Delaunay rasterization[============================----------------------] 56% (4 threads)
Delaunay rasterization[============================----------------------] 57% (4 threads)
Delaunay rasterization[=============================---------------------] 58% (4 threads)
Delaunay rasterization[=============================---------------------] 59% (4 threads)
Delaunay rasterization[==============================--------------------] 60% (4 threads)
Delaunay rasterization[==============================--------------------] 61% (4 threads)
Delaunay rasterization[===============================-------------------] 62% (4 threads)
Delaunay rasterization[===============================-------------------] 63% (4 threads)
Delaunay rasterization[================================------------------] 64% (4 threads)
Delaunay rasterization[================================------------------] 65% (4 threads)
Delaunay rasterization[=================================-----------------] 66% (4 threads)
Delaunay rasterization[=================================-----------------] 67% (4 threads)
Delaunay rasterization[==================================----------------] 68% (4 threads)
Delaunay rasterization[==================================----------------] 69% (4 threads)
Delaunay rasterization[===================================---------------] 70% (4 threads)
Delaunay rasterization[===================================---------------] 71% (4 threads)
Delaunay rasterization[====================================--------------] 72% (4 threads)
Delaunay rasterization[====================================--------------] 73% (4 threads)
Delaunay rasterization[=====================================-------------] 74% (4 threads)
Delaunay rasterization[=====================================-------------] 75% (4 threads)
Delaunay rasterization[======================================------------] 76% (4 threads)
Delaunay rasterization[======================================------------] 77% (4 threads)
Delaunay rasterization[=======================================-----------] 78% (4 threads)
Delaunay rasterization[=======================================-----------] 79% (4 threads)
Delaunay rasterization[========================================----------] 80% (4 threads)
Delaunay rasterization[========================================----------] 81% (4 threads)
Delaunay rasterization[=========================================---------] 82% (4 threads)
Delaunay rasterization[=========================================---------] 83% (4 threads)
Delaunay rasterization[==========================================--------] 84% (4 threads)
Delaunay rasterization[==========================================--------] 85% (4 threads)
Delaunay rasterization[===========================================-------] 86% (4 threads)
Delaunay rasterization[===========================================-------] 87% (4 threads)
Delaunay rasterization[============================================------] 88% (4 threads)
Delaunay rasterization[============================================------] 89% (4 threads)
Delaunay rasterization[=============================================-----] 90% (4 threads)
Delaunay rasterization[=============================================-----] 91% (4 threads)
Delaunay rasterization[==============================================----] 92% (4 threads)
Delaunay rasterization[==============================================----] 93% (4 threads)
Delaunay rasterization[===============================================---] 94% (4 threads)
Delaunay rasterization[===============================================---] 95% (4 threads)
Delaunay rasterization[================================================--] 96% (4 threads)
Delaunay rasterization[================================================--] 97% (4 threads)
Delaunay rasterization[=================================================-] 98% (4 threads)
Delaunay rasterization[=================================================-] 99% (4 threads)
Delaunay rasterization[==================================================] 100% (4 threads)
col <- height.colors(30)
plot(chm, col = col)

# Project CHM to British National Grid
chm   <- project(chm, "EPSG:27700")

# Need to convert chm to data frame for later plotting with ggplot
chm_df <- as.data.frame(chm, xy = TRUE, na.rm = TRUE)

Buffers

# Create buffer zones
buffer_50 <- st_buffer(points, dist = 50)
buffer_200 <- st_buffer(points, dist = 200)

# Assign CRS
buffer_50 <- st_transform(buffer_50, 27700) # Set to BNG EPSG:27700
buffer_200 <- st_transform(buffer_200, 27700)

crs(buffer_50)
[1] "PROJCRS[\"OSGB36 / British National Grid\",\n    BASEGEOGCRS[\"OSGB36\",\n        DATUM[\"Ordnance Survey of Great Britain 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.\"],\n        BBOX[49.75,-9.01,61.01,2.01]],\n    ID[\"EPSG\",27700]]"
crs(buffer_200)
[1] "PROJCRS[\"OSGB36 / British National Grid\",\n    BASEGEOGCRS[\"OSGB36\",\n        DATUM[\"Ordnance Survey of Great Britain 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.\"],\n        BBOX[49.75,-9.01,61.01,2.01]],\n    ID[\"EPSG\",27700]]"
ggplot() +
  geom_sf(data = buffer_200, bg = "transparent", col = "blue") +
  geom_sf(data = buffer_50, bg = "transparent") +
  geom_sf(data = boundary, bg = "transparent", col = "red") +
  geom_sf(data = sssi, bg = "transparent", col = "forestgreen") +
  geom_sf(data = points, bg = "transparent", col = "black") +
  coord_sf()

# Check transects and buffers align
ggplot() +
  geom_sf(data = transects[transects$Number == 1, ], color = "green", fill = NA, size = 1) +
  geom_sf(data = points[points$Number == 1, ], color = "black", size = 2) +
  geom_sf(data = buffer_50[buffer_50$Number == 1, ], color = "black", size = 2, fill = NA) +
  coord_sf(datum = NA)

Extracting Environmental Data

Canopy Height Diversity

# Creating custom forest structural diversity function
fsd_metrics <- function(z) {
  z <- z[z > 0]  # remove ground/no returns
  if (length(z) == 0) return(NULL)
  hist_vals <- hist(z, breaks = seq(0, max(z) + 1, by = 1), plot = FALSE) # Create histogram of height values
  prob <- hist_vals$counts / sum(hist_vals$counts) # Creates probability dist for entropy
  list(
    zmean = mean(z),
    zmax = max(z),
    zmin = min(z),
    zentropy = entropy::entropy(prob, unit = "log2"),
    zsd = sd(z),
    zskew = e1071::skewness(z))}

Per transect

# Calculate metrics per transect

# Ensure CRS match
# Convert LAS CRS to an sf-style CRS object
las_crs <- st_crs(las@crs)
# Reprojecting transects to match LAS
transects_crs <- st_transform(transects, crs = las_crs)

metrics_by_transect <- vector("list", length = nrow(transects_crs))

for (i in seq_len(nrow(transects_crs))) {
  poly <- transects_crs[i, ]
  las_clip <- clip_roi(las, poly)
  if (!lidR::is.empty(las_clip)) 
    {metrics_by_transect[[i]] <- fsd_metrics(las_clip$Z)} 
  else {metrics_by_transect[[i]] <- NA }}

# Create a dataframe of the metrics
metrics_df <- do.call(rbind, metrics_by_transect)
metrics_df <- as.data.frame(metrics_df) # Ensures it is a dataframe
metrics_df$transect_id <- transects_crs$Number

# Joins the metrics to the transect shapefile
height_transects <- cbind(transects_crs, metrics_df)

height_transects <- height_transects %>%
  rename(transect = Number) %>%
  select(-transect_id)

Per 50m Buffer

buffer50_crs <- st_transform(buffer_50, crs = las_crs)

metrics_buffer50 <- vector("list", length = nrow(buffer50_crs))

for (i in seq_len(nrow(buffer50_crs))) {
  poly <- buffer50_crs[i, ]
  las_clip <- clip_roi(las, poly)
  if (!lidR::is.empty(las_clip)) 
    {metrics_buffer50[[i]] <- fsd_metrics(las_clip$Z)}
  else {metrics_buffer50[[i]] <- NA }}

# Create a dataframe of the metrics
metrics50_df <- do.call(rbind, metrics_buffer50)
metrics50_df <- as.data.frame(metrics50_df) # Ensures it is a dataframe
metrics50_df$transect_id <- buffer50_crs$Number

# View Metrics  
metrics50_df
      zmean     zmax    zmin zentropy      zsd       zskew transect_id
1  11.69396 22.65396 0.00292 4.178257 6.295246  -0.5846616          12
2  7.734654 22.61012 0.00125 3.998712 5.933527   0.1929425           1
3  9.925839 25.02056   7e-04 4.247489 5.636467  -0.1895035           7
4  8.739816 27.13721 0.00049 4.250572 7.498247   0.3883114          29
5  11.45561 21.66573 0.00133 4.199989 5.595754  -0.5731387           9
6  12.57806 22.22637 0.00118 4.052151 5.381456  -0.9565064          11
7  5.957389 25.62101 0.00053 3.756337 5.306312     1.03431          19
8  8.503838 19.27454 0.00164 3.466638 6.274653  -0.3111754           5
9  2.362655 22.25043 0.00011  2.54298 3.267851    2.104486          30
10  7.78207 22.63401 0.00287 3.869971 6.243566   0.1199948          13
11 13.69688  28.9755 0.00042 3.926588 9.909547  -0.3623127          22
12 6.454479 19.65465 0.00123 3.900635 4.781247   0.5131124           4
13 12.08129 24.42534   1e-05 4.264428 7.040076  -0.4690887          17
14 7.851172 22.47729 0.00082 3.907197 6.129368    0.135253          14
15  8.19836 24.70047   8e-05 4.048821 6.451534   0.1818351           3
16 12.22543 30.98228   9e-05 4.370529 8.418593 -0.09380323          24
17 9.137033 21.05166 0.00255 3.983734 4.734016  -0.4528941           6
18 10.81797 21.84281 0.00073 4.047789 5.732416  -0.6232267          16
19 9.253438 22.98122 0.00051 4.128115 5.472973  -0.1874253           8
20 13.35016 28.69416 0.00178 4.289988 9.077425    -0.34593          20
21 8.112863 18.79273 0.00176 3.899887 5.229448  -0.2222557           2
22  7.28093 21.40331 0.00052 3.799095 5.505507  0.07082679          15
23 5.211221 18.06695 0.00014 3.602767 5.042661   0.8682786          10
24 3.857767 12.01605 0.00074 2.844996 3.218059  0.05901197          23
# Joins the metrics to the buffer50 shapefile
height_buffer50 <- cbind(buffer50_crs, metrics50_df)

height_buffer50 <- height_buffer50 %>%
  rename(transect = Number) %>%
  select(-transect_id)

Per 200m Buffer

buffer200_crs <- st_transform(buffer_200, crs = las_crs)

metrics_buffer200 <- vector("list", length = nrow(buffer200_crs))

for (i in seq_len(nrow(buffer200_crs))) {
  poly <- buffer200_crs[i, ]
  las_clip <- clip_roi(las, poly)
  if (!lidR::is.empty(las_clip)) {
    metrics_buffer200[[i]] <- fsd_metrics(las_clip$Z)
  } else {
    metrics_buffer200[[i]] <- NA
  }
}

# Create a dataframe of the metrics
metrics200_df <- do.call(rbind, metrics_buffer200)
metrics200_df <- as.data.frame(metrics200_df) # Ensures it is a dataframe
metrics200_df$transect_id <- buffer200_crs$Number

# View Metrics  
metrics200_df
      zmean     zmax    zmin zentropy      zsd        zskew transect_id
1  9.769828 24.11743 0.00032 4.165618 5.999804   -0.2268212          12
2  7.976485 23.95768   4e-05  4.02793 5.636253   0.05403732           1
3  10.07398 28.51159   3e-05 4.344737  6.83106   0.05609885           7
4  5.766593 27.13721   5e-05 3.789805  6.17586      1.06807          29
5  9.891049 23.33204 0.00016 4.134766 5.890679   -0.3216369           9
6  11.25073 23.96118   1e-05 4.154001 6.049772   -0.5775392          11
7  8.439339 29.03266   3e-05 4.187206 7.883751    0.6552822          19
8  7.929396 22.25854   1e-05 3.812047 5.574366   -0.1725842           5
9  4.383184 25.06477   1e-04 3.292185 5.619263     1.506778          30
10 8.931003 29.61337   4e-05 4.101234 7.720469     0.374498          13
11 12.93989 35.42897   3e-05 4.300604 9.345639   -0.1996642          22
12 7.833888 21.98303   7e-05 3.953395 5.614427 -0.004073779           4
13 8.519331 27.79898   1e-05 4.243981 6.634762    0.3007393          17
14  8.14094 28.39724   8e-05 3.992443 6.446568    0.1706909          14
15 7.863946 24.70047   8e-05 3.970401 5.749355   0.04419555           3
16 11.82736 33.35335   9e-05 4.416261 9.027765   0.05588564          24
17 9.484697  24.5289 0.00016 4.114216 5.417404   -0.3035339           6
18 9.408753 22.64541 0.00032 4.082647 5.998575   -0.2479153          16
19 9.909541 23.33204 0.00014 4.150613 5.863031   -0.2922431           8
20 7.772781 29.05853   2e-05 4.092662 7.938413    0.8775078          20
21 7.527942 21.99996 0.00012  3.96617  5.37596   0.05481689           2
22 8.194641 22.71527   1e-05  3.98311 5.629962  -0.04313547          15
23 4.851119 19.60874 0.00011 3.556516 4.638608    0.8155331          10
24  9.80478 33.46384   6e-05 4.361991 8.142519    0.4357119          23
# Joins the metrics to the buffer200 shapefile
height_buffer200 <- cbind(buffer200_crs, metrics200_df)

height_buffer200 <- height_buffer200 %>%
  rename(transect = Number) %>%
  select(-transect_id)

Canopy Gaps and Cover

Parameters for detection

threshold <- 2 # Looks for gaps with max canopy height of 2m
size <-c (1,1000) # Looks for gap sizes between 1-1000m2
size1 <- c(1,50000) # Testing larger size for better landscape scale measurement
## Tested, size1 better, detects more gaps.

Per transect

# Crop CHM to transects
chm_transects <- raster(mask(crop(chm, vect(transects)), vect(transects)))

# Check alignment
plot(chm_transects, col = viridis(10))
plot(st_geometry(transects), add = TRUE, border = "blue", lwd = 2)

# Aligned

# Detecting forest gaps
gaps_duc <- getForestGaps(chm_layer=chm_transects, threshold=threshold, size=size1)

# Plotting chm
plot(chm_transects, col=viridis(10))
# Plotting gaps
plot(gaps_duc, col="red", add=TRUE, main="Forest Canopy Gap", legend=FALSE)

# Get stats on gaps
gaps_stats <- GapStats(gap_layer=gaps_duc, chm_layer=chm_transects)

# Converting to SF
gaps_sf <- st_as_sf(GapSPDF(gaps_duc)) # Convert to Sf
gaps_sf <- st_transform(gaps_sf, "EPSG:27700") # ensures in correct crs
                   
# Join transect number to gaps
gaps_transects <- st_join(gaps_sf, transects["Number"], join = st_intersects)

# Clear space
rm(gaps_sf)

# Merge Files
gaps <- merge(gaps_transects, gaps_stats, by="gap_id")
gaps <- gaps %>%
  rename(transect = Number)

# Get area of transect for canopy cover calc
st_area(transects)
Units: [m^2]
 [1] 500.0000 496.0833 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
 [9] 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
[17] 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000 500.0000
# Summarise per transect
canopygaps_transects <- gaps %>%
  group_by(transect) %>%
  summarise(
    num_gaps = n(),
    total_gaparea = sum(gap_area),
    mean_gaparea = mean(gap_area),
    min_gaparea = min(gap_area),
    max_gaparea = max(gap_area),
    sd_gaparea = sd(gap_area),
    skew_gaparea = if (n() >= 3) e1071::skewness(gap_area, type = 2) else NA,
    canopy_cover_perc = (500 - total_gaparea) / 500 * 100) %>%
  mutate(cv_area = sd_gaparea / mean_gaparea) %>%
  arrange(transect)
# Rearranging columns for simplicity
canopygaps_transects <- canopygaps_transects[, c(1,2,3,4,5,6,7,8,11,9,10)]

# test one plot
ggplot() +
  geom_tile(data = chm_df, aes(x = x, y = y, fill = Z)) + # tile to save memory
  scale_fill_viridis_c() +
  geom_sf(data = transects[transects$Number == 1, ], color = "green", fill = NA, size = 1) +
  geom_sf(data = gaps_transects[gaps_transects$Number == 1, ], fill = "red", color = NA, alpha = 0.6) +
  coord_sf(xlim = st_bbox(transects[transects$Number == 1, ])[c("xmin", "xmax")],
           ylim = st_bbox(transects[transects$Number == 1, ])[c("ymin", "ymax")], datum = NA) +
  labs(title = "Transect 1 Canopy Gaps") +
  theme_minimal()

Per 50m Buffer

# Crop CHM to buffer
chm_buffer50 <- raster(mask(crop(chm, vect(buffer_50)), vect(buffer_50)))

# Check alignment
plot(chm_buffer50, col = viridis(10))
plot(st_geometry(buffer_50), add = TRUE, border = "blue", lwd = 2)

# Aligned

# Detecting forest gaps
gaps_buffer50 <- getForestGaps(chm_layer=chm_buffer50, threshold=threshold, size=size1)

# Plotting chm
plot(chm_buffer50, col=viridis(10))
# Plotting gaps
plot(gaps_buffer50, col="red", add=TRUE, main="Forest Canopy Gap", legend=FALSE)

# Get stats on gaps
gaps_stats50 <- GapStats(gap_layer=gaps_buffer50, chm_layer=chm_buffer50)

# Converting to SF
gaps50_sf <- st_as_sf(GapSPDF(gaps_buffer50)) # Convert to Sf
gaps50_sf <- st_transform(gaps50_sf, "EPSG:27700") # ensures in correct crs

# Join transect number to gaps
gaps_buffer50 <- st_join(gaps50_sf, buffer_50["Number"], join = st_intersects)

# Clear space
rm(gaps50_sf)

# Merge Files
gaps50 <- merge(gaps_buffer50, gaps_stats50, by="gap_id")
gaps50 <- gaps50 %>%
  rename(transect = Number)

# Get area of buffer_50 for canopy cover calc
st_area(buffer_50)
Units: [m^2]
 [1] 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393
 [9] 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393
[17] 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393 7850.393
# Summarise per buffer
canopygaps_buffer50 <- gaps50 %>%
  group_by(transect) %>%
  summarise(
    num_gaps = n(),
    total_gaparea = sum(gap_area),
    mean_gaparea = mean(gap_area),
    min_gaparea = min(gap_area),
    max_gaparea = max(gap_area),
    sd_gaparea = sd(gap_area),
    skew_gaparea = if (n() >= 3) e1071::skewness(gap_area, type = 2) else NA,
    canopy_cover_perc = (7850.393 - total_gaparea) / 7850.393 * 100) %>%
  mutate(cv_area = sd_gaparea / mean_gaparea) %>%
  arrange(transect)
# Rearranging columns for simplicity
canopygaps_buffer50 <- canopygaps_buffer50[, c(1,2,3,4,5,6,7,8,11,9,10)]

view(canopygaps_buffer50)

# test one plot
ggplot() +
  geom_tile(data = chm_df, aes(x = x, y = y, fill = Z)) + # tile to save memory
  scale_fill_viridis_c() +
  geom_sf(data = buffer_50[buffer_50$Number == 1, ], color = "green", fill = NA, size = 1) +
  geom_sf(data = gaps_buffer50[gaps_buffer50$Number == 1, ], fill = "red", color = NA, alpha = 0.6) +
  coord_sf(xlim = st_bbox(buffer_50[buffer_50$Number == 1, ])[c("xmin", "xmax")],
           ylim = st_bbox(buffer_50[buffer_50$Number == 1, ])[c("ymin", "ymax")], datum = NA) +
  labs(title = "Transect 1 Canopy Gaps Buffer50") +
  theme_minimal()

Per 200m Buffer - ISSUES WITH THIS RANGE DEMONSTRATED IN NEXT SECTION

# Crop CHM to buffer
chm_buffer200 <- raster(mask(crop(chm, vect(buffer_200)), vect(buffer_200)))

# Check alignment
plot(chm_buffer200, col = viridis(10))
plot(st_geometry(buffer_200), add = TRUE, border = "blue", lwd = 2)

# Aligned

# Detecting forest gaps
gaps_duc200 <- getForestGaps(chm_layer=chm_buffer200, threshold=threshold, size=size1)

# Get stats on gaps
gaps_stats200 <- GapStats(gap_layer=gaps_duc200, chm_layer=chm_buffer200)

# Plotting chm
plot(chm_buffer200, col=viridis(10))
# Plotting gaps
plot(gaps_duc200, col="red", add=TRUE, main="Forest Canopy Gap", legend=FALSE)

# Converting to SF
gaps200_sf <- st_as_sf(GapSPDF(gaps_duc200)) # Convert to Sf
gaps200_sf <- st_transform(gaps200_sf, "EPSG:27700") # ensures in correct crs

# Join transect number to gaps
gaps_buffer200 <- st_join(gaps200_sf, buffer_200["Number"], join = st_intersects)


#---------------------------------------------------------------------------------

# Crop gaps to buffer shape
# This catches those gaps which extend into another nearby buffer as some 200m buffers overlap

# Create cropped geometries
cropped_geoms <- pmap(
  list(gap_geom = gaps_buffer200$geometry, buffer_id = gaps_buffer200$Number),
  function(gap_geom, buffer_id) {
    st_intersection(
      gap_geom,
      st_geometry(buffer_200[buffer_200$Number == buffer_id, ]) %>%
        st_set_crs(st_crs(gap_geom)))})

gaps_cropped <- gaps_buffer200 %>%
  st_set_geometry(NULL) %>%                         # strip old geometry
  mutate(geometry = st_sfc(cropped_geoms)) %>%
  st_as_sf(crs = st_crs(gaps_buffer200))

# Check if cropping has worked
ggplot() +
  geom_sf(data = gaps_cropped[gaps_cropped$Number == 16,], fill = "red") +
  geom_sf(data = buffer_200[buffer_200$Number == 16,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 11,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 17,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 12,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 15,], fill = "transparent") +
  coord_sf()

## CROPPING HAS WORKED BUT CANNOT GET STATS ON NEW FILE AS IT IS NOT CORRECT FORMAT


#-------------------------------------------------------------------------------


# Merge Files
gaps200 <- merge(gaps_buffer200, gaps_stats200, by="gap_id")
gaps200 <- gaps200 %>%
  rename(transect = Number)

# Get area of buffer_200 for canopy cover calc
st_area(buffer_200)
Units: [m^2]
 [1] 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3
 [9] 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3
[17] 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3 125606.3
# Summarise per transect
canopygaps_buffer200 <- gaps200 %>%
  group_by(transect) %>%
  summarise(
    num_gaps = n(),
    total_gaparea = sum(gap_area),
    mean_gaparea = mean(gap_area),
    min_gaparea = min(gap_area),
    max_gaparea = max(gap_area),
    sd_gaparea = sd(gap_area),
    skew_gaparea = if (n() >= 3) e1071::skewness(gap_area, type = 2) else NA,
    canopy_cover_perc = (125606.3 - total_gaparea) / 125606.3 * 100) %>%
  mutate(cv_area = sd_gaparea / mean_gaparea) %>%
  arrange(transect)
# Rearranging columns for simplicity
canopygaps_buffer200 <- canopygaps_buffer200[, c(1,2,3,4,5,6,7,8,11,9,10)]

view(canopygaps_buffer200)

# test one plot
ggplot() +
  geom_tile(data = chm_df, aes(x = x, y = y, fill = Z)) + # tile to save memory
  scale_fill_viridis_c() +
  geom_sf(data = buffer_200[buffer_200$Number == 4, ], color = "green", fill = NA, size = 1) +
  geom_sf(data = gaps_buffer200[gaps_buffer200$Number == 4, ], fill = "red", color = NA, alpha = 0.6) +
  coord_sf(xlim = st_bbox(buffer_200[buffer_200$Number == 4, ])[c("xmin", "xmax")],
           ylim = st_bbox(buffer_200[buffer_200$Number == 4, ])[c("ymin", "ymax")], datum = NA) +
  labs(title = "Transect 1 Canopy Gaps Buffer200") +
  theme_minimal()

Issues

Some 200m buffers overlap. Using the normal method, gaps that overlap into a a buffer are entirely included even if only a small part is within the buffer. This means that the buffer area is higher than it should be.

This only occurs in the areas where buffers overlap.

I successfully fixed this by cropping the gaps to gaps_cropped. However, to get stats on this file, it has to be directly from the Getforestgaps() function.

I am unable to get the stats for the 200m buffers correctly.

ggplot() +
  geom_sf(data = gaps_buffer200[gaps_buffer200$Number == 16,], fill = "red") +
  geom_sf(data = buffer_200[buffer_200$Number == 16,], fill = "transparent") +
  coord_sf()

ggplot() +
  geom_sf(data = gaps_buffer200[gaps_buffer200$Number == 16,], fill = "red") +
  geom_sf(data = buffer_200[buffer_200$Number == 16,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 11,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 17,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 12,], fill = "transparent") +
  geom_sf(data = buffer_200[buffer_200$Number == 15,], fill = "transparent") +
  coord_sf()