Site icon R-bloggers

Global sea ice time series visualization sandbox with R, ggplot, and plotly

[This article was first published on Maxwell B. Joseph, 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.

In the past couple of weeks I’ve noticed a flurry of visualizations of global sea ice extent on social media. If you’re like me, and curious to see what the data look like yourself, here’s a bit of R code to fetch and visualize the most recent data.

library(dplyr)
library(ggplot2)
library(viridis)
library(plotly)


# Define helper functions to get and parse sea ice data from NSIDC --------
read_nsidc_csv <- function(file) {
  # reads and loads a csv file of sea ice extent from NSIDC's ftp server
  hem_code <- substr(file, 1, 2)
  if (hem_code == "NH") {
    hemisphere <- "north"
  } else if (hem_code == "SH") {
    hemisphere <- "south"
  } else {
    stop("Filename is not correctly formatted")
  }
  prefix <- paste0("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/",
                    hemisphere,
                   "/daily/data/")
  file_loc <- paste0(prefix, file)
  d <- read.csv(file_loc, stringsAsFactors = FALSE)
  d <- d[-1, ] # this row contains the units

  maybe_as_numeric <- function(x) {
    numeric_x <- suppressWarnings(as.numeric(x))
    if (!all(is.na(numeric_x))) x <- numeric_x
    x
  }

  # convert numeric columns to numbers and return data.frame
  lapply(d, maybe_as_numeric) %>%
    as.data.frame(stringsAsFactors = FALSE)
}


get_seaice_data <- function() {
  # gets seaice data from northern and southern hemispheres and joins to one df
  northern_hem <- full_join(read_nsidc_csv("NH_seaice_extent_final_v2.csv"),
                              read_nsidc_csv("NH_seaice_extent_nrt_v2.csv")) %>%
    mutate(hem = "Northern hemisphere")
  southern_hem <- full_join(read_nsidc_csv("SH_seaice_extent_final_v2.csv"),
                            read_nsidc_csv("SH_seaice_extent_nrt_v2.csv")) %>%
    mutate(hem = "Southern hemisphere")
  all_data <- full_join(northern_hem, southern_hem)
  all_data %>%
    mutate(date = ISOdate(Year, Month, Day),
           day_of_year = as.numeric(format(date, "%j")))
}


# load data from NSIDC's ftp server
seaice <- get_seaice_data()



# Visualize ---------------------------------------------------------------
seaice %>%
  ggplot(aes(x = day_of_year, y = Extent, color = Year, group = Year)) +
  geom_line() +
  xlab("Day of year") +
  ylab("Sea ice extent (million sq km)") +
  scale_color_viridis(direction = -1) +
  geom_text(aes(x = day_of_year, y = Extent, label = Year),
            data = subset(seaice, date == max(date)),
            nudge_x = c(15, -17)) +
  ggtitle(paste("Global sea ice extent:",
                 min(format(seaice$date, "%Y-%m-%d")),
                 "to",
                 max(format(seaice$date, "%Y-%m-%d")))) +
  facet_wrap(~ hem)
  ggplotly()

Click here to view a larger interactive plot on RPubs

Disclaimer: no part of this post is sponsored by the National Snow and Ice Data Center, which serves the data. This post just shows how to fetch and plot the data.

To leave a comment for the author, please follow the link and comment on their blog: Maxwell B. Joseph.

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.