## ---- include = FALSE---------------------------------------------------------

knitr::opts_chunk$set(warning = FALSE, 
                      message = FALSE,
                      fig.pos = "!H",
                      echo = TRUE,
                      eval = FALSE)          # this vignette by default does not evaluate the code snippets 


## ---- eval = TRUE-------------------------------------------------------------

require(IceSat2R)
require(magrittr)
require(sf)

sf::sf_use_s2(use_s2 = FALSE)                        # disable 's2' in this vignette
mapview::mapviewOptions(leafletHeight = '600px', 
                        leafletWidth = '700px')      # applies to all leaflet maps


## -----------------------------------------------------------------------------
#  
#  #....................
#  # winter (2020, 2021)
#  #....................
#  
#  start_date_w = "2020-12-15"
#  end_date_w = "2021-02-15"
#  
#  rgt_winter = time_specific_orbits(date_from = start_date_w,
#                                    date_to = end_date_w,
#                                    RGT_cycle = NULL,
#                                    download_method = 'curl',
#                                    threads = parallel::detectCores(),
#                                    verbose = TRUE)
#  
#  # ICESAT-2 orbits: 'Earliest-Date' is '2018-10-13'  'Latest-Date' is '2022-06-21'
#  # -----------------------------------------------------
#  # The .zip file of 'RGT_cycle_9' will be downloaded ...
#  # -----------------------------------------------------
#  #   % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
#  #                                  Dload  Upload   Total   Spent    Left  Speed
#  # 100  140M  100  140M    0     0  3662k      0  0:00:39  0:00:39 --:--:-- 2904k
#  # The downloaded .zip file will be extracted in the '/tmp/RtmpPleeLI/RGT_cycle_9' directory .
#  # Download and unzip the RGT_cycle_9 .zip file: Elapsed time: 0 hours and 0 minutes and 41 sec
#  # 138 .kml files will be processed ...
#  # Parallel processing of 138 .kml files using  8  threads starts ...
#  # The 'description' column of the output data will be processed ...
#  # The temproary files will be removed ...
#  # Processing of cycle 'RGT_cycle_9': Elapsed time: 0 hours and 1 minutes and 9 seconds.
#  # -----------------------------------------------------
#  # The .zip file of 'RGT_cycle_10' will be downloaded ...
#  # -----------------------------------------------------
#  #   % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
#  #                                  Dload  Upload   Total   Spent    Left  Speed
#  # 100  140M  100  140M    0     0  3795k      0  0:00:37  0:00:37 --:--:-- 3841k
#  # The downloaded .zip file will be extracted in the '/tmp/RtmpPleeLI/RGT_cycle_10' directory .
#  # Download and unzip the RGT_cycle_10 .zip file: Elapsed time: 0 hours and 0 minutes and 40 sec
#  # 824 .kml files will be processed ...
#  # Parallel processing of 824 .kml files using  8  threads starts ...
#  # The 'description' column of the output data will be processed ...
#  # The temproary files will be removed ...
#  # Processing of cycle 'RGT_cycle_10': Elapsed time: 0 hours and 7 minutes and 53 seconds.
#  # Total Elapsed time: 0 hours and 10 minutes and 25 seconds.
#  

## -----------------------------------------------------------------------------
#  
#  rgt_winter
#  
#  # Simple feature collection with 91390 features and 14 fields
#  # Geometry type: POINT
#  # Dimension:     XY
#  # Bounding box:  xmin: -179.9029 ymin: -87.66478 xmax: 179.9718 ymax: 87.32805
#  # CRS:           4326
#  # First 10 features:
#  #    Name timestamp begin  end  altitudeMode tessellate extrude visibility drawOrder icon  RGT
#  # 1            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 2            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 3            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 4            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 5            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 6            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 7            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 8            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 9            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  # 10           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1250
#  

## -----------------------------------------------------------------------------
#  
#  #..............
#  # summer (2021)
#  #..............
#  
#  start_date_s = "2021-06-15"
#  end_date_s = "2021-08-15"
#  
#  rgt_summer = time_specific_orbits(date_from = start_date_s,
#                                    date_to = end_date_s,
#                                    RGT_cycle = NULL,
#                                    download_method = 'curl',
#                                    threads = parallel::detectCores(),
#                                    verbose = TRUE)
#  
#  # ICESAT-2 orbits: 'Earliest-Date' is '2018-10-13'  'Latest-Date' is '2022-06-21'
#  # -----------------------------------------------------
#  # The .zip file of 'RGT_cycle_11' will be downloaded ...
#  # -----------------------------------------------------
#  #   % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
#  #                                  Dload  Upload   Total   Spent    Left  Speed
#  # 100  140M  100  140M    0     0  3440k      0  0:00:41  0:00:41 --:--:-- 4227k
#  # The downloaded .zip file will be extracted in the '/tmp/RtmpPleeLI/RGT_cycle_11' directory .
#  # Download and unzip the RGT_cycle_11 .zip file: Elapsed time: 0 hours and 0 minutes and 44 sec
#  # 142 .kml files will be processed ...
#  # Parallel processing of 142 .kml files using  8  threads starts ...
#  # The 'description' column of the output data will be processed ...
#  # The temproary files will be removed ...
#  # Processing of cycle 'RGT_cycle_11': Elapsed time: 0 hours and 1 minutes and 30 seconds.
#  # -----------------------------------------------------
#  # The .zip file of 'RGT_cycle_12' will be downloaded ...
#  # -----------------------------------------------------
#  #   % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
#  #                                  Dload  Upload   Total   Spent    Left  Speed
#  # 100  140M  100  140M    0     0  3204k      0  0:00:44  0:00:44 --:--:-- 2333k
#  # The downloaded .zip file will be extracted in the '/tmp/RtmpPleeLI/RGT_cycle_12' directory .
#  # Download and unzip the RGT_cycle_12 .zip file: Elapsed time: 0 hours and 0 minutes and 47 sec
#  # 820 .kml files will be processed ...
#  # Parallel processing of 820 .kml files using  8  threads starts ...
#  # The 'description' column of the output data will be processed ...
#  # The temproary files will be removed ...
#  # Processing of cycle 'RGT_cycle_12': Elapsed time: 0 hours and 8 minutes and 20 seconds.
#  # Total Elapsed time: 0 hours and 11 minutes and 22 seconds.
#  

