carlos-alberto-silva / rGEDI

rGEDI: An R Package for NASA's Global Ecosystem Dynamics Investigation (GEDI) Data Visualization and Processing.
157 stars 67 forks source link

Possibly incorrect column numbers dz bins in plotPAVDProfile #24

Closed iosefa closed 3 years ago

iosefa commented 3 years ago

Hi,

I encountered unexpected results when I plotted the plans area volume density profile using the plotPAVDProfile function. Looking at the code, it seems like the function is referencing values by column number:

ref https://github.com/carlos-alberto-silva/rGEDI/blob/7b3316bd6d56265cfc83cc8d0b7fcd36319b9194/R/plotPAVDProfile.R#L59

However, these don't look to be the correct columns. As a result, the plotted figure shows pavd values between 0 and what looks like the elevation of the highest return.

The output of getLevel2BPAVDProfile looks like (using data described below):

> str(level2BPAVDProfile)
Classes ‘data.table’ and 'data.frame':  7 obs. of  41 variables:
 $ beam              : chr  "BEAM0010" "BEAM0010" "BEAM0010" "BEAM0010" ...
 $ shot_number       :integer64 77980203000498040 77980203200498041 77980203400498042 77980203600498043 77980203800498044 77980204000498045 77980204200498046 
 $ algorithmrun_flag : int  1 1 1 1 0 1 1
 $ l2b_quality_flag  : int  1 1 1 0 0 1 1
 $ delta_time        : num  73306625 73306625 73306625 73306625 73306625 ...
 $ lat_lowestmode    : num  -16.3 -16.3 -16.3 -16.3 -16.3 ...
 $ lon_lowestmode    : num  46.8 46.8 46.8 46.8 46.8 ...
 $ elev_highestreturn: num  157 158 158 161 162 ...
 $ elev_lowestmode   : num  143 143 145 156 148 ...
 $ height_lastbin    : num  -3.65 -2.64 -1.85 -2.57 -9999 ...
 $ height_bin0       : num  14.3 15.3 14.5 14.7 -9999 ...
 $ pavd_z0_5m        : num  1.87e-01 2.17e-01 2.43e-01 3.08e-03 -1.00e+04 ...
 $ pavd_z5_10m       : num  2.07e-01 2.25e-01 2.59e-01 1.54e-03 -1.00e+04 ...

This should work as a reproducible example:

library("rGEDI")

ymax <- -16.3188638251292
xmin <- 46.805069035532
ymin <- -16.3216145717765
xmax <- 46.8097723784352

daterange <- c("2019-01-01", "2022-01-01")

gedi_list <- gedifinder(
  "GEDI02_B",
  ymax,
  xmin,
  ymin,
  xmax,
  version="001",
  daterange=daterange
)

# here are the results of `gedifinder`:
# > gedi_list
# [1] "https://e4ftl01.cr.usgs.gov/GEDI/GEDI02_B.001/2020.04.28/GEDI02_B_2020119094224_O07798_T02281_02_001_01.h5"
# [2] "https://e4ftl01.cr.usgs.gov/GEDI/GEDI02_B.001/2019.06.09/GEDI02_B_2019160175158_O02778_T02281_02_001_01.h5"

# assuming that I downloaded the first result into the directory `data/rs/GEDI/`:
filepath <- "data/rs/GEDI/GEDI02_B_2020119094224_O07798_T02281_02_001_01.h5"

level2B <-  readLevel2B(level2Bpath=filepath)

level2BPAVDProfile <- getLevel2BPAVDProfile(level2B) %>%
  clipLevel2BPAVDProfile(xmin, xmax, ymin, ymax)

# looking at the data structure (same as above snippet)
str(level2BPAVDProfile)

# plotting the figure
plotPAVDProfile(level2BPAVDProfile, beam = "BEAM0010")

The plot looks like:

Screen Shot 2021-01-25 at 16 35 14

However, if I modify ref https://github.com/carlos-alberto-silva/rGEDI/blob/7b3316bd6d56265cfc83cc8d0b7fcd36319b9194/R/plotPAVDProfile.R#L59 such that it reads:

  dft <- data.table::melt(level2BPAVDProfile_sub[, c(2, 9, 11, 12 : 41)], id.vars = c("shot_number", "elev_lowestmode", "height_bin0"), variable.name = "pavd", value.name = "value")

I get a figure that looks more like I would expect:

Screen Shot 2021-01-25 at 16 45 05
iosefa commented 3 years ago

I forgot to mention 2 other things:

carlos-alberto-silva commented 3 years ago

Fixed

iosefa commented 3 years ago

thank you! :) And thank you for building this fantastic library!