Workshops

AGU 2022 GEDI Workshop

Workshop Icebreaker Activity

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

  1. 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.

  2. 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.

  3. 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:

Launch Rstudio Binder

To download this code, please refer to the workshop github:

https://github.com/IGSA-SICCS/GEDI_workshop

Tutorial:

SECTION 0: PREPARE R ENVIRONMENT

0.1: Libraries

#Define which libraries to use
library(raster)       ## spatial tools
library(devtools)     ## tools for building and instaling 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)      ## nice plot tools
library(RColorBrewer) ## color palettes
library(rasterVis)    ## raster viz tools
library(viridis)      ## color palettes

0.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 5: DOWNLOAD GEDI DATA FOR A REGION OF INTEREST

5.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")

5.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")

5.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."

5.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."

5.5: Use gedi_download to download GEDI data (DO NOT RUN – takes a while and we have some files downloades already)

#download data product to their output data directory
#This functions will download heavy files (>1GB each) that can take long time
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 2.1 Read GEDI data (DO NOT RUN)

#read the level 1B, 2A, and 2B for a single orbit
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)

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"))

###START 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

#read GEDI L2A data into data.table
level2AM <- getLevel2AM(gedilevel2A)
##
  |
  |                                                                      |   0%
  |
  |=========                                                             |  12%
  |
  |==================                                                    |  25%
  |
  |==========================                                            |  38%
  |
  |===================================                                   |  50%
  |
  |============================================                          |  62%
  |
  |====================================================                  |  75%
  |
  |=============================================================         |  88%
  |
  |======================================================================| 100%
head(level2AM, 3) ## Check the data structure
##        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
##    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
##    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
##     rh10 rh11 rh12 rh13 rh14 rh15 rh16 rh17 rh18 rh19 rh20 rh21 rh22 rh23 rh24
## 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 9.19 9.56 9.90
## 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 1.90 2.09 2.27
## 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 5.11 5.86 6.53
##     rh25  rh26  rh27  rh28  rh29  rh30  rh31  rh32  rh33  rh34  rh35  rh36
## 1: 10.23 10.57 10.87 11.13 11.39 11.65 11.91 12.14 12.32 12.55 12.74 12.92
## 2:  2.50  2.69  2.87  3.10  3.36  3.58  3.88  4.14  4.48  4.85  5.23  5.64
## 3:  7.09  7.54  7.92  8.33  8.70  9.11  9.45  9.78 10.12 10.46 10.87 11.32
##     rh37  rh38  rh39  rh40  rh41  rh42  rh43  rh44  rh45  rh46  rh47  rh48
## 1: 13.11 13.30 13.45 13.63 13.78 13.93 14.08 14.23 14.38 14.53 14.68 14.83
## 2:  6.05  6.42  6.83  7.24  7.69  8.14  8.59  9.00  9.34  9.63  9.90 10.16
## 3: 11.76 12.17 12.59 12.96 13.30 13.59 13.89 14.19 14.45 14.72 15.01 15.28
##     rh49  rh50  rh51  rh52  rh53  rh54  rh55  rh56  rh57  rh58  rh59  rh60
## 1: 14.98 15.13 15.24 15.39 15.54 15.69 15.84 15.95 16.10 16.25 16.36 16.51
## 2: 10.38 10.61 10.83 11.05 11.32 11.58 11.91 12.25 12.66 13.00 13.33 13.63
## 3: 15.54 15.80 16.06 16.32 16.58 16.88 17.14 17.37 17.59 17.78 17.97 18.19
##     rh61  rh62  rh63  rh64  rh65  rh66  rh67  rh68  rh69  rh70  rh71  rh72
## 1: 16.62 16.77 16.88 16.99 17.11 17.22 17.33 17.44 17.56 17.67 17.74 17.85
## 2: 13.89 14.16 14.38 14.60 14.83 15.01 15.20 15.43 15.61 15.80 15.99 16.21
## 3: 18.38 18.56 18.79 19.01 19.20 19.42 19.65 19.83 20.02 20.21 20.39 20.54
##     rh73  rh74  rh75  rh76  rh77  rh78  rh79  rh80  rh81  rh82  rh83  rh84
## 1: 17.93 18.04 18.15 18.23 18.34 18.41 18.53 18.60 18.71 18.83 18.90 19.01
## 2: 16.40 16.58 16.77 16.96 17.14 17.33 17.52 17.70 17.89 18.08 18.26 18.41
## 3: 20.69 20.84 20.99 21.14 21.29 21.44 21.59 21.74 21.89 22.08 22.23 22.41
##     rh85  rh86  rh87  rh88  rh89  rh90  rh91  rh92  rh93  rh94  rh95  rh96
## 1: 19.12 19.20 19.31 19.42 19.54 19.65 19.80 19.91 20.06 20.21 20.39 20.58
## 2: 18.60 18.75 18.90 19.09 19.24 19.42 19.65 19.83 20.10 20.39 20.81 21.37
## 3: 22.60 22.79 22.97 23.16 23.35 23.53 23.76 23.98 24.21 24.43 24.69 24.95
##     rh97  rh98  rh99 rh100
## 1: 20.84 21.10 21.52 22.19
## 2: 21.93 22.41 22.90 23.72
## 3: 25.25 25.66 26.26 27.05
dim(level2AM)
## [1] 3403  112
#convert shot_number from integer to character
level2AM$shot_number <- paste0(level2AM$shot_number)