## -----------------------------------------------------------------------------
#  
#  rgt_summer
#  
#  # Simple feature collection with 89965 features and 14 fields
#  # Geometry type: POINT
#  # Dimension:     XY
#  # Bounding box:  xmin: -179.9861 ymin: -87.66478 xmax: 179.9393 ymax: 87.32805
#  # CRS:           4326
#  # First 10 features:
#  #    Name timestamp begin  end  altitudeMode tessellate extrude visibility drawOrder icon  RGT
#  # 1            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 2            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 3            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 4            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 5            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 6            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 7            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 8            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 9            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  # 10           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1256
#  

## -----------------------------------------------------------------------------
#  
#  #.......................................
#  # extract the 'Greenland Ice Sheet'
#  # glacier from the 'Natural Earth' data
#  #.......................................
#  
#  ne_glaciers = system.file('data_files', 'ne_10m_glaciated_areas.RDS', package = "IceSat2R")
#  ne_obj = readRDS(file = ne_glaciers)
#  
#  greenl_sh = subset(ne_obj, name == "Greenland Ice Sheet")
#  greenl_sh
#  

## -----------------------------------------------------------------------------
#  
#  greenl_sh_east = greenl_sh[2, ]
#  # mapview::mapview(greenl_sh_east, legend = F)
#  
#  #.....................................................
#  # create the bounding of the selected area because
#  # it's required for the 'OpenAltimetry' functions this
#  #  will increase the size of the initial east area
#  #.....................................................
#  
#  bbx_greenl_sh_east = sf::st_bbox(obj = greenl_sh_east)
#  sfc_bbx_greenl_sh_east = sf::st_as_sfc(bbx_greenl_sh_east)
#  # mapview::mapview(sfc_bbx_greenl_sh_east, legend = F)
#  

## -----------------------------------------------------------------------------
#  
#  #..............................................
#  # intersection with the computed "winter" RGT's
#  #..............................................
#  
#  inters_winter = sf::st_intersects(x = sfc_bbx_greenl_sh_east,
#                                    y = sf::st_geometry(rgt_winter),
#                                    sparse = TRUE)
#  
#  #.....................
#  # matched (RGT) tracks
#  #.....................
#  
#  df_inters_winter = data.frame(inters_winter)
#  rgt_subs_winter = rgt_winter[df_inters_winter$col.id, , drop = FALSE]
#  rgt_subs_winter
#  
#  # Simple feature collection with 1079 features and 14 fields
#  # Geometry type: POINT
#  # Dimension:     XY
#  # Bounding box:  xmin: -53.06014 ymin: 61.26295 xmax: -19.23143 ymax: 80.40264
#  # CRS:           4326
#  # First 10 features:
#  #     Name timestamp begin  end  altitudeMode tessellate extrude visibility drawOrder icon  RGT
#  # 117           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1251
#  # 207           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1252
#  # 208           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1252
#  # 209           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1252
#  # 210           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1252
#  # 211           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1252
#  # 212           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1252
#  # 977           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 978           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 979           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  

## -----------------------------------------------------------------------------
#  
#  #..............................................
#  # intersection with the computed "summer" RGT's
#  #..............................................
#  
#  inters_summer = sf::st_intersects(x = sfc_bbx_greenl_sh_east,
#                                    y = sf::st_geometry(rgt_summer),
#                                    sparse = TRUE)
#  
#  #.....................
#  # matched (RGT) tracks
#  #.....................
#  
#  df_inters_summer = data.frame(inters_summer)
#  rgt_subs_summer = rgt_summer[df_inters_summer$col.id, , drop = FALSE]
#  rgt_subs_summer
#  
#  # Simple feature collection with 1066 features and 14 fields
#  # Geometry type: POINT
#  # Dimension:     XY
#  # Bounding box:  xmin: -53.06014 ymin: 61.26295 xmax: -19.23143 ymax: 80.40264
#  # CRS:           4326
#  # First 10 features:
#  #      Name timestamp begin  end  altitudeMode tessellate extrude visibility drawOrder icon  RGT
#  # 407            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 408            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 409            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 410            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 411            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 412            <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1260
#  # 1062           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1267
#  # 1063           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1267
#  # 1064           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1267
#  # 1065           <NA>  <NA> <NA> clampToGround         -1       0          1        NA <NA> 1267
#  

## -----------------------------------------------------------------------------
#  
#  #...............................................
#  # compute the unique RGT's for summer and winter
#  #...............................................
#  
#  unq_rgt_winter = unique(rgt_subs_winter$RGT)
#  unq_rgt_summer = unique(rgt_subs_summer$RGT)
#  dif_rgt = setdiff(unique(rgt_subs_winter$RGT), unique(rgt_subs_summer$RGT))
#  
#  cat(glue::glue("Number of RGT winter: {length(unq_rgt_winter)}"), '\n')
#  # Number of RGT winter: 230
#  
#  cat(glue::glue("Number of RGT summer: {length(unq_rgt_summer)}"), '\n')
#  # Number of RGT summer: 227
#  
#  cat(glue::glue("Difference in RGT: {length(dif_rgt)}"), '\n')
#  # Difference in RGT: 3
#  
#  
#  #......................
#  # find the intersection
#  #......................
#  
#  inters_rgts = intersect(unq_rgt_winter, unq_rgt_summer)
#  
#  
#  #...................................................
#  # create the subset data RGT's for summer and winter
#  #...................................................
#  
#  idx_w = which(rgt_subs_winter$RGT %in% inters_rgts)
#  subs_rgt_winter = rgt_subs_winter[idx_w, , drop = F]
#  
#  idx_s = which(rgt_subs_summer$RGT %in% inters_rgts)
#  subs_rgt_summer = rgt_subs_summer[idx_s, , drop = F]
#  

