---
title: "ICESat-2 Mission Orbits"
output: rmarkdown::html_vignette
always_allow_html: true
vignette: >
  %\VignetteIndexEntry{ICESat-2 Mission Orbits}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r, include = FALSE}
knitr::opts_chunk$set(fig.width = 4,
                      fig.height = 4, 
                      fig.align = "center",
                      fig.pos = "!H",
                      warning = FALSE, 
                      message = FALSE,
                      echo = TRUE,
                      eval = TRUE)
```

<br>

This *first* vignette demonstrates how to download and process *time specific orbits*. We'll use one of the *Reference Ground Track (RGT) cycles* and merge it with other data sources with the purpose to visualize specific areas. 

We'll load one of the latest which is *"RGT_cycle_14"* (from *December 22, 2021* to *March 23, 2022*). The documentation of the *"RGT_cycle_14"* data includes more details on how a user can come to the same data format for any of the RGT Cycles.

<br>

```{r}

pkgs = c('IceSat2R', 'magrittr', 'mapview', 'sf', 'rnaturalearth', 
         'data.table', 'DT', 'stargazer')
load_pkgs = lapply(pkgs, require, character.only = TRUE)  # load required R packages

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

#.............................
# load the 'RGT_cycle_14' data
#.............................

data(RGT_cycle_14)

res_rgt_many = sf::st_as_sf(x = RGT_cycle_14, coords = c('longitude', 'latitude'), crs = 4326)
res_rgt_many

```

<br>

## ICESat-2 and Countries intersection

<br>

We'll proceed to merge the orbit geometry points with the countries data of the *rnaturalearth* R package (1:110 million scales) and for this purpose, we keep only the *"sovereignt"* and *"sov_a3"* columns,

<br>

```{r}

cntr = rnaturalearth::ne_countries(scale = 110, type = 'countries', returnclass = 'sf')
cntr = cntr[, c('sovereignt', 'sov_a3')]
cntr

```

<br>

We then merge the orbit points with the country geometries and specify also *"left = TRUE"* to keep also observations that do not intersect with the *rnaturalearth* countries data,

<br>

```{r}

dat_both = suppressMessages(sf::st_join(x = res_rgt_many,
                                        y = cntr, 
                                        join = sf::st_intersects, 
                                        left = TRUE))
dat_both

```

<br>

The unique number of RGT's for *"RGT_cycle_14"* are

<br>

```{r}

length(unique(dat_both$RGT))

```


<br>

We observe that from *December 22, 2021* to *March 23, 2022*,

<br>


```{r}

df_tbl = data.frame(table(dat_both$sovereignt), stringsAsFactors = F)
colnames(df_tbl) = c('country', 'Num_IceSat2_points')

df_subs = dat_both[, c('RGT', 'sovereignt')]
df_subs$geometry = NULL
df_subs = data.table::data.table(df_subs, stringsAsFactors = F)
colnames(df_subs) = c('RGT', 'country')
df_subs = split(df_subs, by = 'country')
df_subs = lapply(df_subs, function(x) {
  unq_rgt = sort(unique(x$RGT))
  items = ifelse(length(unq_rgt) < 5, length(unq_rgt), 5)
  concat = paste(unq_rgt[1:items], collapse = '-')
  iter_dat = data.table::setDT(list(country = unique(x$country), 
                                    Num_RGTs = length(unq_rgt), 
                                    first_5_RGTs = concat))
  iter_dat
})

df_subs = data.table::rbindlist(df_subs)

df_tbl = merge(df_tbl, df_subs, by = 'country')
df_tbl = df_tbl[order(df_tbl$Num_IceSat2_points, decreasing = T), ]

```

<br>

```{r}

DT_dtbl = DT::datatable(df_tbl, rownames = FALSE)

```

<br>

```{r, echo = FALSE}

DT_dtbl

```

<br>

all RGT's (1387 in number) intersect with *"Antarctica"* and almost all with *"Russia"*.

<br>

## 'Onshore' and 'Offshore' Points ICESat-2 coverage

<br>

The **onshore** and **offshore** number of ICESat-2 points and percentages for the *"RGT_cycle_14"* equal to

<br>

```{r} 

num_sea = sum(is.na(dat_both$sovereignt))
num_land = sum(!is.na(dat_both$sovereignt))

perc_sea = round(num_sea / nrow(dat_both), digits = 4) * 100.0
perc_land = round(num_land / nrow(dat_both), digits = 4) * 100.0

dtbl_land_sea = data.frame(list(percentage = c(perc_sea, perc_land),
                                Num_Icesat2_points = c(num_sea, num_land)))

row.names(dtbl_land_sea) = c('sea', 'land')


```

<br>

```{r, results = 'asis'}