#filtering for quality_flag == 1
level2AM_qual <- level2AM[quality_flag == 1]
dim(level2AM_qual)
## [1] 3323  112

2.4 Visualize location of quality filtered GEDI footprints

#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") %>%
  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.4 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", "GEDI Footprints"))

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 raw waveforms
par(mfrow = c(1,1), mar=c(4,4,1,1), cex.axis = 1.5)

plot(wf_forest, relative=FALSE, polygon=TRUE, type="l", lwd=2, col="#FDE725",
    xlab="Waveform Amplitude", ylab="Elevation (m)", main="Raw Waveform: Forest", xlim=c(200,500), ylim=c(200,260))

plot(wf_edge, relative=FALSE, polygon=TRUE, type="l", lwd=2, col="#21908D",
    xlab="Waveform Amplitude", ylab="Elevation (m)", main="Raw Waveform: Edge", xlim=c(200,500), ylim=c(200,260), add = TRUE)

plot(wf_grass, relative=FALSE, polygon=TRUE, type="l", lwd=2, col="#440154",
    xlab="Waveform Amplitude", ylab="Elevation (m)", main="Raw Waveform: Grass", xlim=c(200,500), ylim=c(200,260), add = TRUE)

3.2 Plot GEDIL2A Elevation and Height Metrics

#plot waveform with GEDIL2A relative height metrics
plotWFMetrics(level1b = gedilevel1B, level2a = gedilevel2A, shot_number = shotnum_forest, rh=c(25, 50, 75, 90))
##
  |
  |                                                                      |   0%
  |
  |=========                                                             |  12%
  |
  |==================                                                    |  25%
  |
  |==========================                                            |  38%
  |
  |===================================                                   |  50%
  |
  |============================================                          |  62%
  |
  |====================================================                  |  75%
  |
  |=============================================================         |  88%
  |
  |======================================================================| 100%

plotWFMetrics(gedilevel1B, gedilevel2A, shotnum_edge, rh=c(25, 50, 75, 90))
##
  |
  |                                                                      |   0%
  |
  |=========                                                             |  12%
  |
  |==================                                                    |  25%
  |
  |==========================                                            |  38%
  |
  |===================================                                   |  50%
  |
  |============================================                          |  62%
  |
  |====================================================                  |  75%
  |
  |=============================================================         |  88%
  |
  |======================================================================| 100%

plotWFMetrics(gedilevel1B, gedilevel2A, shotnum_grass, rh=c(25, 50, 75, 90))
##
  |
  |                                                                      |   0%
  |
  |=========                                                             |  12%
  |
  |==================                                                    |  25%
  |
  |==========================                                            |  38%
  |
  |===================================                                   |  50%
  |
  |============================================                          |  62%
  |
  |====================================================                  |  75%
  |
  |=============================================================         |  88%
  |
  |======================================================================| 100%

Lets use the same elevation range to compare the previous shots

#plot waveform with GEDIL2A relative height metrics but with a standard scale

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.43 45.38 27.49649 8.093237
#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.2403696

4.1 Compute Descriptive Height Stats

#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 : 23, 15, 345, 4  (nrow, ncol, ncell, nlayers)
## resolution : 0.01, 0.01  (x, y)
## extent     : -74.9999, -74.8499, 0.04716706, 0.2771671  (xmin, xmax, ymin, ymax)
## crs        : NA
## source     : memory
## names      :        min,        max,       mean,         sd
## min values :  3.4300001, 19.9799995,  8.5405555,  0.4242643
## max values :   32.58000,   45.38000,   35.14636,   11.57263
#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 = TRUE,
                      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 definition
## 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 definition
## 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 definition
## 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 definition
## 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 definition
## 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 definition
## 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 definition
## 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 definition

4.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", reverse = TRUE,
                      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" )
## 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 definition
## 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 definition
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## 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 definition
## 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 definition
## 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 definition
## 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 definition
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## 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 definition
## 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 definition
#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
deforested_poly <- readWKT(gsub('%20', ' ', wkt_deforested))
forest_poly <- readWKT(gsub('%20', ' ', wkt_forest))


# 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() %>%
  addPolygons(data = deforested_poly, color = 'red', fill = NA, group = 'Deforested' ) %>%
  addPolygons(data = forest_poly, color = 'blue', fill = NA, group = 'Forest' ) %>%
  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))  %>%
    addLegend(pal = paiPal, values = pairastvalues,
            position = 'topleft', title= "Height (m)") %>%
  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"))

Now we will compare the profiles and max. height (rh100) in both polygons

rh100Pts <- rbind.data.frame(cbind(deforested_pts, label = 'Deforested', clr = "#F8766D"),
                                   cbind(forest_pts, label = 'Forest', clr = '#B79F00'))

ggplot(data = rh100Pts[, c('shot_number', 'rh100', 'label', 'clr')],
       aes(x = label, y = label, 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 = "Height", title = 'Max. heigth', color = 'Area')

plotRh100 <- plot_rh_profiles(rh100Pts, idcolumn = 'label', col_vector = rh100Pts$clr)