## -----------------------------------------------------------------------------
#  
#  #...................................................
#  # function to find which of the time specific RGT's
#  # match the OpenAltimetry Track-ID's
#  #...................................................
#  
#  verify_RGTs = function(rgt_subs, bbx_aoi, verbose = FALSE) {
#  
#    tracks_dates_RGT = list()
#  
#    for (item in 1:nrow(rgt_subs)) {
#  
#      dat_item = rgt_subs[item, , drop = F]
#      Date = as.Date(dat_item$Date_time)
#  
#      op_tra = getTracks(minx = as.numeric(bbx_aoi['xmin']),
#                         miny = as.numeric(bbx_aoi['ymin']),
#                         maxx = as.numeric(bbx_aoi['xmax']),
#                         maxy = as.numeric(bbx_aoi['ymax']),
#                         date = as.character(Date),
#                         outputFormat = 'csv',
#                         download_method = 'curl',
#                         verbose = FALSE)
#  
#      date_obj = dat_item$Date_time
#      tim_rgt = glue::glue("Date: {date_obj} Time specific RGT: '{dat_item$RGT}'")
#  
#      if (nrow(op_tra) > 0) {
#  
#        inters_trc = intersect(op_tra$track, dat_item$RGT)
#  
#        if (length(inters_trc) > 0) {
#  
#          iter_op_trac = paste(op_tra$track, collapse = ', ')
#          if (verbose) cat(glue::glue("Row: {item} {tim_rgt}  OpenAltimetry: '{iter_op_trac}'"), '\n')
#  
#          num_subl = glue::glue("{as.character(Date)}_{inters_trc}")
#          tracks_dates_RGT[[num_subl]] = data.table::setDT(list(date = as.character(Date),
#                                                                RGT = inters_trc))
#        }
#      }
#      else {
#        if (verbose) cat(glue::glue("{tim_rgt} without an OpenAltimetry match!"), '\n')
#      }
#    }
#  
#    tracks_dates_RGT = data.table::rbindlist(tracks_dates_RGT)
#    return(tracks_dates_RGT)
#  }
#  

## -----------------------------------------------------------------------------
#  
#  #.........................................................
#  # we keep the relevant columns and remove duplicated
#  # Dates and RGTs to iterate over each pair of observations
#  #.........................................................
#  
#  #.......
#  # winter
#  #.......
#  
#  subs_rgt_w_trc = subs_rgt_winter[, c('RGT', 'Date_time')]
#  subs_rgt_w_trc$Date_time = as.character(as.Date(subs_rgt_w_trc$Date_time))
#  dups_w = which(duplicated(sf::st_drop_geometry(subs_rgt_w_trc)))
#  subs_rgt_w_trc = subs_rgt_w_trc[-dups_w, ]
#  
#  ver_trc_winter = verify_RGTs(rgt_subs = subs_rgt_w_trc,
#                               bbx_aoi = bbx_greenl_sh_east,
#                               verbose = TRUE)
#  
#  colnames(ver_trc_winter) = c('date_winter', 'RGT')
#  
#  #.......
#  # summer
#  #.......
#  
#  subs_rgt_s_trc = subs_rgt_summer[, c('RGT', 'Date_time')]
#  subs_rgt_s_trc$Date_time = as.character(as.Date(subs_rgt_s_trc$Date_time))
#  dups_s = which(duplicated(sf::st_drop_geometry(subs_rgt_s_trc)))
#  subs_rgt_s_trc = subs_rgt_s_trc[-dups_s, ]
#  
#  ver_trc_summer = verify_RGTs(rgt_subs = subs_rgt_s_trc,
#                               bbx_aoi = bbx_greenl_sh_east,
#                               verbose = TRUE)
#  
#  colnames(ver_trc_summer) = c('date_summer', 'RGT')
#  
#  #.............
#  # merge by RGT
#  #.............
#  
#  rgts_ws = merge(ver_trc_winter, ver_trc_summer, by = 'RGT')
#  rgts_ws
#  
#  #       RGT date_summer date_winter
#  #   1:    2  2021-06-23  2020-12-24
#  #   2:   10  2021-06-24  2020-12-24
#  #   3:   11  2021-06-24  2020-12-24
#  #   4:   17  2021-06-24  2020-12-25
#  #   5:   18  2021-06-24  2020-12-25
#  #  ---
#  # 214: 1366  2021-06-22  2020-12-22
#  # 215: 1367  2021-06-22  2020-12-22
#  # 216: 1373  2021-06-22  2020-12-23
#  # 217: 1374  2021-06-22  2020-12-23
#  # 218: 1382  2021-06-23  2020-12-23
#  

## ---- echo = F----------------------------------------------------------------
#  
#  #.........................................................................................
#  # For reproducibility purposes save only the 'date' and 'RGT' matches of winter and summer
#  #.........................................................................................
#  
#  data.table::fwrite(x = rgts_ws, file = system.file('data_files', 'vignette_data', 'winter_summer_inters.csv', package = "IceSat2R"), row.names = F)
#  

## -----------------------------------------------------------------------------
#  
#  greenl_grid = degrees_to_global_grid(minx = as.numeric(bbx_greenl_sh_east['xmin']),
#                                       maxx = as.numeric(bbx_greenl_sh_east['xmax']),
#                                       maxy = as.numeric(bbx_greenl_sh_east['ymax']),
#                                       miny = as.numeric(bbx_greenl_sh_east['ymin']),
#                                       degrees = 4.5,
#                                       square_geoms = TRUE,
#                                       crs_value = 4326,
#                                       verbose = TRUE)
#  

## -----------------------------------------------------------------------------
#  
#  inters_init = sf::st_intersects(sf::st_geometry(greenl_sh_east), greenl_grid)
#  inters_init = data.frame(inters_init)
#  inters_init = inters_init$col.id
#  
#  greenl_grid_subs = greenl_grid[inters_init, , drop = F]
#  
#  # Simple feature collection with 27 features and 1 field
#  # Geometry type: POLYGON
#  # Dimension:     XY
#  # Bounding box:  xmin: -53.10888 ymin: 60.11961 xmax: -17.10888 ymax: 82.61961
#  # CRS:           EPSG:4326
#  # First 10 features:
#  #                          geometry             area
#  # 1  POLYGON ((-53.10888 60.1196... 116709.30 [km^2]
#  # 2  POLYGON ((-48.60888 60.1196... 116709.30 [km^2]
#  # 3  POLYGON ((-44.10888 60.1196... 116709.30 [km^2]
#  # 9  POLYGON ((-53.10888 64.6196...  98928.06 [km^2]
#  # 10 POLYGON ((-48.60888 64.6196...  98928.06 [km^2]
#  # 11 POLYGON ((-44.10888 64.6196...  98928.06 [km^2]
#  # 12 POLYGON ((-39.60888 64.6196...  98928.06 [km^2]
#  # 13 POLYGON ((-35.10888 64.6196...  98928.06 [km^2]
#  # 14 POLYGON ((-30.60888 64.6196...  98928.06 [km^2]
#  # 15 POLYGON ((-26.10888 64.6196...  98928.06 [km^2]
#  

