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.


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
## Simple feature collection with 131765 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -179.9986 ymin: -87.66742 xmax: 179.9984 ymax: 87.3305
## Geodetic CRS:  WGS 84
## First 10 features:
##    day_of_year       Date hour minute second RGT                      geometry
## 1          356 2021-12-22    7     57     49   1 POINT (-0.1318472 0.02795893)
## 2          356 2021-12-22    7     58     49   1   POINT (-0.5162124 3.868758)
## 3          356 2021-12-22    7     59     49   1    POINT (-0.901809 7.709809)
## 4          356 2021-12-22    8      0     49   1    POINT (-1.289879 11.55065)
## 5          356 2021-12-22    8      1     49   1    POINT (-1.681755 15.39082)
## 6          356 2021-12-22    8      2     49   1     POINT (-2.078916 19.2299)
## 7          356 2021-12-22    8      3     49   1    POINT (-2.483051 23.06748)
## 8          356 2021-12-22    8      4     49   1    POINT (-2.896146 26.90316)
## 9          356 2021-12-22    8      5     49   1      POINT (-3.3206 30.73662)
## 10         356 2021-12-22    8      6     49   1    POINT (-3.759374 34.56754)


ICESat-2 and Countries intersection


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,


cntr = rnaturalearth::ne_countries(scale = 110, type = 'countries', returnclass = 'sf')
cntr = cntr[, c('sovereignt', 'sov_a3')]
cntr
## Simple feature collection with 177 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## CRS:           +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## First 10 features:
##             sovereignt sov_a3                       geometry
## 0          Afghanistan    AFG MULTIPOLYGON (((61.21082 35...
## 1               Angola    AGO MULTIPOLYGON (((16.32653 -5...
## 2              Albania    ALB MULTIPOLYGON (((20.59025 41...
## 3 United Arab Emirates    ARE MULTIPOLYGON (((51.57952 24...
## 4            Argentina    ARG MULTIPOLYGON (((-65.5 -55.2...
## 5              Armenia    ARM MULTIPOLYGON (((43.58275 41...
## 6           Antarctica    ATA MULTIPOLYGON (((-59.57209 -...
## 7               France    FR1 MULTIPOLYGON (((68.935 -48....
## 8            Australia    AU1 MULTIPOLYGON (((145.398 -40...
## 9              Austria    AUT MULTIPOLYGON (((16.97967 48...


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,


dat_both = suppressMessages(sf::st_join(x = res_rgt_many,
                                        y = cntr, 
                                        join = sf::st_intersects, 
                                        left = TRUE))
dat_both
## Simple feature collection with 131765 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -179.9986 ymin: -87.66742 xmax: 179.9984 ymax: 87.3305
## Geodetic CRS:  WGS 84
## First 10 features:
##    day_of_year       Date hour minute second RGT   sovereignt sov_a3
## 1          356 2021-12-22    7     57     49   1         <NA>   <NA>
## 2          356 2021-12-22    7     58     49   1         <NA>   <NA>
## 3          356 2021-12-22    7     59     49   1        Ghana    GHA
## 4          356 2021-12-22    8      0     49   1 Burkina Faso    BFA
## 5          356 2021-12-22    8      1     49   1         Mali    MLI
## 6          356 2021-12-22    8      2     49   1         Mali    MLI
## 7          356 2021-12-22    8      3     49   1         Mali    MLI
## 8          356 2021-12-22    8      4     49   1      Algeria    DZA
## 9          356 2021-12-22    8      5     49   1      Algeria    DZA
## 10         356 2021-12-22    8      6     49   1      Morocco    MAR
##                         geometry
## 1  POINT (-0.1318472 0.02795893)
## 2    POINT (-0.5162124 3.868758)
## 3     POINT (-0.901809 7.709809)
## 4     POINT (-1.289879 11.55065)
## 5     POINT (-1.681755 15.39082)
## 6      POINT (-2.078916 19.2299)
## 7     POINT (-2.483051 23.06748)
## 8     POINT (-2.896146 26.90316)
## 9       POINT (-3.3206 30.73662)
## 10    POINT (-3.759374 34.56754)


The unique number of RGT’s for “RGT_cycle_14” are


length(unique(dat_both$RGT))
## [1] 1387


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


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), ]


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



all RGT’s (1387 in number) intersect with “Antarctica” and almost all with “Russia”.


‘Onshore’ and ‘Offshore’ Points ICESat-2 coverage


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


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


stargazer::stargazer(dtbl_land_sea,
                     type = 'html',
                     summary = FALSE, 
                     rownames = TRUE, 
                     header = FALSE, 
                     table.placement = 'h', 
                     title = 'Land and Sea Proportions')
Land and Sea Proportions
percentage Num_Icesat2_points
sea 67.070 88,369
land 32.930 43,396


Global glaciated areas and ICESat-2 coverage


We can also observe the ICESat-2 “RGT_cycle_14” coverage based on the 1 to 10 million large scale Natural Earth Glaciated Areas data,


data(ne_10m_glaciated_areas)


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


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
## Simple feature collection with 68 features and 5 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.99993 xmax: 180 ymax: 82.96573
## Geodetic CRS:  WGS 84
## First 10 features:
##     recnum scalerank      featurecla                    name min_zoom
## 143    143         3 Glaciated areas    Mount Brown Icefield      2.1
## 148    148         5 Glaciated areas    Braithwaite Icefield      5.0
## 152    152         3 Glaciated areas         Hooker Icefield      2.1
## 206    206         5 Glaciated areas       Homathko Icefield      5.0
## 214    214         6 Glaciated areas Clachnacudainn Icefield      5.7
## 215    215         6 Glaciated areas         Albert Icefield      5.7
## 228    228         3 Glaciated areas        Plateau Icefield      2.1
## 230    230         5 Glaciated areas      Pemberton Icefield      5.0
## 256    256         3 Glaciated areas        Cambria Icefiled      2.1
## 273      0         3 Glaciated areas          Lyell Icefield      2.1
##                           geometry
## 143 MULTIPOLYGON (((-118.4066 5...
## 148 MULTIPOLYGON (((-119.9303 5...
## 152 MULTIPOLYGON (((-117.8572 5...
## 206 MULTIPOLYGON (((-124.6489 5...
## 214 MULTIPOLYGON (((-118.0284 5...
## 215 MULTIPOLYGON (((-117.6752 5...
## 228 MULTIPOLYGON (((-123.8453 5...
## 230 MULTIPOLYGON (((-123.3869 5...
## 256 MULTIPOLYGON (((-129.661 56...
## 273 MULTIPOLYGON (((-117.2649 5...


and we’ll visualize the subset using the mapview package,


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