Neal Swayze

GEDI Data Processing and Visualization - Post 1

The NASA Global Ecosystem Dynamics Investigation (aka GEDI) is a recent NASA project using high resolution laser ranging observations to model the three dimensional structure of the earth. The GEDI instrument makes precise measurements of forest canopy height, canopy vertical structure, and surface elevation.

The GEDI sensor sends out pulses of laser light, which travel from the sensor to the earths surface, interact with the ground/vegetation, then a certain proportion of the light pulse returns and is captured by the sensors telescope. The time between the emission and reception of the light pulse is then used to calculate a waveform (distribution). These waveforms are further processed to estimate an elevation, ground surface, and a variety of canopy metrics, including relative height, percent cover, foliage height diversity, and plant area index.

The goal of the GEDI project will be to investigate data fusions for “scaling-up” GEDI footprint information to landscape scales using additional continuous remote sensing data sources. The landscape scale structure metrics will then be incorporated into habitat modeling for wildlife species of management and conservation interest, and for identifying diversity “hotspots” of these species.

Part of investigating new data comes with the task of visualizing abstract data. Sure, tables and graphs are typical, but what about something a bit more visually appealing?

Downloading and Processing GEDI DATA

### Install rGEDI package
install.packages("rGEDI")

### Load Nessesary Packages
library(rGEDI)
library(lidR)
library(data.table)
library(dplyr)
library(sp)
library(sf)

### Set the working directory
rootdir <- ("D:/gedi_example")
dir.create(file.path(rootdir, "gedi_downloads"))
download_dir <- file.path(rootdir, "/gedi_downloads")
dir.create(file.path(rootdir, "gedi_processed"))
processed_data_dir <- file.path(rootdir, "/gedi_processed")

### Define Study area boundary box coordinates
ul_lat<- -44.0654
lr_lat<- -44.17246
ul_lon<- -13.76913
lr_lon<- -13.67646

### Specifying the date range
daterange=c("2019-07-01","2019-09-23")

### Find the desired granules
gLevel2A <- gedifinder(product="GEDI02_A",ul_lat, ul_lon, lr_lat, lr_lon,version="002",daterange=daterange)
gLevel2A

### Download the desired granule file
gediDownload(filepath=gLevel2A, outdir=download_dir)


### Generating a list of granule files to process from the GEDI downloads folder
setwd(download_dir)
list <- list.files(setwd(download_dir), pattern = "*.h5")
list

### Define the desired beams you want to retain - (these are full power beams)
target_beams <- c("BEAM0101", "BEAM0110", "BEAM1000", "BEAM1011")

### For loop to read in and process the GEDI data
for (i in seq(list)){

  ### Get the first object in the list
  desired_granule <- list[[i]]
  desired_granule
  
  ### Get an output name
  out_name <- tools::file_path_sans_ext(desired_granule)
  out_name_shp <- paste0(out_name, ".shp")
  out_name_csv <- paste0(out_name, ".csv")
  
  ### Read in the .h5 file, convert to dataframe object
  setwd(download_dir)
  file = readLevel2A(desired_granule)
  file = getLevel2AM(file)
  file = as.data.frame(file)

  ### Filter the GEDI dataframe to contain desired information
  file = subset(file, select = c(shot_number, beam, degrade_flag, quality_flag,
                                 sensitivity, solar_elevation, lat_lowestmode,
                                 lon_lowestmode, elev_lowestmode, rh0, rh5, 
                                 rh10, rh20, rh30, rh40,rh50, rh60, rh65, 
                                 rh70, rh75, rh80, rh85, rh90, rh95, rh96, 
                                 rh97, rh98, rh99, rh100))
  
  file = filter(file, ((solar_elevation < 0) & 
                         (degrade_flag == 0) & 
                         (quality_flag == 1) & 
                         (beam %in% target_beams) & 
                         (sensitivity >= 0.95) ) )
  
  ### Convert the gedi dataframe to a spatial points object, reproject to EPSG 5070
  file_spdf = SpatialPointsDataFrame(cbind(file$lon_lowestmode,file$lat_lowestmode), data=file)
  proj4string(file_spdf) = CRS("+init=epsg:4326")
  file_spdf <- spTransform(file_spdf, CRS("+init=epsg:5070"))
  proj4string(file_spdf) = CRS("+init=epsg:5070")

  ### Write out a shapefile object
  out_shp <- st_as_sf(file_spdf)
  setwd(processed_data_dir)
  st_write(out_shp, out_name_shp, append=FALSE)
  
  ### Convert back to dataframe, write out .csv file
  file <- as.data.frame(file_spdf)
  setwd(processed_data_dir)
  fwrite(file, file = out_name_csv, sep = ",")
  
  ### Clean up your working environment to save memory
  rm(file, out_name, desired_granule)
  gc()
}

#title

Visualizing GEDI Data

#title

title

Generating Point Cloud from GEDI footprints

### Read in the processed CSV file
setwd(processed_data_dir)
gedi_data <- fread("GEDI02_A_2019226215413_O03805_04_T02285_02_003_01_V002.csv")
gedi_data <- as.data.frame(gedi_data)
gedi_data

### Creating a las file for the extracted X, Y, and Z, writing out file
x = round(gedi_data$coords.x1, digits = 3)
y = round(gedi_data$coords.x2, digits = 3)
z = round(gedi_data$elev_lowestmode, digits = 3)

### Write out the las file using the rlas package
lasdata = data.frame(X = x, Y = y, Z = z)
lasheader = header_create(lasdata)
file = file.path(processed_data_dir, "gedi_true_elevation.las")
setwd(processed_data_dir)
write.las(file, lasheader, lasdata)

#title

#title