## -----------------------------------------------------------------------------
#  
#  #............
#  # winter join
#  #............
#  
#  subs_join_w = sf::st_join(x = greenl_grid_subs,
#                            y = subs_rgt_w_trc,
#                            join = sf::st_intersects,
#                            left = FALSE)
#  
#  subs_join_w = sf::st_as_sfc(unique(sf::st_geometry(subs_join_w)), crs = 4326)
#  
#  subs_join_w
#  
#  # Geometry set for 9 features
#  # Geometry type: POLYGON
#  # Dimension:     XY
#  # Bounding box:  xmin: -53.10888 ymin: 60.11961 xmax: -17.10888 ymax: 82.61961
#  # CRS:           EPSG:4326
#  # First 5 geometries:
#  # POLYGON ((-53.10888 60.11961, -48.60888 60.1196...
#  # POLYGON ((-48.60888 60.11961, -44.10888 60.1196...
#  # POLYGON ((-44.10888 60.11961, -39.60888 60.1196...
#  # POLYGON ((-26.10888 73.61961, -21.60888 73.6196...
#  # POLYGON ((-21.60888 73.61961, -17.10888 73.6196...
#  

## -----------------------------------------------------------------------------
#  #............
#  # summer join
#  #............
#  
#  subs_join_s = sf::st_join(x = greenl_grid_subs,
#                            y = subs_rgt_s_trc,
#                            join = sf::st_intersects,
#                            left = FALSE)
#  
#  subs_join_s = sf::st_as_sfc(unique(sf::st_geometry(subs_join_s)), crs = 4326)
#  
#  subs_join_s
#  
#  # Geometry set for 9 features
#  # Geometry type: POLYGON
#  # Dimension:     XY
#  # Bounding box:  xmin: -53.10888 ymin: 60.11961 xmax: -17.10888 ymax: 82.61961
#  # CRS:           EPSG:4326
#  # First 5 geometries:
#  # POLYGON ((-53.10888 60.11961, -48.60888 60.1196...
#  # POLYGON ((-48.60888 60.11961, -44.10888 60.1196...
#  # POLYGON ((-44.10888 60.11961, -39.60888 60.1196...
#  # POLYGON ((-26.10888 73.61961, -21.60888 73.6196...
#  # POLYGON ((-21.60888 73.61961, -17.10888 73.6196...
#  

## ---- eval = T, echo = FALSE--------------------------------------------------

subs_join_s = readRDS(file = system.file('data_files', 'vignette_data', 'Greenland_grid_cells.RDS', package = "IceSat2R"))


## -----------------------------------------------------------------------------
#  
#  identical(subs_join_w, subs_join_s)
#  # [1] TRUE
#  

## ---- eval = TRUE-------------------------------------------------------------

mapview::mapview(subs_join_s, legend = F)


## ---- eval = TRUE-------------------------------------------------------------

join_geoms = 1:5
subs_join_reduced = subs_join_s[join_geoms]

mapview::mapview(subs_join_reduced, legend = F)


## ---- eval = T, echo = F, results = 'asis'------------------------------------

pth_subs_RGTs = system.file('data_files', 'vignette_data', 'Greenland_RGTs_grid_cells.csv', package = "IceSat2R")
dtbl_RGTs = data.table::fread(pth_subs_RGTs, header = T, stringsAsFactors = F)
dtbl_RGTs = dtbl_RGTs[order(dtbl_RGTs$RGT, decreasing = F), ]

stargazer::stargazer(dtbl_RGTs,
                     summary = FALSE, 
                     rownames = FALSE, 
                     header = FALSE, 
                     float = FALSE,
                     table.placement = 'h', 
                     title = 'Greenland Grid Cells')


# #.................................................................
# # I came to these RGT's after processing the next for-loop for
# #   - all 9 Greenland geoms 
# #   - all RGT's (of the 'rgts_ws' variable)
# #.................................................................
# 
# require(magrittr)
# logs_out_dtbl = data.table::rbindlist(logs_out)         # the output of the 'logs_out' list is required!
# logs_out_dtbl$index = names(dat_out_w)                  # the output of the 'dat_out_w' list is required!
# logs_out_dtbl = logs_out_dtbl[order(logs_out_dtbl$N_rows_winter, decreasing = T), ]
# logs_out_dtbl_subs = subset(logs_out_dtbl, N_rows_winter > 0 & N_rows_summer > 0)
# vec_grid = as.vector(unlist(lapply(strsplit(logs_out_dtbl_subs$index, '_'), function(x) x[3])))
# logs_out_dtbl_subs$greenland_cell = vec_grid
# spl_rgt = split(logs_out_dtbl_subs, by = 'greenland_cell')
# spl_subs = spl_rgt[as.character(1:5)] %>%
#   lapply(function(x) x[1:4, ]) %>%
#   data.table::rbindlist()