stargazer::stargazer(dtbl_land_sea,
                     type = 'html',
                     summary = FALSE, 
                     rownames = TRUE, 
                     header = FALSE, 
                     table.placement = 'h', 
                     title = 'Land and Sea Proportions')

```

<br>

## Global glaciated areas and ICESat-2 coverage

<br>

We can also observe the ICESat-2 *"RGT_cycle_14"* coverage based on the 1 to 10 million large scale [Natural Earth Glaciated Areas](https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-glaciated-areas/) data,

<br>

```{r}

data(ne_10m_glaciated_areas)

```

<br>

We'll restrict the processing to the major polar glaciers (that have a name included),

<br>

```{r}

ne_obj_subs = subset(ne_10m_glaciated_areas, !is.na(name))
ne_obj_subs = sf::st_make_valid(x = ne_obj_subs)      # check validity of geometries
ne_obj_subs

```

<br>

and we'll visualize the subset using the *mapview* package,

<br>

```{r}

mpv = mapview::mapview(ne_obj_subs, 
                       color = 'cyan', 
                       col.regions = 'blue', 
                       alpha.regions = 0.5, 
                       legend = FALSE)
mpv
```

<br>

We will see which orbits of the ICESat-2 *"RGT_cycle_14"* intersect with these major polar glaciers,

<br>

```{r}

res_rgt_many$id_rgt = 1:nrow(res_rgt_many)       # include 'id' for fast subsetting

dat_glac_sf = suppressMessages(sf::st_join(x = ne_obj_subs,
                                           y = res_rgt_many, 
                                           join = sf::st_intersects))

dat_glac = data.table::data.table(sf::st_drop_geometry(dat_glac_sf), stringsAsFactors = F)
dat_glac = dat_glac[complete.cases(dat_glac), ]              # keep non-NA observations
dat_glac

```

<br>

We'll split the merged data by the *'name'* of the glacier,

<br>

```{r}

dat_glac_name = split(x = dat_glac, by = 'name')

sum_stats_glac = lapply(dat_glac_name, function(x) {
  
  dtbl_glac = x[, .(name_glacier = unique(name), 
                    Num_unique_Dates = length(unique(Date)),
                    Num_unique_RGTs = length(unique(RGT)))]
  dtbl_glac
})

sum_stats_glac = data.table::rbindlist(sum_stats_glac)
sum_stats_glac = sum_stats_glac[order(sum_stats_glac$Num_unique_RGTs, decreasing = T), ]

```

<br>

The next table shows the total number of days and RGTs for each one of the major polar glaciers,

<br>

```{r, results = 'asis'}

stargazer::stargazer(sum_stats_glac, 
                     type = 'html',
                     summary = FALSE, 
                     rownames = FALSE, 
                     header = FALSE, 
                     table.placement = 'h', 
                     title = 'Days and RGTs')

```

<br>

We can restrict to one of the glaciers to visualize the ICESat-2 *"RGT_cycle_14"* coverage over this specific area (*'Southern Patagonian Ice Field'*),

<br>

```{r}

sample_glacier = 'Southern Patagonian Ice Field'
dat_glac_smpl = dat_glac_name[[sample_glacier]]

```

<br>

```{r, results = 'asis'}

cols_display = c('name', 'day_of_year', 'Date', 'hour', 'minute', 'second', 'RGT')

stargazer::stargazer(dat_glac_smpl[, ..cols_display],
                     type = 'html',
                     summary = FALSE, 
                     rownames = FALSE, 
                     header = FALSE, 
                     table.placement = 'h', 
                     title = 'Southern Patagonian Ice Field')

```

<br>

and we gather the intersected RGT coordinates points with the selected glacier,

<br>

```{r}

subs_rgts = subset(res_rgt_many, id_rgt %in% dat_glac_smpl$id_rgt)

set.seed(1)
samp_colrs = sample(x = grDevices::colors(distinct = TRUE), 
                    size = nrow(subs_rgts))
subs_rgts$color = samp_colrs

```

<br>

```{r}

ne_obj_subs_smpl = subset(ne_obj_subs, name == sample_glacier)

mpv_glacier = mapview::mapview(ne_obj_subs_smpl, 
                               color = 'cyan', 
                               col.regions = 'blue', 
                               alpha.regions = 0.5, 
                               legend = FALSE)

mpv_RGTs = mapview::mapview(subs_rgts,
                            color = subs_rgts$color,
                            alpha.regions = 0.0,
                            lwd = 6,
                            legend = FALSE)
```

<br>

and visualize both the glacier and the subset of the intersected RGT coordinate points (of the different Days) in the same map. The clickable map and point popups include more information,

<br>


```{r}

lft = mpv_glacier + mpv_RGTs
lft

```

<br>
