AGU 2022 GEDI Workshop
Studying forest structure from space: An introduction to GEDI lidar data pre-processing and analysis in ecological applications
Code authors: Melissa Rose, Ivan Gonzalez, Jenna Keany \(\\\)\(\\\)
Technical contributors: Emma Reich, Shelby Sunquist, Patrick Jantz
The purpose of this workshop is to assist interested Ecologists, Biologists, and Geologists in accessing NASA’s Global Ecosystem Dynamics Investigation (GEDI) near-global lidar data and applying it to their research. We will review recent applications and discuss future directions for GEDI data. Using open-access R software, we will demonstrate how participants may easily access, manipulate, explore, assess, and apply GEDI data in their study areas. We will develop modules (download > cleaning > visualization, etc) and share code, inputs, and outputs to prevent bottlenecks in case of technical issues. This strategy allows attendees to participate in every module with products derived from previous stages even if previous steps/results weren’t completed. Participants may be required to install open-source software prior to the workshop. This workshop is organized primarily by PhD students in the Ecological and Environmental Informatics program at Northern Arizona University.
Objectives
The primary objective of this workshop is to disseminate information on the use and functionality of spaceborne lidar (light detection and ranging) from GEDI, including: how the instrument functions, spatial resolution and coverage, available data products, where and how to access the data, and uses of lidar data in a variety of earth science subjects.
Participants will learn how to download GEDI data for an area of interest, process the data, and apply it to an ecological question in R.
Presenters will provide example code and intermediate data products so that participants can follow along and explore results in the event of technical challenges. \(\\\)\(\\\)\(\\\)\(\\\)
To launch an interactive coding environment:
To download this code, please refer to the workshop github:
Tutorial:
SECTION 0: PREPARE R ENVIRONMENT
0.1: Libraries
Do not need to run this chunk when running in Binder
#Install libraries if not already installed
#libs <- c('raster', 'devtools', 'sf', 'sp', 'data.table', 'dplyr', 'lidR',
# 'httr', 'tidyr', 'curl', 'hdf5r', 'leaflet', 'leaflet.extras',
# 'leafsync', 'ggplot2', 'rgeos', 'RColorBrewer', 'rasterVis',
# 'viridis', 'rdgal', 'rhdf5', 'dismo', 'sdm', 'randomForest', 'hdf5r')
#sapply(libs, FUN = function(x){
# if (!x %in% installed.packages()){
# install.packages(x)
# }
#})
#sdm::installAll()#Define which libraries to use
library(raster) ## spatial tools
library(devtools) ## tools for building and installing libraries
library(sf) ## spatial tools
library(sp) ## spatial tools
library(data.table) ## data / table manipulation
library(lidR) ## lidar manipulation tools
library(httr) ## methods for web files interaction
library(tidyr) ## data / table manipulation
library(curl) ## download files from the web
library(hdf5r) ## reading hdf5 files
library(leaflet) ## color palettes
library(leafsync) ## interactive maps
library(ggplot2) ## nice plot tools
library(rgeos) ## geometric for vector data
library(RColorBrewer) ## color palettes
library(rasterVis) ## raster viz tools
library(viridis) ## color palettes
library(rgdal) ## readOGR: Read OGR vector maps into Spatial objects
library(rhdf5) ## R Interface to HDF5
library(dismo) ## Species Distribution Modeling
library(sdm) ## Species Distribution Modelling
library(randomForest) ## For sdm methd 'rf'
library(hdf5r) ## To show h5 file structure
#library(rGEDI) ## GEDI functions in R0.2: Directory and Subdirectories
#define folder paths
wd <- getwd()
gL1B_dir <- './data/GEDI01_B.002/'
gL2A_dir <- './data/GEDI02_A.002/'
gL2B_dir <- './data/GEDI02_B.002/'0.3: Source Scripts
#source rGEDI functions in scripts directory
source('./scripts/gedi_functions.R')SECTION 1: DOWNLOAD GEDI DATA FOR A REGION OF INTEREST
1.1: Define region of interest and date range
#ROI: Colombia
#specify bounding box coordinates for ROI
ll_lon <- -75 #lower left longitude
ll_lat <- 0 #lower left latitude
ur_lon <- -74 #upper right longitude
ur_lat <- 1 #upper right latitude
#convert into bounding box coordinates
bbox <- paste(ll_lon, ll_lat, ur_lon, ur_lat, sep = ',')
#specify the date range to filter orbits by
daterange <- c("2020-05-01","2020-07-31")1.2: Visualize region of interest
#view footprint locations for entire orbit on map
leaflet() %>%
addRectangles(
lng1=ll_lon, lat1=ll_lat,
lng2=ur_lon, lat2=ur_lat,
color = "grey",
fillColor = "transparent") %>%
addScaleBar(options = list(imperial = FALSE)) %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addLegend(colors = "grey", labels = "Region of Interest") 1.3: Use gedi_finder to find overlapping GEDI orbits
#create a list of all overlapping orbits for each GEDI data product
gL1B <- gedi_finder('GEDI01_B.002', bbox) ## [1] "Found 157 GEDI01_B.002 overlapping orbits found."gL2A <- gedi_finder('GEDI02_A.002', bbox)## [1] "Found 157 GEDI02_A.002 overlapping orbits found."gL2B <- gedi_finder('GEDI02_B.002', bbox)## [1] "Found 157 GEDI02_B.002 overlapping orbits found."1.4 Use gedi_temp_filter to subset orbit list within date range
#subset orbit list using daterange
gL1B <- gedi_temp_filter(gL1B, 'GEDI01_B.002', daterange)## [1] "Found 21 GEDI01_B.002 overlapping orbits from 2020-05-01 to 2020-07-31."gL2A <- gedi_temp_filter(gL2A, 'GEDI02_A.002', daterange)## [1] "Found 21 GEDI02_A.002 overlapping orbits from 2020-05-01 to 2020-07-31."gL2B <- gedi_temp_filter(gL2B, 'GEDI02_B.002', daterange)## [1] "Found 21 GEDI02_B.002 overlapping orbits from 2020-05-01 to 2020-07-31."Stop here ——————————- 1.5: Use gedi_download to download GEDI data (DO NOT RUN)
#download data product to their output data directory
# gedi_download(filepath=gL1B,outdir=gL1B_dir)
# gedi_download(filepath=gL2A,outdir=gL2A_dir)
# gedi_download(filepath=gL2B,outdir=gL2B_dir)SECTION 2: CLIP GEDI DATA PRODUCTS TO REGION OF INTEREST AND FILTER FOR QUALITY SHOTS 2.1 Read GEDI data (DO NOT RUN)
#read the level 1B, 2A, and 2B for a single orbit
## Remove the "#" symbol in the next lines in order to run them
# gedilevel1B <- readLevel1B(level1Bpath = paste0(gL1B_dir,"GEDI01_B_2020160005619_O08428_03_T01337_02_005_01_V002.h5"))
# gedilevel2A <- readLevel2A(level2Apath = paste0(gL2A_dir,"GEDI02_A_2020160005619_O08428_03_T01337_02_003_01_V002.h5"))
# gedilevel2B <- readLevel2B(level2Bpath = paste0(gL2B_dir,"GEDI02_B_2020160005619_O08428_03_T01337_02_003_01_V002.h5"))2.2 Clip data within bbox coordinates (DO NOT RUN)
## Remove the "#" symbol in the next lines in order to run them
# gedilevel1B <- clipLevel1B(gedilevel1B, ll_lon, ur_lon, ll_lat, ur_lat, output = paste0(gL1B_dir,"GEDI01_B_2020160005619_O08428_03_T01337_02_005_01_V002_clip.h5"))
# gedilevel2A <- clipLevel2A(gedilevel2A, ll_lon, ur_lon, ll_lat, ur_lat, output = paste0(gL2A_dir,"GEDI02_A_2020160005619_O08428_03_T01337_02_003_01_V002_clip.h5"))
# gedilevel2B <- clipLevel2B(gedilevel2B, ll_lon, ur_lon, ll_lat, ur_lat, output = paste0(gL2B_dir,"GEDI02_B_2020160005619_O08428_03_T01337_02_003_01_V002_clip.h5"))Restart here ——————————-
Read clipped files into R
gedilevel1B <- readLevel1B(level1Bpath = paste0(gL1B_dir,"GEDI01_B_2020160005619_O08428_03_T01337_02_005_01_V002_clip.h5"))
gedilevel2A <- readLevel2A(level2Apath = paste0(gL2A_dir,"GEDI02_A_2020160005619_O08428_03_T01337_02_003_01_V002_clip.h5"))
gedilevel2B <- readLevel2B(level2Bpath = paste0(gL2B_dir,"GEDI02_B_2020160005619_O08428_03_T01337_02_003_01_V002_clip.h5"))2.3 Use GEDI L2A quality flag to select only quality shots
#view GEDI L2A h5 layers
h5groups2A <- hdf5r::list.groups(gedilevel2A@h5[['BEAM0000']])
h5groups2A## [1] "ancillary" "geolocation"
## [3] "land_cover_data" "rx_1gaussfit"
## [5] "rx_1gaussfit/ancillary" "rx_assess"
## [7] "rx_assess/ancillary" "rx_processing_a1"
## [9] "rx_processing_a1/ancillary" "rx_processing_a2"
## [11] "rx_processing_a2/ancillary" "rx_processing_a3"
## [13] "rx_processing_a3/ancillary" "rx_processing_a4"
## [15] "rx_processing_a4/ancillary" "rx_processing_a5"
## [17] "rx_processing_a5/ancillary" "rx_processing_a6"
## [19] "rx_processing_a6/ancillary"#view GEDI L2A h5 layers
h5datasets2A <- hdf5r::list.datasets(gedilevel2A@h5[['BEAM0000']])
head(h5datasets2A, 100)## [1] "ancillary/l2a_alg_count" "beam"
## [3] "channel" "degrade_flag"
## [5] "delta_time" "digital_elevation_model"
## [7] "digital_elevation_model_srtm" "elev_highestreturn"
## [9] "elev_lowestmode" "elevation_bias_flag"
## [11] "elevation_bin0_error" "energy_total"
## [13] "geolocation/elev_highestreturn_a1" "geolocation/elev_highestreturn_a2"
## [15] "geolocation/elev_highestreturn_a3" "geolocation/elev_highestreturn_a4"
## [17] "geolocation/elev_highestreturn_a5" "geolocation/elev_highestreturn_a6"
## [19] "geolocation/elev_lowestmode_a1" "geolocation/elev_lowestmode_a2"
## [21] "geolocation/elev_lowestmode_a3" "geolocation/elev_lowestmode_a4"
## [23] "geolocation/elev_lowestmode_a5" "geolocation/elev_lowestmode_a6"
## [25] "geolocation/elev_lowestreturn_a1" "geolocation/elev_lowestreturn_a2"
## [27] "geolocation/elev_lowestreturn_a3" "geolocation/elev_lowestreturn_a4"
## [29] "geolocation/elev_lowestreturn_a5" "geolocation/elev_lowestreturn_a6"
## [31] "geolocation/elevation_1gfit" "geolocation/elevs_allmodes_a1"
## [33] "geolocation/elevs_allmodes_a2" "geolocation/elevs_allmodes_a3"
## [35] "geolocation/elevs_allmodes_a4" "geolocation/elevs_allmodes_a5"
## [37] "geolocation/elevs_allmodes_a6" "geolocation/energy_lowestmode_a1"
## [39] "geolocation/energy_lowestmode_a2" "geolocation/energy_lowestmode_a3"
## [41] "geolocation/energy_lowestmode_a4" "geolocation/energy_lowestmode_a5"
## [43] "geolocation/energy_lowestmode_a6" "geolocation/lat_highestreturn_a1"
## [45] "geolocation/lat_highestreturn_a2" "geolocation/lat_highestreturn_a3"
## [47] "geolocation/lat_highestreturn_a4" "geolocation/lat_highestreturn_a5"
## [49] "geolocation/lat_highestreturn_a6" "geolocation/lat_lowestmode_a1"
## [51] "geolocation/lat_lowestmode_a2" "geolocation/lat_lowestmode_a3"
## [53] "geolocation/lat_lowestmode_a4" "geolocation/lat_lowestmode_a5"
## [55] "geolocation/lat_lowestmode_a6" "geolocation/lat_lowestreturn_a1"
## [57] "geolocation/lat_lowestreturn_a2" "geolocation/lat_lowestreturn_a3"
## [59] "geolocation/lat_lowestreturn_a4" "geolocation/lat_lowestreturn_a5"
## [61] "geolocation/lat_lowestreturn_a6" "geolocation/latitude_1gfit"
## [63] "geolocation/lats_allmodes_a1" "geolocation/lats_allmodes_a2"
## [65] "geolocation/lats_allmodes_a3" "geolocation/lats_allmodes_a4"
## [67] "geolocation/lats_allmodes_a5" "geolocation/lats_allmodes_a6"
## [69] "geolocation/lon_highestreturn_a1" "geolocation/lon_highestreturn_a2"
## [71] "geolocation/lon_highestreturn_a3" "geolocation/lon_highestreturn_a4"
## [73] "geolocation/lon_highestreturn_a5" "geolocation/lon_highestreturn_a6"
## [75] "geolocation/lon_lowestmode_a1" "geolocation/lon_lowestmode_a2"
## [77] "geolocation/lon_lowestmode_a3" "geolocation/lon_lowestmode_a4"
## [79] "geolocation/lon_lowestmode_a5" "geolocation/lon_lowestmode_a6"
## [81] "geolocation/lon_lowestreturn_a1" "geolocation/lon_lowestreturn_a2"
## [83] "geolocation/lon_lowestreturn_a3" "geolocation/lon_lowestreturn_a4"
## [85] "geolocation/lon_lowestreturn_a5" "geolocation/lon_lowestreturn_a6"
## [87] "geolocation/longitude_1gfit" "geolocation/lons_allmodes_a1"
## [89] "geolocation/lons_allmodes_a2" "geolocation/lons_allmodes_a3"
## [91] "geolocation/lons_allmodes_a4" "geolocation/lons_allmodes_a5"
## [93] "geolocation/lons_allmodes_a6" "geolocation/num_detectedmodes_a1"
## [95] "geolocation/num_detectedmodes_a2" "geolocation/num_detectedmodes_a3"
## [97] "geolocation/num_detectedmodes_a4" "geolocation/num_detectedmodes_a5"
## [99] "geolocation/num_detectedmodes_a6" "geolocation/quality_flag_a1"#read select GEDI L2A layers into data.table
level2AM <- getLevel2AM(gedilevel2A)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%head(level2AM, 5)## beam shot_number degrade_flag quality_flag delta_time sensitivity
## 1: BEAM0000 84280000300247456 0 1 76817129 0.9446448
## 2: BEAM0000 84280000300247457 0 1 76817129 0.9768031
## 3: BEAM0000 84280000300247458 0 1 76817129 0.9755690
## 4: BEAM0000 84280000300247459 0 1 76817129 0.9304212
## 5: BEAM0000 84280000300247460 0 1 76817129 0.9657014
## sensitivity_a2 surface_flag water_persistence urban_proportion
## 1: 0.9773547 1 0 0
## 2: 0.9768031 1 0 0
## 3: 0.9755690 1 0 0
## 4: 0.9676956 1 0 0
## 5: 0.9657014 1 0 0
## solar_elevation lat_lowestmode lon_lowestmode elev_highestreturn
## 1: -41.60833 0.2771671 -74.99973 215.9003
## 2: -41.60883 0.2767456 -74.99944 217.4226
## 3: -41.60933 0.2763244 -74.99914 221.4247
## 4: -41.60983 0.2759028 -74.99884 221.2567
## 5: -41.61033 0.2754821 -74.99855 227.2185
## elev_lowestmode rh0 rh1 rh2 rh3 rh4 rh5 rh6 rh7 rh8 rh9
## 1: 193.7073 -1.27 -0.85 -0.48 -0.11 0.22 0.59 1.00 1.45 1.98 2.54
## 2: 193.6978 -3.28 -2.46 -1.94 -1.56 -1.30 -1.08 -0.85 -0.67 -0.48 -0.33
## 3: 194.3747 -4.03 -3.21 -2.57 -2.01 -1.56 -1.19 -0.85 -0.52 -0.26 0.03
## 4: 194.4308 -3.62 -3.21 -2.91 -2.65 -2.42 -2.24 -2.09 -1.94 -1.79 -1.64
## 5: 196.0961 -4.25 -3.51 -2.98 -2.42 -1.68 -0.70 0.11 0.97 1.90 2.72
## rh10 rh11 rh12 rh13 rh14 rh15 rh16 rh17 rh18 rh19 rh20 rh21
## 1: 3.02 3.51 3.96 4.40 4.85 5.38 5.90 6.50 7.17 7.77 8.33 8.78
## 2: -0.14 0.00 0.18 0.33 0.52 0.67 0.85 1.00 1.19 1.34 1.53 1.71
## 3: 0.29 0.59 0.89 1.15 1.49 1.86 2.39 2.91 3.36 3.77 4.14 4.59
## 4: -1.53 -1.41 -1.30 -1.19 -1.08 -0.97 -0.89 -0.78 -0.70 -0.59 -0.52 -0.41
## 5: 3.36 4.07 4.85 5.67 6.72 7.92 9.15 10.42 11.47 12.14 12.77 13.37
## rh22 rh23 rh24 rh25 rh26 rh27 rh28 rh29 rh30 rh31 rh32 rh33
## 1: 9.19 9.56 9.90 10.23 10.57 10.87 11.13 11.39 11.65 11.91 12.14 12.32
## 2: 1.90 2.09 2.27 2.50 2.69 2.87 3.10 3.36 3.58 3.88 4.14 4.48
## 3: 5.11 5.86 6.53 7.09 7.54 7.92 8.33 8.70 9.11 9.45 9.78 10.12
## 4: -0.33 -0.22 -0.14 -0.07 0.03 0.11 0.22 0.29 0.37 0.48 0.56 0.67
## 5: 13.97 14.45 14.87 15.24 15.61 15.99 16.40 16.77 17.14 17.52 17.89 18.23
## rh34 rh35 rh36 rh37 rh38 rh39 rh40 rh41 rh42 rh43 rh44 rh45
## 1: 12.55 12.74 12.92 13.11 13.30 13.45 13.63 13.78 13.93 14.08 14.23 14.38
## 2: 4.85 5.23 5.64 6.05 6.42 6.83 7.24 7.69 8.14 8.59 9.00 9.34
## 3: 10.46 10.87 11.32 11.76 12.17 12.59 12.96 13.30 13.59 13.89 14.19 14.45
## 4: 0.74 0.85 0.97 1.04 1.15 1.27 1.38 1.53 1.64 1.79 1.94 2.09
## 5: 18.53 18.79 19.05 19.27 19.46 19.65 19.83 20.02 20.21 20.39 20.62 20.81
## rh46 rh47 rh48 rh49 rh50 rh51 rh52 rh53 rh54 rh55 rh56 rh57
## 1: 14.53 14.68 14.83 14.98 15.13 15.24 15.39 15.54 15.69 15.84 15.95 16.10
## 2: 9.63 9.90 10.16 10.38 10.61 10.83 11.05 11.32 11.58 11.91 12.25 12.66
## 3: 14.72 15.01 15.28 15.54 15.80 16.06 16.32 16.58 16.88 17.14 17.37 17.59
## 4: 2.27 2.50 2.72 2.98 3.32 3.69 4.14 4.63 5.15 5.71 6.23 6.76
## 5: 20.99 21.18 21.37 21.55 21.74 21.89 22.08 22.23 22.37 22.52 22.71 22.86
## rh58 rh59 rh60 rh61 rh62 rh63 rh64 rh65 rh66 rh67 rh68 rh69
## 1: 16.25 16.36 16.51 16.62 16.77 16.88 16.99 17.11 17.22 17.33 17.44 17.56
## 2: 13.00 13.33 13.63 13.89 14.16 14.38 14.60 14.83 15.01 15.20 15.43 15.61
## 3: 17.78 17.97 18.19 18.38 18.56 18.79 19.01 19.20 19.42 19.65 19.83 20.02
## 4: 7.24 7.73 8.25 8.89 9.86 11.58 13.37 14.94 15.95 16.66 17.26 17.70
## 5: 23.01 23.16 23.31 23.46 23.61 23.76 23.91 24.06 24.24 24.39 24.54 24.69
## rh70 rh71 rh72 rh73 rh74 rh75 rh76 rh77 rh78 rh79 rh80 rh81
## 1: 17.67 17.74 17.85 17.93 18.04 18.15 18.23 18.34 18.41 18.53 18.60 18.71
## 2: 15.80 15.99 16.21 16.40 16.58 16.77 16.96 17.14 17.33 17.52 17.70 17.89
## 3: 20.21 20.39 20.54 20.69 20.84 20.99 21.14 21.29 21.44 21.59 21.74 21.89
## 4: 18.15 18.53 18.86 19.20 19.50 19.80 20.10 20.36 20.62 20.88 21.14 21.37
## 5: 24.84 24.99 25.10 25.25 25.40 25.55 25.70 25.85 26.00 26.15 26.34 26.48
## rh82 rh83 rh84 rh85 rh86 rh87 rh88 rh89 rh90 rh91 rh92 rh93
## 1: 18.83 18.90 19.01 19.12 19.20 19.31 19.42 19.54 19.65 19.80 19.91 20.06
## 2: 18.08 18.26 18.41 18.60 18.75 18.90 19.09 19.24 19.42 19.65 19.83 20.10
## 3: 22.08 22.23 22.41 22.60 22.79 22.97 23.16 23.35 23.53 23.76 23.98 24.21
## 4: 21.59 21.81 22.04 22.26 22.49 22.71 22.90 23.12 23.35 23.57 23.79 24.02
## 5: 26.63 26.82 26.97 27.16 27.31 27.49 27.68 27.83 28.02 28.20 28.39 28.58
## rh94 rh95 rh96 rh97 rh98 rh99 rh100
## 1: 20.21 20.39 20.58 20.84 21.10 21.52 22.19
## 2: 20.39 20.81 21.37 21.93 22.41 22.90 23.72
## 3: 24.43 24.69 24.95 25.25 25.66 26.26 27.05
## 4: 24.28 24.54 24.84 25.18 25.55 26.04 26.82
## 5: 28.80 29.03 29.25 29.51 29.85 30.26 31.12#convert shot_number from integer to character
level2AM$shot_number <- paste0(level2AM$shot_number)
dim(level2AM)## [1] 3403 116#filtering for quality_flag == 1
level2AM_qual <- level2AM[quality_flag == 1]
dim(level2AM_qual)## [1] 3323 1162.4 Visualize location of quality_flag filtered GEDI footprints
#view footprint locations on map
leaflet() %>%
addRectangles(
lng1=ll_lon, lat1=ll_lat,
lng2=ur_lon, lat2=ur_lat,
color = "grey",
fillColor = "transparent") %>%
addCircleMarkers(level2AM[quality_flag == 0, lon_lowestmode],
level2AM[quality_flag == 0, lat_lowestmode],
radius = 1,
opacity = 0.5,
color = "red",
popup = paste0("shot_number: ", level2AM$shot_number)) %>%
addCircleMarkers(level2AM_qual$lon_lowestmode,
level2AM_qual$lat_lowestmode,
radius = 1,
opacity = 0.5,
color = "blue",
popup = paste0("shot_number: ", level2AM_qual$shot_number)) %>%
setView(median(level2AM_qual$lon_lowestmode),
median(level2AM_qual$lat_lowestmode),
zoom = 10) %>%
addScaleBar(options = list(imperial = FALSE)) %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addLegend(colors = c("grey", "blue", "red"), labels = c("Region of Interest", "Quality GEDI Footprints", "Non-Quality GEDI Footprints"))2.5 Filter beams to select only Power Beams
#view beams based on their type
#BEAM0000,BEAM0001,BEAM0010,BEAM0011 are Coverage Beams.
#BEAM0101,BEAM0110,BEAM1000,BEAM1011 are Full Power Beams
level2AM_qual$beamType <- ifelse(level2AM_qual$beam %in% c('BEAM0000','BEAM0001','BEAM0010','BEAM0011'), 'Coverage', 'Power')
#view footprint locations on map
leaflet() %>%
addCircleMarkers(level2AM_qual$lon_lowestmode,
level2AM_qual$lat_lowestmode,
radius = 1,
opacity = 0.5,
color = c('red', 'blue')[as.numeric(as.factor(level2AM_qual$beamType))],
popup = paste0(level2AM_qual$beamType, ' beam - ', level2AM_qual$beam,
'<br>', 'shot_number: ', level2AM_qual$shot_number )) %>%
addScaleBar(options = list(imperial = FALSE)) %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addLegend(colors = c( "blue", "red"),
labels = c("Power beam", "Coverage beam"))#plot power and coverage shots ID from previous interactive map
plotVertProfiles <- plot_rh_profiles(level2AM_qual[shot_number %in%
c("84280500300243078", "84280300300243183")], # power and coverage
idcolumn = 'beamType', col_vector = c('blue', 'red'))
#subset data frame to only include Power Beam shots
dim(level2AM_qual)## [1] 3323 117level2AM_qual <- level2AM_qual[beamType == 'Power']
dim(level2AM_qual)## [1] 1496 1172.6 Further filter shots using region-specific quality parameters
#using a set of other region-specific quality parameters
dim(level2AM_qual)## [1] 1496 117level2AM_qual <- level2AM_qual[sensitivity_a2 > 0.95]
level2AM_qual <- level2AM_qual[rh100 <= 50]
level2AM_qual <- level2AM_qual[surface_flag == 1]
level2AM_qual <- level2AM_qual[water_persistence < 10]
level2AM_qual <- level2AM_qual[urban_proportion < 50]
dim(level2AM_qual)## [1] 1491 1172.7 Visualize final filtered GEDI L2A dataset
#view footprint locations on map
leaflet() %>%
addCircleMarkers(level2AM_qual$lon_lowestmode,
level2AM_qual$lat_lowestmode,
radius = 1,
opacity = 0.5,
color = 'blue',
popup = paste0('shot_number: ', level2AM_qual$shot_number )) %>%
addScaleBar(options = list(imperial = FALSE)) %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addLegend(colors = "blue",
labels = "Final Filtered Shots")2.8 Select individual footprints along a structurally heterogeneous gradient
sample_shots <- c("84280500300242755", "84280500300242756", "84280500300242757", "84280500300242758", "84280500300242759")
sample_region <- level2AM_qual[shot_number %in% sample_shots]
sample_region$sortid <- 1:nrow(sample_region)
#view footprint in sample region
colorRamp <- colorNumeric( palette = "viridis", domain = sample_region$sortid, reverse = TRUE)
pointColors <- colorRamp(sample_region$sortid)
leaflet() %>%
addRectangles(
lng1=min(sample_region$lon_lowestmode)-0.00025, lat1=min(sample_region$lat_lowestmode)-0.00025,
lng2=max(sample_region$lon_lowestmode)+0.00025, lat2=max(sample_region$lat_lowestmode)+0.00025,
color = "grey",
fillColor = "transparent") %>%
addCircleMarkers(sample_region$lon_lowestmode,
sample_region$lat_lowestmode,
radius = 25,
weight = 5,
fillColor = "transparent",
color = pointColors,
popup = paste0("shot_number: ", sample_region$shot_number)) %>%
addScaleBar(options = list(imperial = FALSE)) %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addLegend(colors = c("grey", "blue"), labels = c("Sample Region", "Selected Shots")) Now lets see their vertical profiles
plotVertProfiles <- plot_rh_profiles(sample_region, idcolumn = 'shot_number', col_vector = pointColors)
SECTION 3: Visualize GEDI Data Products
3.1 Plot GEDIL1B Full Waveform
#select shot_numbers for different cover types
shotnum_forest <- "84280500300242755"
shotnum_edge <- "84280500300242757"
shotnum_grass <- "84280500300242759"
sample_shotnums <- c(shotnum_forest, shotnum_edge, shotnum_grass)
#get waveform data for each selected footprint
wf_forest <- getLevel1BWF(gedilevel1B, shot_number=shotnum_forest)
wf_edge <- getLevel1BWF(gedilevel1B, shot_number=shotnum_edge)
wf_grass <- getLevel1BWF(gedilevel1B, shot_number=shotnum_grass)
#check shot's data structure
head(wf_forest@dt)## rxwaveform elevation
## 1: 205.6654 282.0606
## 2: 206.3967 281.9114
## 3: 207.3341 281.7621
## 4: 208.4659 281.6128
## 5: 209.7527 281.4636
## 6: 211.0228 281.3143#plot waveforms
par(mfrow = c(1,1), mar=c(4,4,1,1), cex.axis = 1.5)
plot(wf_forest, relative=TRUE, polygon=TRUE, type="l", lwd=2, col="#FDE725",
xlab="Waveform Amplitude(%)", ylab="Elevation (m)", main="GEDIL1B rxwaveform: Forest", ylim=c(200,260))
plot(wf_edge, relative=TRUE, polygon=TRUE, type="l", lwd=2, col="#21908D",
xlab="Waveform Amplitude(%)", ylab="Elevation (m)", main="GEDIL1B rxwaveform: Edge", ylim=c(200,260), add = TRUE)
plot(wf_grass, relative=TRUE, polygon=TRUE, type="l", lwd=2, col="#440154",
xlab="Waveform Amplitude(%)", ylab="Elevation (m)", main="GEDIL1B rxwaveform: Grass", ylim=c(200,260), add = TRUE)
3.2 Plot GEDIL2A Elevation and Height Metrics
#plot waveform with GEDIL2A relative height metrics
plotWFMetrics(gedilevel1B, gedilevel2A, shotnum_forest, rh=c(25, 50, 75, 90), customylim = c(200, 240), main = 'Forest', add = TRUE, colBG = "#FDE725")##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%
plotWFMetrics(gedilevel1B, gedilevel2A, shotnum_edge, rh=c(25, 50, 75, 90), customylim = c(200, 240), main = 'Edge', add = TRUE, colBG = "#21908D")##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%
plotWFMetrics(gedilevel1B, gedilevel2A, shotnum_grass, rh=c(25, 50, 75, 90), customylim = c(200, 240), main = 'grass', add = TRUE, colBG = "#440154")##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%
3.3 Plot GEDIL2B Plant Area Index (PAI) and Plant Area Volume Density (PAVD)
#get canopy cover and vertical profile metrics
level2BVPM <- getLevel2BVPM(gedilevel2B)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%#check the data structure and 3 first values
head(level2BVPM, 3) ## beam shot_number algorithmrun_flag l2b_quality_flag delta_time
## 1: BEAM0000 84280000300247456 1 1 76817129
## 2: BEAM0000 84280000300247457 1 1 76817129
## 3: BEAM0000 84280000300247458 1 1 76817129
## sensitivity solar_elevation latitude_lastbin latitude_bin0 longitude_bin0
## 1: 0.9446448 -41.60833 0.2771651 0.2771795 -74.99974
## 2: 0.9768031 -41.60883 0.2767428 0.2767589 -74.99945
## 3: 0.9755690 -41.60933 0.2763217 0.2763395 -74.99915
## longitude_lastbin elev_highestreturn elev_lowestmode rh100 pai
## 1: -74.99973 215.9003 193.7073 2219 4.205245
## 2: -74.99944 217.4226 193.6978 2371 2.380235
## 3: -74.99914 221.4247 194.3747 2704 2.826008
## fhd_normal omega pgap_theta cover
## 1: 3.001456 1 0.1213672 0.8760028
## 2: 3.095614 1 0.3031006 0.6948133
## 3: 3.235472 1 0.2423805 0.7553517# check the dimensions
dim(level2BVPM)## [1] 3403 19#filtering for quality_flag == 1
level2BVPM <- level2BVPM[l2b_quality_flag == 1]
dim(level2BVPM)## [1] 3320 19#removing NA values
level2BVPM$pai[level2BVPM$pai==-9999]<-NA # assigning NA to -9999
dim(level2BVPM)## [1] 3320 19#get pai profile metrics
level2BPAIProfile <- getLevel2BPAIProfile(gedilevel2B)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%dim(level2BPAIProfile)## [1] 3403 41#filtering for quality_flag == 1
level2BPAIProfile <- level2BPAIProfile[l2b_quality_flag == 1]
dim(level2BPAIProfile)## [1] 3320 41#convert shot_number from integer to character
level2BPAIProfile$shot_number <- paste0(level2BPAIProfile$shot_number)
#select the previous shots to examine their PAI
smallLevel2BPAIProfile <- level2BPAIProfile[level2BPAIProfile$shot_number %in% sample_shotnums, ]
sample_colors <- c("#FDE725", "#21908D", "#440154")
#checking the data structure
head(smallLevel2BPAIProfile)[, 1:18]## beam shot_number algorithmrun_flag l2b_quality_flag delta_time
## 1: BEAM0101 84280500300242755 1 1 76817130
## 2: BEAM0101 84280500300242757 1 1 76817130
## 3: BEAM0101 84280500300242759 1 1 76817130
## lat_lowestmode lon_lowestmode elev_highestreturn elev_lowestmode
## 1: 0.2362527 -74.99871 238.1628 211.9363
## 2: 0.2354130 -74.99811 233.5918 216.8920
## 3: 0.2345715 -74.99752 222.5892 218.8159
## height_lastbin height_bin0 pai_z0_5m pai_z5_10m pai_z10_15m pai_z15_20m
## 1: -11.120045 26.354006 4.50259972 4.127542 2.777871 1.67275357
## 2: -18.360970 16.714741 1.45967424 1.211199 0.574635 0.01757622
## 3: -6.554091 3.788747 0.07969099 0.000000 0.000000 0.00000000
## pai_z20_25m pai_z25_30m pai_z30_35m
## 1: 0.895469 0.006062514 0
## 2: 0.000000 0.000000000 0
## 3: 0.000000 0.000000000 0#plot some shots
PAIplot <- plotPAIshots(smallLevel2BPAIProfile, col_vector = sample_colors, cutUntil0 = TRUE, idcolumn = 'shot_number')
#plot Level2B PAI Profile
gPAIprofile <- plotPAIProfile(level2BPAIProfile, beam="BEAM0101", elev=TRUE)
#get pavd profile metrics
level2BPAVDProfile <- getLevel2BPAVDProfile(gedilevel2B)##
|
| | 0%
|
|========= | 12%
|
|================== | 25%
|
|========================== | 38%
|
|=================================== | 50%
|
|============================================ | 62%
|
|==================================================== | 75%
|
|============================================================= | 88%
|
|======================================================================| 100%dim(level2BPAVDProfile)## [1] 3403 41#filtering for quality_flag == 1
level2BPAVDProfile <- level2BPAVDProfile[l2b_quality_flag == 1]
dim(level2BPAVDProfile)## [1] 3320 41#convert shot_number from integer to character
level2BPAVDProfile$shot_number <- paste0(level2BPAVDProfile$shot_number)
#select the previous shots to examine their PAI
smallLevel2BPVDProfile <- level2BPAVDProfile[as.character(level2BPAVDProfile$shot_number) %in% sample_shotnums, ]
sample_colors <- c("#FDE725", "#21908D", "#440154")
#checking the data structure
head(smallLevel2BPVDProfile)[, 1:18]## beam shot_number algorithmrun_flag l2b_quality_flag delta_time
## 1: BEAM0101 84280500300242755 1 1 76817130
## 2: BEAM0101 84280500300242757 1 1 76817130
## 3: BEAM0101 84280500300242759 1 1 76817130
## lat_lowestmode lon_lowestmode elev_highestreturn elev_lowestmode
## 1: 0.2362527 -74.99871 238.1628 211.9363
## 2: 0.2354130 -74.99811 233.5918 216.8920
## 3: 0.2345715 -74.99752 222.5892 218.8159
## height_lastbin height_bin0 pavd_z0_5m pavd_z5_10m pavd_z10_15m pavd_z15_20m
## 1: -11.120045 26.354006 0.07501164 0.172472835 0.2454788 0.1882402
## 2: -18.360970 16.714741 0.04969503 0.088503920 0.1193623 0.0574635
## 3: -6.554091 3.788747 0.01593820 0.007969098 0.0000000 0.0000000
## pavd_z20_25m pavd_z25_30m pavd_z30_35m
## 1: 0.166669115 0.0895469 0.0006062514
## 2: 0.001757622 0.0000000 0.0000000000
## 3: 0.000000000 0.0000000 0.0000000000#plot some shots
PAVDplot <- plotPAVDshots(smallLevel2BPVDProfile, col_vector = sample_colors, cutUntil0 = TRUE, idcolumn = 'shot_number')
#plot Level2B PAVD Profile
gPAVDprofile<-plotPAVDProfile(level2BPAVDProfile, beam="BEAM0101", elev=TRUE)
SECTION 4: Create Raster Layers
4.1 Compute Descriptive Statistics
#compute a series of statistics for GEDI relative height metrics
rh100metrics_st <- polyStatsLevel2AM(level2AM_qual, func=SetOfMetrics(rh100))
print(rh100metrics_st) # Check the results## min max mean sd
## 1: 3.77 45.38 31.1662 5.485769#compute a series of statistics for GEDI canopy cover metrics
cover_metrics_st <- polyStatsLevel2BVPM(level2BVPM, func=SetOfMetrics(cover))
print(cover_metrics_st) # Check the results## min max mean sd
## 1: 0.000909579 0.9775591 0.7713407 0.24036964.1 Compute Descriptive Height Statistics
#compute a series of statistics for GEDI relative height metrics
rh100metrics <- gridStatsLevel2AM(level2AM_qual, func=SetOfMetrics(rh100), res=0.010)
rh100metrics # the results here are rasters## class : RasterBrick
## dimensions : 18, 12, 216, 4 (nrow, ncol, ncell, nlayers)
## resolution : 0.01, 0.01 (x, y)
## extent : -74.9999, -74.8799, 0.05794554, 0.2379455 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : min, max, mean, sd
## min values : 3.770000, 27.450001, 13.343846, 2.528121
## max values : 33.889999, 45.380001, 36.730000, 8.796759#plot rh100 stats
rh100maps <- levelplot(rh100metrics,
layout=c(4, 1),
margin=FALSE,
xlab = "Longitude (degree)", ylab = "Latitude (degree)",
colorkey=list(
space='right',
labels=list(at=seq(0, 46, 2), font=4),
axis.line=list(col='black'),
width=1),
par.settings=list(
strip.border=list(col='gray'),
strip.background=list(col='gray'),
axis.line=list(col='gray')
),
scales=list(draw=TRUE),
col.regions=viridis,
at=seq(0, 46, len=101),
names.attr=c("rh100 min","rh100 max","rh100 mean", "rh100 sd"))
rh100maps
#plot rh100 stats over the map
rh100metrics@crs <- CRS("+proj=longlat +datum=WGS84 +no_defs")
sample_level2AM_qual <- level2AM_qual[sample(size = 100, nrow(level2AM_qual)), ]
rhrastvalues <- unique(as.numeric(rh100metrics@data@values))
rhPal<- colorNumeric(palette = "viridis", reverse = FALSE,
domain = rhrastvalues, na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(rh100metrics[['mean']], colors = rhPal, opacity = .7, group = "Mean") %>%
addRasterImage(rh100metrics[['min']], colors = rhPal, opacity = .7, group = "Min") %>%
addRasterImage(rh100metrics[['max']], colors = rhPal, opacity = .7, group = "Max") %>%
addRasterImage(rh100metrics[['sd']], colors = rhPal, opacity = .7, group = "SD") %>%
addLegend(pal = rhPal, values = rhrastvalues,
position = 'topleft',
title= ""#, opacity = .3
#, labFormat = labelFormat(transform = function(x) sort(x, decreasing = TRUE))
) %>%
addCircleMarkers(data = sample_level2AM_qual, lng = sample_level2AM_qual$lon_lowestmode,
lat = sample_level2AM_qual$lat_lowestmode,
popup = paste0('RH100:', round(sample_level2AM_qual$rh100, 2),'<br>',
'Shot:', sample_level2AM_qual$shot_number),
radius = .2, group = 'GEDI',
color = ~rhPal(sample_level2AM_qual$rh100)) %>%
addLayersControl(position = 'topleft',
overlayGroups = c('Mean', 'Min', 'Max', 'SD', 'GEDI'),
options = layersControlOptions(collapsed = FALSE),
baseGroups = c("OpenStreetMap", "Esri.WorldImagery")) %>%
addProviderTiles( "Esri.WorldImagery", group = "Esri.WorldImagery" )## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition4.2 Compute Descriptive Canopy Cover Stats
#compute a series of statistics for GEDI plant area index
pai_metrics <- gridStatsLevel2BVPM(level2BVPM = level2BVPM, func=SetOfMetrics(pai), res=0.010)
fhd_metrics <- gridStatsLevel2BVPM(level2BVPM = level2BVPM, func=SetOfMetrics(fhd_normal), res=0.010)
cover_metrics <- gridStatsLevel2BVPM(level2BVPM = level2BVPM, func=SetOfMetrics(cover), res=0.010)#plot pai stats
pai_maps <- levelplot(pai_metrics,
layout=c(4, 1),
margin=FALSE,
xlab = "Longitude (degree)", ylab = "Latitude (degree)",
colorkey=list(
space='right',
labels=list(at=seq(0, 8, 0.2), font=4),
axis.line=list(col='black'),
width=1),
par.settings=list(
strip.border=list(col='gray'),
strip.background=list(col='gray'),
axis.line=list(col='gray')
),
scales=list(draw=TRUE),
col.regions=viridis,
at=seq(0, 8, len=101),
names.attr=c("PAI min","PAI max","PAI mean", "PAI sd"))
pai_maps
#plot pai stats over the map
pai_metrics@crs <- CRS("+proj=longlat +datum=WGS84 +no_defs")
sample_level2BPAI <- level2BPAIProfile[sample(size = 100, nrow(level2BPAIProfile)), ]
pairastvalues <- unique(as.numeric(pai_metrics@data@values))
paiPal<- colorNumeric(palette = "viridis",
domain = pairastvalues, na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(pai_metrics[['mean']], colors = rhPal, opacity = .7, group = "Mean") %>%
addRasterImage(pai_metrics[['min']], colors = rhPal, opacity = .7, group = "Min") %>%
addRasterImage(pai_metrics[['max']], colors = rhPal, opacity = .7, group = "Max") %>%
addRasterImage(pai_metrics[['sd']], colors = rhPal, opacity = .7, group = "SD") %>%
addLegend(pal = paiPal, values = pairastvalues,
position = 'topleft',
title= ""#, opacity = .
) %>%
addCircleMarkers(data = sample_level2BPAI, lng = sample_level2BPAI$lon_lowestmode,
lat = sample_level2BPAI$lat_lowestmode,
popup = paste0('PAI-5m: ' , round(sample_level2BPAI$pai_z0_5m, 2),'<br>',
'Shot: ', sample_level2BPAI$shot_number),
radius = .2, group = 'GEDI',
color = ~paiPal(sample_level2BPAI$pai_z0_5m)) %>%
addLayersControl(position = 'topleft',
overlayGroups = c('Mean', 'Min', 'Max', 'SD', 'GEDI'),
options = layersControlOptions(collapsed = FALSE),
baseGroups = c("OpenStreetMap", "Esri.WorldImagery")) %>%
addProviderTiles( "Esri.WorldImagery", group = "Esri.WorldImagery" )#plot fhd stats
fhd_maps <- levelplot(fhd_metrics,
layout=c(2, 2),
margin=FALSE,
xlab = "Longitude (degree)", ylab = "Latitude (degree)",
colorkey=list(
space='right',
labels=list(at=seq(0, 3.5, 0.25), font=4),
axis.line=list(col='black'),
width=1),
par.settings=list(
strip.border=list(col='gray'),
strip.background=list(col='gray'),
axis.line=list(col='gray')
),
scales=list(draw=TRUE),
col.regions=viridis,
at=seq(0, 3.5, len=101),
names.attr=c("FHD min","FHD max","FHD mean", "FHD sd"))
fhd_maps
#plot cover stats
cover_maps <- levelplot(cover_metrics,
layout=c(2, 2),
margin=FALSE,
xlab = "Longitude (degree)", ylab = "Latitude (degree)",
colorkey=list(
space='right',
labels=list(at=seq(0, 1, 0.05), font=4),
axis.line=list(col='black'),
width=1),
par.settings=list(
strip.border=list(col='gray'),
strip.background=list(col='gray'),
axis.line=list(col='gray')
),
scales=list(draw=TRUE),
col.regions=viridis,
at=seq(0, 1, len=101),
names.attr=c("CCF min","CCF max","CCF mean", "CCF sd"))
cover_maps
SECTION 5: Compare two vegetation covers
5.1 Define polygon areas
We are definig two contrasting polygons, humand grasslands and forests. For this you can read polygons in WKT format (text) or from vecor/raster files using functions like raster::raster(“mytif.tif”) or rgdal::readOGR(‘mygpkg.gpkg’)
# Provide study areas
wkt_deforested <- "POLYGON((-75.001718%200.23935,-74.998442%200.238726,-74.991579%200.237478,-74.986899%200.236073,-74.981439%200.234513,-74.983311%200.231079,-74.980659%200.22749,-74.975355%200.22234,-74.966463%200.225149,-74.965371%200.229831,-74.969115%200.234513,-74.977695%200.241223,-74.977383%200.244188,-74.969895%200.240287,-74.967243%200.23701,-74.958976%200.234357,-74.958352%200.228114,-74.962564%200.22156,-74.966463%200.220936,-74.967555%200.217971,-74.966775%200.214538,-74.966307%200.213913,-74.96116%200.20892,-74.955076%200.209232,-74.957884%200.204862,-74.96038%200.200961,-74.962252%200.195343,-74.966151%200.196591,-74.967711%200.201117,-74.967087%200.205799,-74.965527%200.209856,-74.970675%200.217035,-74.973764%200.219715,-74.977586%200.221665,-74.982343%200.22346,-74.983045%200.226113,-74.986165%200.229,-74.991157%200.231341,-74.989861%200.230941,-74.99232%200.231292,-74.99388%200.231994,-74.994582%200.234608,-74.995791%200.236481,-74.998794%200.237066,-75.000822%200.237027,-75.001718%200.23935))"
wkt_forest <- "POLYGON((-74.927009%200.104278,-74.892711%200.149222,-74.830976%200.101157,-74.874627%200.049971,-74.927009%200.104278))"
# Convert WKT into R polygons
forest_poly <- readWKT(gsub('%20', ' ', wkt_forest))
deforested_poly <- readWKT(gsub('%20', ' ', wkt_deforested))
## You can load your own polygons loading shapefiles:
# forest_poly <- rgdal::readOGR('/path/to', 'myforestlayer.shp')# Filter coordinates
deforested_pts <- level2AM_qual[na.omit(raster::extract(deforested_poly, level2AM_qual[, c('lon_lowestmode', 'lat_lowestmode')]))$point.ID,]
forest_pts <- level2AM_qual[na.omit(raster::extract(forest_poly, level2AM_qual[, c('lon_lowestmode', 'lat_lowestmode')]))$point.ID,]
ptsPal<- colorNumeric(palette = "viridis", reverse = TRUE,
domain = c(deforested_pts$rh100, forest_pts$rh100), na.color = "transparent")
# Mapping the areas and GEDI points
leaflet() %>%
addCircleMarkers(data = deforested_pts,
lng = deforested_pts$lon_lowestmode,
lat = deforested_pts$lat_lowestmode,
radius = 1, opacity = 0.5,
color = ~ptsPal(rh100),
popup = paste0("shot_number: ", deforested_pts$shot_number)) %>%
addCircleMarkers(data = forest_pts,
forest_pts$lon_lowestmode,
forest_pts$lat_lowestmode,
radius = 1, opacity = 0.5,
color = ~ptsPal(rh100),
popup = paste0("shot_number: ", forest_pts$shot_number)) %>%
addPolygons(data = deforested_poly, color = 'red', fill = NA, group = 'Deforested' ) %>%
addPolygons(data = forest_poly, color = 'blue', fill = NA, group = 'Forest' ) %>%
addLegend(pal = paiPal, values = pairastvalues,
position = 'topleft', title= "Height (m)") %>%
addScaleBar(options = list(imperial = FALSE)) %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addLegend(colors = c("blue", "red"), labels = c("Forest", "Grasslands"))Explore how the energy curves looks in both polygons
rh100Pts <- rbind.data.frame(cbind(deforested_pts, label = 'Deforested', clr = "#F8766D"),
cbind(forest_pts, label = 'Forest', clr = '#00BFC4'))
plotRh100 <- plot_rh_profiles(rh100Pts, idcolumn = 'shot_number', col_vector = rh100Pts$clr, bunch = TRUE)
Now we will compare the profiles and max. height (rh100) in both polygons
ggp <- ggplot(data = rh100Pts[, c('shot_number', 'rh100', 'label', 'clr')],
aes(x = label, y = rh100, color = label)) +
geom_boxplot(fill=0.8, outlier.size=0.5, color = 'black') +
geom_point(position=position_jitter(width=0.3), alpha=0.8) +
guides(colour = guide_legend(override.aes = list(alpha=1))) +
labs(x = "Area", y = "Heigh (m)t", title = 'Max. heigth', color = 'Area')
ggp
SECTION 6: Model applications
Species Distribution Model fake species (created artificial presence/absence points in QGIS)
sp_locations <- readOGR(dsn= "data/fake_sp.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/megreich/Documents/Emma/NAU/GEDI_workshop/workshop/data/fake_sp.shp", layer: "fake_sp"
## with 40 features
## It has 2 fields
## Integer64 fields read as strings: idhead(sp_locations)## id Presence
## 1 1 1
## 2 2 1
## 3 3 1
## 4 4 1
## 5 5 1
## 6 6 1plot(sp_locations[sp_locations$Presence == 0,],col='red',pch=16)
plot(sp_locations[sp_locations$Presence == 1,],col='blue',pch=16)
Add environmental variable raster brick and then add GEDI rasters to it
bioclim_data <- brick(paste0(wd,"/data/bioclim_data_AOI2.tif"))
#change names of layers
names(bioclim_data) <- c("MeanAnnualTemp", "TempAnnualRange", "MeanAnnualPrecip")
#take a look
plot(bioclim_data$MeanAnnualTemp)
plot(bioclim_data$TempAnnualRange)
plot(bioclim_data$MeanAnnualPrecip)
max_canopyht <- rh100metrics$max
#look at extents
bioclim_data## class : RasterBrick
## dimensions : 6, 4, 24, 3 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -75, -74.83333, 0.04166667, 0.2916667 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : bioclim_data_AOI2.tif
## names : MeanAnnualTemp, TempAnnualRange, MeanAnnualPrecip
## min values : 258, 111, 2784
## max values : 260, 113, 2902rh100metrics## class : RasterBrick
## dimensions : 18, 12, 216, 4 (nrow, ncol, ncell, nlayers)
## resolution : 0.01, 0.01 (x, y)
## extent : -74.9999, -74.8799, 0.05794554, 0.2379455 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : min, max, mean, sd
## min values : 3.770000, 27.450001, 13.343846, 2.528121
## max values : 33.889999, 45.380001, 36.730000, 8.796759#resample so that we can combine rasters
#rh100_resample<-resample(rh100metrics$max, bioclim_data, method="ngb")
bioclim_resample<-resample(bioclim_data, max_canopyht, method="ngb")
#add the GEDI max canopy height layer to the Bioclim raster stack
preds <- addLayer(bioclim_resample, max_canopyht)
plot(preds$max)
Using the sdm package to run a Linear model and Random Forest model
#get data into format for sdm package and add 60 random background points
d<- sdmData(Presence~ ., sp_locations, predictors=preds) ## Loading required package: gbm## Warning: package 'gbm' was built under R version 4.0.2## Loaded gbm 2.1.8## Loading required package: tree## Warning: package 'tree' was built under R version 4.0.2## Loading required package: mda## Warning: package 'mda' was built under R version 4.0.2## Loading required package: class## Loaded mda 0.5-2## Loading required package: mgcv## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:raster':
##
## getData## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.## Loading required package: glmnet## Warning: package 'glmnet' was built under R version 4.0.2## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 4.1-1## Loading required package: earth## Warning: package 'earth' was built under R version 4.0.2## Loading required package: Formula## Warning: package 'Formula' was built under R version 4.0.2## Loading required package: plotmo## Warning: package 'plotmo' was built under R version 4.0.2## Loading required package: plotrix## Warning: package 'plotrix' was built under R version 4.0.2## Loading required package: TeachingDemos## Warning: package 'TeachingDemos' was built under R version 4.0.2## Loading required package: rJava## Warning: package 'rJava' was built under R version 4.0.2## Loading required package: RSNNS## Warning: package 'RSNNS' was built under R version 4.0.2## Loading required package: Rcpp## Loading required package: ranger##
## Attaching package: 'ranger'## The following object is masked from 'package:randomForest':
##
## importance## Loading required package: rpart## Loading required package: kernlab## Warning: package 'kernlab' was built under R version 4.0.2##
## Attaching package: 'kernlab'## The following object is masked from 'package:ggplot2':
##
## alpha## The following objects are masked from 'package:raster':
##
## buffer, rotated#what models are available in the sdm package?
getmethodNames()## $bioclim
## [1] "bioclim" "Bioclim"
##
## $bioclim.dismo
## [1] "bioclim.dismo" "BioclimDismo" "bioclimDismo" "bioclim.dis"
## [5] "bioclimD"
##
## $brt
## [1] "brt" "BRT" "gbm" "GBM"
##
## $cart
## [1] "cart" "CART" "tree"
##
## $domain.dismo
## [1] "domain.dismo" "DomainDismo" "domainDismo" "domain.dis" "domainD"
##
## $fda
## [1] "fda" "FDA"
##
## $gam
## [1] "gam" "GAM"
##
## $glm
## [1] "glm" "GLM" "lm"
##
## $glmnet
## [1] "glmnet" "GLMNET" "glmelastic" "glmlasso"
##
## $glmpoly
## [1] "glmpoly" "glmpolynomial" "glmp"
##
## $mahal.dismo
## [1] "mahal.dismo" "MahalDismo" "mahalanoubisDismo"
## [4] "mahal.dis" "mahalD"
##
## $mars
## [1] "mars" "MARS" "earth"
##
## $maxent
## [1] "maxent" "Maxent" "MAXENT" "entropy" "maxentropy"
##
## $maxlike
## [1] "maxlike" "MaxLike"
##
## $mda
## [1] "mda" "MDA"
##
## $mlp
## [1] "mlp" "MLP" "nnet.mlp" "nnetMLP"
##
## $ranger
## [1] "ranger" "rangerRF" "rangerForest"
##
## $rbf
## [1] "rbf" "RBF" "nnet.rbf" "nnetRbf"
##
## $rf
## [1] "rf" "RF" "randomForest" "rforest"
##
## $rpart
## [1] "rpart" "RPART"
##
## $svm
## [1] "svm" "SVM" "ksvm"#Let's try 2: a generalized linear model and a random forest model
m1 <- sdm(Presence~.,data=d,methods=c('glm', 'rf'))## Loading required package: parallelm1## class : sdmModels
## ========================================================
## number of species : 1
## number of modelling methods : 2
## names of modelling methods : glm, rf
## ------------------------------------------
## model run success percentage (per species) :
## ------------------------------------------
## method Presence
## ----------------------
## glm : 100 %
## rf : 100 %
##
## ###################################################################
## model performance (per species), using training test dataset:
## -------------------------------------------------------------------------------
##
## ## species : Presence
## =========================
##
## methods : AUC | COR | TSS | Deviance
## -------------------------------------------------------------------------
## glm : 0.98 | 0.85 | 0.9 | 0.38
## rf : 0.98 | 0.9 | 0.9 | 0.46#evaluate with ROC curves
roc(m1)
#variable importance Plot
getVarImp(m1)##
## The variable importance for all the models are combined (averaged)...## Relative Variable Importance List
## =============================================================
## method : Permutation based on two metrics (Pearson Correlation and AUC)
## number of variables : 4
## variable names : MeanAnnualTemp, TempAnnualRange, MeanAnnualPrecip, max, ...
## number of models : 2
## =============================================================
## Summary of relative variable importance
## ----------------------------------------------
## Based on Correlation metric:
## ----------------------------------------------
## MeanAnnualTemp : [] (3.5 %)
## TempAnnualRange : ***[***--] (14 %)
## MeanAnnualPrecip : *******************[**-] (44.2 %)
## max : ] (2 %)
## =============================================================
## Based on AUC metric:
## ----------------------------------------------
## MeanAnnualTemp [: ***---] (6.1 %)
## TempAnnualRange : [******-----] (13.6 %)
## MeanAnnualPrecip : ***************[**-] (35.9 %)
## max [*] (1.3 %)
## =============================================================#create outputs folder for prediction files
if(!file.exists("outputs")) {dir.create("outputs")}
#predict
p1 <- predict(m1,newdata=preds, filename="./outputs/sdm_prediction", overwrite=T) ## Warning in .rasterFromRasterFile(grdfile, band = band, objecttype, ...): size of
## values file does not match the number of cells (given the data type)
## Warning in .rasterFromRasterFile(grdfile, band = band, objecttype, ...): size of
## values file does not match the number of cells (given the data type)plot(p1)