## -----------------------------------------------------------------------------
#  
#  #...............................................
#  # keep a subset of RGTs and Greenland Grid cells
#  #...............................................
#  
#  RGTs = c(33, 41, 56, 94, 108, 284,
#           421, 422, 437, 460, 475,
#           521, 544, 658, 681, 787,
#           794, 1290, 1366, 1373)
#  
#  #......................
#  # update the input data
#  #......................
#  
#  rgts_ws_reduced = subset(rgts_ws, RGT %in% RGTs)
#  rgts_ws_RGT = rgts_ws_reduced$RGT
#  
#  
#  #..........................................................
#  # we'll loop over the RGT's for the Date start and Date end
#  #..........................................................
#  
#  # length(unique(rgts_ws_reduced$RGT)) == nrow(rgts_ws_reduced)   # check for duplicated RGT's
#  rgts_ws_reduced$date_summer = as.Date(rgts_ws_reduced$date_summer)
#  rgts_ws_reduced$date_winter = as.Date(rgts_ws_reduced$date_winter)
#  min_max_dates = apply(rgts_ws_reduced[, -1], 2, function(x) c(min(x), max(x)))
#  
#  start_w = as.character(min_max_dates[1, 'date_winter'])
#  end_w = as.character(min_max_dates[2, 'date_winter'])
#  
#  start_s = as.character(min_max_dates[1, 'date_summer'])
#  end_s = as.character(min_max_dates[2, 'date_summer'])
#  
#  dat_out_w = dat_out_s = logs_out = list()
#  LEN = length(subs_join_reduced)
#  LEN_rgt = length(rgts_ws_RGT)
#  
#  t_start = proc.time()
#  for (idx_grid in 1:LEN) {
#  
#    geom_iter = subs_join_reduced[idx_grid]
#    bbx_iter = sf::st_bbox(obj = geom_iter, crs = 4326)
#  
#    for (j in 1:LEN_rgt) {
#  
#      message("Greenland Geom: ",
#              idx_grid, "/",
#              LEN, "   RGT-index: ",
#              j, "/",
#              LEN_rgt, "\r",
#              appendLF = FALSE)
#  
#      utils::flush.console()
#  
#      track_i = rgts_ws_RGT[j]
#      name_iter = glue::glue("geom_idx_{idx_grid}_RGT_{track_i}")
#  
#      iter_dat_winter = get_level3a_data(minx = as.numeric(bbx_iter['xmin']),
#                                         miny = as.numeric(bbx_iter['ymin']),
#                                         maxx = as.numeric(bbx_iter['xmax']),
#                                         maxy = as.numeric(bbx_iter['ymax']),
#                                         startDate = start_w,
#                                         endDate = end_w,
#                                         trackId = track_i,
#                                         beamName = NULL,        # return data of all 6 beams
#                                         product = 'atl06',
#                                         client = 'portal',
#                                         outputFormat = 'csv',
#                                         verbose = FALSE)
#  
#      iter_dat_summer = get_level3a_data(minx = as.numeric(bbx_iter['xmin']),
#                                         miny = as.numeric(bbx_iter['ymin']),
#                                         maxx = as.numeric(bbx_iter['xmax']),
#                                         maxy = as.numeric(bbx_iter['ymax']),
#                                         startDate = start_s,
#                                         endDate = end_s,
#                                         trackId = track_i,
#                                         beamName = NULL,        # return data of all 6 beams
#                                         product = 'atl06',
#                                         client = 'portal',
#                                         outputFormat = 'csv',
#                                         verbose = FALSE)
#  
#      NROW_w = nrow(iter_dat_winter)
#      NROW_s = nrow(iter_dat_summer)
#  
#      if (NROW_s > 0 & NROW_w > 0) {
#        iter_logs = list(RGT = track_i,
#                         N_rows_winter = NROW_w,
#                         N_rows_summer = NROW_s)
#  
#        logs_out[[name_iter]] = data.table::setDT(iter_logs)
#        dat_out_w[[name_iter]] = iter_dat_winter
#        dat_out_s[[name_iter]] = iter_dat_summer
#      }
#    }
#  }
#  IceSat2R:::compute_elapsed_time(time_start = t_start)
#  # Elapsed time: 0 hours and 6 minutes and 15 seconds.

## -----------------------------------------------------------------------------
#  
#  logs_out_dtbl = data.table::rbindlist(logs_out)
#  logs_out_dtbl$index = names(dat_out_w)
#  

## ---- eval = TRUE, echo = FALSE-----------------------------------------------

pth_logs = system.file('data_files', 'vignette_data', 'LOGs.csv', package = "IceSat2R")
logs_out_dtbl = data.table::fread(file = pth_logs, stringsAsFactors = F, header = T)


## ---- eval = TRUE, results = 'asis'-------------------------------------------

logs_out_dtbl = logs_out_dtbl[order(logs_out_dtbl$N_rows_winter, decreasing = T), ]

stargazer::stargazer(logs_out_dtbl, 
                     summary = FALSE, 
                     rownames = FALSE, 
                     header = FALSE, 
                     float = FALSE,
                     table.placement = 'h', 
                     title = 'LOGs')


## -----------------------------------------------------------------------------
#  
#  #................................
#  # we pick one with approx. same
#  # rows for both summer and winter
#  #................................
#  
#  # names(dat_out_w)
#  Greenland_Geom_index = 4
#  RGT = 1290
#  
#  sublist_name = glue::glue("geom_idx_{Greenland_Geom_index}_RGT_{RGT}")
#  
#  #...............
#  # winter sublist
#  #...............
#  
#  w_subs = dat_out_w[[sublist_name]]
#  
#  w_subs
#  
#  
#  #               date segment_id longitude latitude      h_li atl06_quality_summary track_id beam
#  #      1: 2020-12-17     566420 -22.43278 78.11961  553.5237                     0     1290 gt1l
#  #      2: 2020-12-17     566421 -22.43294 78.11943  553.5601                     0     1290 gt1l
#  #      3: 2020-12-17     566422 -22.43310 78.11926  553.5809                     0     1290 gt1l
#  #      4: 2020-12-17     566423 -22.43326 78.11908  553.5331                     0     1290 gt1l
#  #      5: 2020-12-17     566424 -22.43341 78.11891  553.5400                     0     1290 gt1l
#  #     ---
#  # 152558: 2020-12-17     591901 -25.70488 73.62039 1640.8195                     0     1290 gt3r
#  # 152559: 2020-12-17     591902 -25.70496 73.62021 1650.2201                     0     1290 gt3r
#  # 152560: 2020-12-17     591903 -25.70505 73.62003 1657.0470                     0     1290 gt3r
#  # 152561: 2020-12-17     591904 -25.70514 73.61985 1665.8877                     0     1290 gt3r
#  # 152562: 2020-12-17     591905 -25.70523 73.61968 1675.1165                     0     1290 gt3r
#  

## -----------------------------------------------------------------------------
#  
#  #...............
#  # summer sublist
#  #...............
#  
#  s_subs = dat_out_s[[sublist_name]]
#  
#  s_subs
#  
#  #               date segment_id longitude latitude      h_li atl06_quality_summary track_id beam
#  #      1: 2021-06-17     566420 -22.43270 78.11961  554.1851                     0     1290 gt1l
#  #      2: 2021-06-17     566421 -22.43286 78.11943  554.2042                     0     1290 gt1l
#  #      3: 2021-06-17     566422 -22.43302 78.11925  554.1635                     0     1290 gt1l
#  #      4: 2021-06-17     566423 -22.43318 78.11908  554.1546                     0     1290 gt1l
#  #      5: 2021-06-17     566424 -22.43334 78.11890  554.2347                     0     1290 gt1l
#  #     ---
#  # 152048: 2021-06-17     591901 -25.70511 73.62039 1638.9667                     0     1290 gt3r
#  # 152049: 2021-06-17     591902 -25.70521 73.62022 1646.6494                     0     1290 gt3r
#  # 152050: 2021-06-17     591903 -25.70530 73.62004 1654.2173                     0     1290 gt3r
#  # 152051: 2021-06-17     591904 -25.70539 73.61986 1663.3259                     0     1290 gt3r
#  # 152052: 2021-06-17     591905 -25.70548 73.61968 1673.2574                     0     1290 gt3r
#  

