Site icon R-bloggers

Greenland ice thickness

Michael
[This article was first published on r.iresmi.net, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Meltwater in crevasses in Greenland – CC-BY-NC by NASA’s Marshall Space Flight Center

Day 11 of 30DayMapChallenge: « Arctic » (previously).

We’ll use the Greenland 5 km DEM, Ice Thickness, and Bedrock Elevation Grids (J. Bamber 2001) from J. L. Bamber, Layberry, and Gogineni (2001) and Layberry and Bamber (2001). Download here (after registration).

< section id="data" class="level1">

Data

The data needs some wrangling as the format is not straightforward: it’s a wrapped fixed width ASCII file (check the user guide). We need to make one row out of every 31 lines of the file, reverse the order of the lines and give the correct projection and extent.

library(terra)
library(readr)
library(dplyr)
library(tidyr)

thick <- read_fwf("thick_5km_corrected",
                  guess_max = 1e4) |> 
  mutate(row = ceiling(row_number() / 31)) |> 
  group_by(row) |> 
  group_modify(~ as_tibble(as.vector(t(as.matrix(.x))))) |> 
  ungroup() |> 
  mutate(name = rep(paste0("x", 1:310), 561)) |> 
  drop_na(value) |> 
  pivot_wider(values_from = value,
              names_from = name) |>
  arrange(desc(row)) |> 
  select(-row) |> 
  as.matrix() |> 
  rast(crs = "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs") 

ext(thick) = c(-800000, 700000, -3400000, -600000)
< section id="map" class="level1">

Map

Here is the raw map in a polar stereographic projection:

thick |> 
  plot(main = "Greenland Ice thickness",
       col = map.pal("magma"))
Figure 1: Greenland Ice thickness in a polar stereographic projection. data: Bamber J., 2021. NASA National Snow and Ice Data Center

And on an interactive map after reprojection:

library(leaflet)

# native resolution is 5 km, so at 66° N it's about 0.1 degree
# -75 - -15 = 65 ; 65 / 0.1 = 650 pixels wide
thick_wgs84 <- thick |>
  project(rast(nrows = 250, ncols = 650, 
               xmin = -75, xmax = -10, 
               ymin = 60, ymax = 85,
               crs = "EPSG:4326")) |> 
  subst(x = _, 0, NA)

range_m <- c(0, max(values(thick_wgs84), na.rm = TRUE))
pal <- colorNumeric("magma", domain = range_m, na.color = "#ffffff00")
# fix the legend order:
pal_rev <- colorNumeric("magma", domain = range_m, na.color = "#ffffff00", reverse = TRUE)

leaflet() |> 
  addTiles() |> 
  addRasterImage(thick_wgs84, 
                 colors = pal,
                 opacity = 0.5,
                 attribution = "Bamber, 2021 (NASA National Snow and Ice Data Center)") |> 
  addLegend(pal = pal_rev,
            title = "Greenland<br />Ice thickness (m)",
            values = range_m,
            labFormat = labelFormat(transform = function(x) sort(x + 500, decreasing = TRUE)))
To leave a comment for the author, please follow the link and comment on their blog: r.iresmi.net.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Exit mobile version