## -----------------------------------------------------------------------------
#  
#  cols_keep = c('date', 'segment_id', 'longitude', 'latitude', 'h_li', 'beam')
#  
#  w_subs_hq = subset(w_subs, atl06_quality_summary == 0)
#  w_subs_hq = w_subs_hq[, ..cols_keep]
#  colnames(w_subs_hq) = glue::glue("{cols_keep}_winter")
#  
#  s_subs_hq = subset(s_subs, atl06_quality_summary == 0)
#  s_subs_hq = s_subs_hq[, ..cols_keep]
#  colnames(s_subs_hq) = glue::glue("{cols_keep}_summer")
#  
#  sw_hq_merg = merge(x = w_subs_hq,
#                     y = s_subs_hq,
#                     by.x = c('segment_id_winter', 'beam_winter'),
#                     by.y = c('segment_id_summer', 'beam_summer'))
#  
#  sw_hq_merg$dif_height = sw_hq_merg$h_li_winter - sw_hq_merg$h_li_summer
#  
#  summary(sw_hq_merg$dif_height)
#  
#  #     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
#  # -42.6563  -0.4222  -0.0854  -0.1567   0.0537  26.4173
#  

## -----------------------------------------------------------------------------
#  
#  cols_viz = c('segment_id_winter', 'beam_winter', 'h_li_winter', 'h_li_summer')
#  ws_vis = sw_hq_merg[, ..cols_viz]

## ---- eval = TRUE, echo = FALSE-----------------------------------------------

ws_vis = readRDS(file = system.file('data_files', 'vignette_data', 'plotly_beams.RDS', package = "IceSat2R"))


## ---- eval = TRUE-------------------------------------------------------------

ws_vis_mlt = reshape2::melt(ws_vis, id.vars = c('segment_id_winter', 'beam_winter'))
ws_vis_mlt = data.table::data.table(ws_vis_mlt, stringsAsFactors = F)
ws_vis_mlt_spl = split(ws_vis_mlt, by = 'beam_winter')
# BEAMS = names(ws_vis_mlt_spl)       # plot all beams


## ---- eval = T----------------------------------------------------------------

#...................................
# function to plot each subplot beam
#...................................

plotly_beams = function(spl_data, 
                        beam,
                        left_width,
                        left_height,
                        right_width,
                        right_height) {
  
  subs_iter = spl_data[[beam]]
  
  cat(glue::glue("Plot for Beam '{beam}' will be created ..."), '\n')
  
  #......................
  # plot for all segments
  #......................

  fig_lines = plotly::plot_ly(data = subs_iter,
                              x = ~segment_id_winter,
                              y = ~value,
                              color = ~variable,
                              colors = c("blue", "red"), 
                              line = list(width = 2),
                              text = ~glue::glue("land-ice-height: {value}  Segment-id: {segment_id_winter}"),
                              hoverinfo = "text",
                              width = left_width,
                              height = left_height) %>%
    
    plotly::layout(xaxis = list(gridcolor = "grey", showgrid = T),
                   yaxis = list(gridcolor = "grey", showgrid = T)) %>%
    
    plotly::add_lines()
    
  #..............................
  # plot for a subset of segments
  #..............................
  
  segm_ids = 588326:588908     # this subset of segments show a big difference betw. summer and winter
  subs_iter_segm = subset(subs_iter, segment_id_winter %in% segm_ids)
  
  fig_lines_splt = plotly::plot_ly(data = subs_iter_segm,
                                   x = ~segment_id_winter,
                                   y = ~value,
                                   color = ~variable,
                                   colors = c("blue", "red"), 
                                   line = list(width = 2),
                                   text = ~glue::glue("land-ice-height: {value}  Segment-id: {segment_id_winter}"),
                                   hoverinfo = "text",
                                   width = right_width,
                                   height = right_height) %>%
    
    plotly::layout(xaxis = list(gridcolor = "grey", showgrid = T),
                   yaxis = list(gridcolor = "grey", showgrid = T)) %>%
    
    plotly::add_lines(showlegend = FALSE)
  
  both_plt = plotly::subplot(list(fig_lines, fig_lines_splt), 
                             nrows = 1, 
                             margin = 0.03, 
                             widths = c(0.7, 0.3)) %>% 
    plotly::layout(title = glue::glue("Beam: '{beam}' ( Segments: from {min(segm_ids)} to {max(segm_ids)} )"))
  # plotly::export(p = both_plt, file = glue::glue('{beam}.png'))
  
  return(both_plt)
}


## ---- eval = TRUE-------------------------------------------------------------

plt_gt1l = plotly_beams(spl_data = ws_vis_mlt_spl, 
                        beam = "gt1l",
                        left_width = 1800,
                        left_height = 800,
                        right_width = 900,
                        right_height = 400)

plt_gt1l


## ---- eval = TRUE-------------------------------------------------------------

plt_gt1r = plotly_beams(spl_data = ws_vis_mlt_spl, 
                        beam = "gt1r",
                        left_width = 1800,
                        left_height = 800,
                        right_width = 900,
                        right_height = 400)
plt_gt1r


## ---- echo = F----------------------------------------------------------------
#  
#  #......................................................
#  # save all images for all beams in a separate directory
#  #......................................................
#  
#  nams_ws = names(dat_out_w)
#  save_summary = save_dat = list()
#  
#  for (nam_iter in nams_ws) {
#  
#    cat("-----------------\n")
#    cat(nam_iter, '\n')
#    cat("-----------------\n")
#  
#    w_subs = dat_out_w[[nam_iter]]
#    s_subs = dat_out_s[[nam_iter]]
#  
#    cols_keep = c('date', 'segment_id', 'longitude', 'latitude', 'h_li', 'beam')
#  
#    w_subs_hq = subset(w_subs, atl06_quality_summary == 0)
#    w_subs_hq = w_subs_hq[, ..cols_keep]
#    colnames(w_subs_hq) = glue::glue("{cols_keep}_winter")
#  
#    s_subs_hq = subset(s_subs, atl06_quality_summary == 0)
#    s_subs_hq = s_subs_hq[, ..cols_keep]
#    colnames(s_subs_hq) = glue::glue("{cols_keep}_summer")
#  
#    sw_hq_merg = merge(x = w_subs_hq,
#                       y = s_subs_hq,
#                       by.x = c('segment_id_winter', 'beam_winter'),
#                       by.y = c('segment_id_summer', 'beam_summer'))
#  
#    if (nrow(sw_hq_merg) > 0) {
#      sw_hq_merg$dif_height = sw_hq_merg$h_li_winter - sw_hq_merg$h_li_summer
#  
#      save_dat[[nam_iter]] = sw_hq_merg
#      save_summary[[nam_iter]] = data.table::setDT(list(name_iter = nam_iter,
#                                                        min = min(sw_hq_merg$dif_height),
#                                                        mean = mean(sw_hq_merg$dif_height),
#                                                        median = median(sw_hq_merg$dif_height),
#                                                        max = max(sw_hq_merg$dif_height),
#                                                        N_rows = nrow(sw_hq_merg)))
#  
#      #.......................................
#      # save the plots for visual verification
#      #.......................................
#  
#      cols_viz = c('segment_id_winter', 'beam_winter', 'h_li_winter', 'h_li_summer')
#      ws_vis = sw_hq_merg[, ..cols_viz]
#  
#      ws_vis_mlt = reshape2::melt(ws_vis, id.vars = c('segment_id_winter', 'beam_winter'))
#      ws_vis_mlt = data.table::data.table(ws_vis_mlt, stringsAsFactors = F)
#      ws_vis_mlt_spl = split(ws_vis_mlt, by = 'beam_winter')
#  
#      dir_save = file.path('all_beams_all_RGTs', nam_iter)           # !! create the 'all_beams_all_RGTs' directory first
#      if (!dir.exists(dir_save)) dir.create(dir_save)
#  
#      BEAMS = names(ws_vis_mlt_spl)       # plot all beams
#  
#      for (beam in BEAMS) {
#  
#        subs_iter = ws_vis_mlt_spl[[beam]]
#  
#        cat(glue::glue("Plot for Beam '{beam}' will be saved ..."), '\n')
#  
#        #......................
#        # plot for all segments
#        #......................
#  
#        fig_lines = plotly::plot_ly(data = subs_iter,
#                                    x = ~segment_id_winter,
#                                    y = ~value,
#                                    color = ~variable,
#                                    colors = c("blue", "red"),
#                                    line = list(width = 2),
#                                    text = ~glue::glue("land-ice-height: {value}  Segment-id: {segment_id_winter}"),
#                                    hoverinfo = "text",
#                                    width = 1800,
#                                    height = 1000) %>%
#  
#          plotly::layout(xaxis = list(gridcolor = "grey", showgrid = T),
#                         yaxis = list(gridcolor = "grey", showgrid = T)) %>%
#  
#          plotly::add_lines()
#  
#        plotly::export(p = fig_lines, file = file.path(dir_save, glue::glue('{beam}.png')))
#      }
#    }
#    else {
#     message(glue::glue("Empty data table after merging for idx and RGT: '{nam_iter}'"))
#    }
#  }
#  
#  save_summary = data.table::rbindlist(save_summary)
#  save_summary = save_summary[order(save_summary$max, decreasing = T), ]
#  save_summary
#  

## -----------------------------------------------------------------------------
#  
#  Greenland_Geom_index = 2
#  RGT = 33
#  
#  sublist_name = glue::glue("geom_idx_{Greenland_Geom_index}_RGT_{RGT}")
#  
#  w_subs = dat_out_w[[sublist_name]]
#  s_subs = dat_out_s[[sublist_name]]
#  
#  cols_keep = c('date', 'segment_id', 'longitude', 'latitude', 'h_li', 'beam')
#  
#  w_subs_hq = subset(w_subs, atl06_quality_summary == 0)
#  w_subs_hq = w_subs_hq[, ..cols_keep]
#  colnames(w_subs_hq) = glue::glue("{cols_keep}_winter")
#  
#  s_subs_hq = subset(s_subs, atl06_quality_summary == 0)
#  s_subs_hq = s_subs_hq[, ..cols_keep]
#  colnames(s_subs_hq) = glue::glue("{cols_keep}_summer")
#  
#  sw_hq_merg = merge(x = w_subs_hq,
#                     y = s_subs_hq,
#                     by.x = c('segment_id_winter', 'beam_winter'),
#                     by.y = c('segment_id_summer', 'beam_summer'))
#  

## -----------------------------------------------------------------------------
#  
#  #...............................
#  # compute the pair-wise distance
#  #...............................
#  
#  sw_hq_merg$dist_dif = geodist::geodist(x = sw_hq_merg[, c('longitude_winter', 'latitude_winter')],
#                                         y = sw_hq_merg[, c('longitude_summer', 'latitude_summer')],
#                                         paired = TRUE,
#                                         measure = 'geodesic')
#  #..............
#  # split by beam
#  #..............
#  
#  spl_beam = split(sw_hq_merg, by = 'beam_winter')
#  
#  #...............................................
#  # compute the summary of the distance to observe
#  # the beams with the lowest distance difference
#  #...............................................
#  
#  sm_stats = lapply(spl_beam, function(x) {
#    summary(x$dist_dif)
#  })
#  
#  # $gt1l
#  #     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
#  # 0.006829 1.568755 2.107793 2.083547 2.654482 3.911961
#  # $gt1r
#  #      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
#  # 0.0008427 0.2483260 0.5219864 0.5766648 0.8521821 2.0041663
#  # $gt2l
#  #      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
#  # 0.0002195 0.3001418 0.6112493 0.6753225 0.9860298 2.2545978
#  # $gt2r
#  #    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  #   3.170   4.697   5.356   5.300   5.931   7.170
#  # $gt3l
#  #    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  #   1.197   2.679   3.115   3.155   3.621   5.069
#  # $gt3r
#  #    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#  #   3.881   5.066   5.488   5.498   5.949   7.114
#  

## -----------------------------------------------------------------------------
#  
#  #......................................
#  # keep only the 'gt1r' and 'gt2l' beams
#  #......................................
#  
#  sw_hq_merg_beams = subset(sw_hq_merg, beam_winter %in% c('gt1r', 'gt2l'))
#  
#  
#  #.......................................................
#  # keep only a pair of coordinates and rename the columns
#  #.......................................................
#  
#  sw_hq_merg_beams = sw_hq_merg_beams[, c('segment_id_winter', 'beam_winter', 'longitude_winter',
#                                          'latitude_winter', 'h_li_winter', 'h_li_summer')]
#  colnames(sw_hq_merg_beams) = c('segment_id', 'beam', 'longitude',
#                                 'latitude', 'h_li_winter', 'h_li_summer')
#  
#  #       segment_id beam longitude latitude h_li_winter h_li_summer
#  #    1:     351810 gt1r -44.10891 63.26390    2724.035    2725.023
#  #    2:     351811 gt1r -44.10895 63.26408    2724.004    2725.155
#  #    3:     351812 gt1r -44.10900 63.26426    2724.445    2725.298
#  #    4:     351813 gt1r -44.10904 63.26444    2724.530    2725.461
#  #    5:     351814 gt1r -44.10908 63.26462    2725.084    2725.650
#  #   ---
#  # 8700:     359399 gt1r -44.44637 64.61873    2777.249    2777.390
#  # 8701:     359400 gt1r -44.44642 64.61891    2777.136    2777.297
#  # 8702:     359401 gt1r -44.44647 64.61908    2777.081    2777.192
#  # 8703:     359402 gt1r -44.44651 64.61926    2776.936    2777.076
#  # 8704:     359403 gt1r -44.44656 64.61944    2776.818    2776.968
#  

## -----------------------------------------------------------------------------
#  
#  sf_aoi = sf::st_as_sf(sw_hq_merg_beams, coords = c('longitude', 'latitude'), crs = 4326)
#  bbx_aoi = sf::st_bbox(sf_aoi)
#  sfc_aoi = sf::st_as_sfc(bbx_aoi)
#  
#  dem_dir = tempdir()
#  dem_dir
#  
#  dem30 = CopernicusDEM::aoi_geom_save_tif_matches(sf_or_file = sfc_aoi,
#                                                   dir_save_tifs = dem_dir,
#                                                   resolution = 30,
#                                                   crs_value = 4326,
#                                                   threads = parallel::detectCores(),
#                                                   verbose = TRUE)
#  
#  if (nrow(dem30$csv_aoi) > 1) {          # create a .VRT file if I have more than 1 .tif files
#    file_out = file.path(dem_dir, 'VRT_mosaic_FILE.vrt')
#  
#    vrt_dem30 = CopernicusDEM::create_VRT_from_dir(dir_tifs = dem_dir,
#                                                   output_path_VRT = file_out,
#                                                   verbose = TRUE)
#  }
#  
#  if (nrow(dem30$csv_aoi) == 1) {         # if I have a single .tif file keep the first index
#    file_out = list.files(dem_dir, pattern = '.tif', full.names = T)[1]
#  }
#  
#  
#  #.......................................................
#  # crop the raster to the bounding box of the coordinates
#  #.......................................................
#  
#  rst_inp = terra::rast(x = file_out)
#  vec_crop = terra::vect(x = sfc_aoi)
#  rst_crop = terra::crop(x = rst_inp,
#                         y = vec_crop,
#                         snap = "out")      # snap = "in" gives NA's
#  
#  #...............................................
#  # we also have to find the closest elevation
#  # value to the 'winter' and 'summer' coordinates
#  # using the raster resolution
#  #...............................................
#  
#  ter_dtbl = data.table::as.data.table(x = rst_crop, xy = TRUE, cells = TRUE)
#  colnames(ter_dtbl) = c("cell", "lon_dem30", "lat_dem30", "dem30")
#  # length(unique(ter_dtbl$cell)) == nrow(ter_dtbl)
#  
#  xy = as.matrix(sw_hq_merg_beams[, c('longitude', 'latitude')])
#  sw_cells = terra::cellFromXY(object = rst_crop, xy = xy)
#  
#  sw_hq_merg_beams$cell = sw_cells
#  # length(unique(sw_hq_merg_beams$cell)) < nrow(sw_hq_merg_beams)
#  
#  merg_cells = merge(x = sw_hq_merg_beams, y = ter_dtbl, by = 'cell')
#  
#  
#  #......................................................................
#  # compute also the difference in distance between the beam measurements
#  # and the DEM coordinates (based on the 30-meter resolution cells)
#  #......................................................................
#  
#  merg_cells$dem_dif_dist = geodist::geodist(x = merg_cells[, c('longitude', 'latitude')],
#                                             y = merg_cells[, c('lon_dem30', 'lat_dem30')],
#                                             paired = TRUE,
#                                             measure = 'geodesic')
#  
#  summary(merg_cells$dem_dif_dist)
#  
#   #   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#   # 0.2356  8.1619 11.5943 11.1490 14.3201 20.5394
#  

## -----------------------------------------------------------------------------
#  
#  cols_viz_dem = c('beam', 'longitude', 'latitude', 'h_li_winter', 'h_li_summer', 'dem30')
#  merg_cells_viz = merg_cells[, ..cols_viz_dem]
#  

## ---- echo = FALSE, eval = FALSE----------------------------------------------
#  
#  merg_cells_viz = readRDS(file = system.file('data_files', 'vignette_data', '3d_plot.RDS', package = "IceSat2R"))
#  

## ---- eval = FALSE------------------------------------------------------------
#  
#  merg_cells_viz_mlt = reshape2::melt(merg_cells_viz, id.vars = c('beam', 'longitude', 'latitude'))
#  merg_cells_viz_mlt = data.table::data.table(merg_cells_viz_mlt, stringsAsFactors = F)
#  colnames(merg_cells_viz_mlt) = c('beam', 'longitude', 'latitude',  'variable', 'height')
#  
#  fig_height = plotly::plot_ly(merg_cells_viz_mlt,
#                               x = ~longitude,
#                               y = ~latitude,
#                               z = ~height,
#                               split = ~beam,           # split by beam
#                               type = 'scatter3d',
#                               mode = 'lines',
#                               color = ~variable,
#                               line = list(width = 10),
#                               width = 1000,
#                               height = 900)
#  fig_height
#  

