Here we calculate sea ice season for the southern hemisphere, after the method of Massom et al. (2013).
The northern hemisphere requires a different approach, so here we set up the files for testing but leave the application for later.
We use the R package raadtools which provides access to a collection of all NSIDC 25 km sea ice concentration for both southern and northern hemispheres. Reading the ice data is provided by the readice function which accepts a vector of dates, other options and returns a Raster object.
First we pre-cache a single file of all daily data for each year after doing some sanity checks. This allows ease of code development and debugging.
library(raadtools) library(dplyr) sfiles <- icefiles(hemisphere = "south") sfiles %>% summarize(count = n(), datemin = min(date), datemax = max(date))
Data are missing from 1987-12-04 to early Jan 1988, Massom et al. use "the climatology" but it seems pertinent to just assume we cannot calculate day of retreat for that season since the mean value can be wildly different from a given year.
Data are two-daily from 1978 to mid 80s, we just interpolate to fill the days.
dp <- "/mnt/acebulk/seaiceseason" ## function to ensure we have clean boundaries (untidy for now) snaptodoy <- function(date, mon = 2, mday = 15, start = TRUE) { ldate <- as.POSIXlt(date) doy <- as.integer(format(ISOdatetime(ldate$year + 1900, mon, mday, 0, 0, 0, tz = "UTC"), "%j")) test <- ldate$yday > doy if (start) { if (test) {ldate$mday <- mday; ldate$mon <- mon - 1; ldate$year <- ldate$year + 1} if (!test) {ldate$mday <- mday; ldate$mon <- mon - 1} } else { if (test) {ldate$mday <- mday; ldate$mon <- mon - 1} if (!test) {ldate$mday <- mday; ldate$mon <- mon + 1; ldate$year <- ldate$year - 1} } as.POSIXct(ldate) } ## function to read a series of time slices, and interpolate a full sequence optionally readIceBrick <- function(dates, hemisphere = "south", interpolate = TRUE) { x <- suppressWarnings(readice(dates, hemisphere = hemisphere, setNA = FALSE)) ## we need to interpolate (or insert monthly clim) if (!interpolate) return(x) if (nlayers(x) == length(dates)) return(x) xmat <- values(x) xmat[xmat > 100] <- 0 dimnames(xmat) <- list(NULL, NULL) ## rebuild the entire data set with linear interpolation xmat2 <- matrix(0, nrow(xmat), length(dates)) t0 <- timedateFrom(getZ(x)) for (i in seq_len(nrow(xmat2))) { xmat2[i,] <- approxfun(t0, xmat[i,], rule = 2)(dates) } x <- setValues(x, xmat2) ##names(x) <- sprintf("layer%0.3i", seq(nlayers(x))) names(x) <- format(dates, "%Y_%m_%d") setZ(x, dates) } calc_ice_season <- function(yfile, threshval = 15) { hemi <- substr(basename(yfile[1]), 1, 5) threshold.value <- threshval ndays <- 5 len <- 15 ## north_ice_1979_02_15_1980_02_14.grd" ice <- brick(yfile) year_n <- nlayers(ice) template <- ice[[1]] * 0 icemat <- values(ice) icemat[icemat > 100] <- 0 ## here we need to get the next day for the interpolation . . . ##icemat[is.na(icemat)] <- 0 adv <- numeric(nrow(icemat)) ret <- numeric(nrow(icemat)) threshold <- icemat >= threshold.value rsum <- .rowSums(threshold, nrow(threshold), ncol(threshold)) ## all values less alllt <- rsum == 0 ## all values greater than threshold allgt <- rsum == ncol(threshold) ## values missing ##miss <- icemat[,1] > 100 ## all the rest visit <- which(!alllt & !allgt) for (ii in seq_along(visit)) { rl <- rle(threshold[visit[ii], ]) ## annual day of advance is the time when the ice ## concentration in a given pixel first exceeds 15% (taken to ## approximate the ice edge) for at least 5 days for (ir in seq_along(rl$lengths)) { if (rl$values[ir] & rl$lengths[ir] >= ndays) { adv[visit[ii]] <- if (ir == 1) 1L else sum(head(rl$lengths, ir - 1)) break; } } if (adv[visit[ii]] == 0) adv[visit[ii]] <- NA ##while day of retreat ## ## is the time when concentration remains below 15% until the end ## ## of the given sea ice year revlengths <- rev(rl$lengths) revvals <- rev(rl$values) for (ri in seq_along(revlengths)) { if (revvals[ri]) { ret[visit[ii]] <- if (ri == 1) length(year_n) else sum(revlengths[ri:length(revlengths)]) break; } } if (ret[visit[ii]] == 0) ret[visit[ii]] <- NA } adv[alllt] <- NA ## adv[miss] <- NA adv[allgt] <- 1 ret[alllt] <- NA ret[allgt] <- length(year_n) list(adv = setValues(template, adv), ret = setValues(template, ret)) }
Collate a single file for each year, in a local folder. Years start at February 15 and end on February 14 in the southern hemisphere, and start August 15 and end August 14 in the northern hemisphere.
hmon <- 2 hemi <- "south" hmday <- 15 allfiles <- sfiles startdate <- snaptodoy(min(allfiles$date), mon = hmon, mday = hmday, start = TRUE) endate <- snaptodoy(max(allfiles$date), mon = hmon, mday = hmday -1, start = FALSE) ## start dates, including the end of the last year bdates <- seq(startdate, endate, by = "1 year") for (iyear in seq_along(head(bdates, - 1))) { #yearn <- (ydates %>% filter(yearN == levels(yearN)[iyear]))$date yearn <- seq(bdates[iyear], bdates[iyear + 1], by = "1 day") tfilename <- file.path(dp, "yearfiles", sprintf("%s_ice_%s_%s.grd", hemi, format(min(yearn), "%Y_%m_%d"), format(max(yearn), "%Y_%m_%d"))) print(iyear) if (!file.exists(tfilename)) { ice <- readIceBrick(yearn, hemisphere = hemi, interpolate = TRUE) writeRaster(ice, tfilename) } }
Ice season is defined by a calculated day of advance and retreat for a pixel corresponding to the development of sea ice cover and its subsequent melting.
The annual day of advance is the time when the ice concentration in a given pixel first exceeds 15% (taken to < approximate the ice edge) for at least 5 days while day of retreat is the time when concentration remains below 15% until the end of the given sea ice year. Massom et al (2013)
To calculate this we arrange an entire year of data in a single 3-dimensional array of X * Y * Time pixels. Where there is a break in the Time dimension each pixel's sequence in time is linearly interpolated to each calendar day. This occurs in the first 12???? years of ice data where there is only one data set every two days, and once in ????1998 where there is 1.5 month gap. (Note that our method diverges from that of Massom et al. who use the monthly climatology in this case).
The method proceeds to threshold the ice values at 15% or greater concentration, then processing the sequence base on the rules in Massom et al. (2103). In short:
library(raadtools) files <- list.files(file.path(dp, "yearfiles"), full.names = TRUE, pattern = "^south.*grd$") advance <- vector("list", length(files)) retreat <- vector("list", length(files)) for (i in seq_along(files)) { ## obj contains 2 rasters "adv" and "ret" obj <- calc_ice_season(files[i]) advance[[i]] <- obj$adv retreat[[i]] <- obj$ret } advance <- brick(stack(advance), filename = file.path(dp, "south_advance.grd")) retreat <- brick(stack(retreat), filename = file.path(dp, "south_retreat.grd"))
Load the advance and retreat layers and calculate overall mean and variability.
advance <- brick(file.path(dp, "south_advance.grd")) retreat <- brick(file.path(dp, "south_retreat.grd")) isd <- calc(retreat - advance, sd, na.rm = TRUE) msd <- calc(retreat - advance, median, na.rm = TRUE) plot(msd, col = viridis::viridis(100), zlim = c(0, 365)) jet.colors_alt <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#FF7F00", "red", "#7F0000")) plot(isd, col = jet.colors_alt(12), breaks = seq(0, 50, by = 5), zlim = c(0, 50))
Massom R, Reid P, Stammerjohn S, Raymond B, Fraser A, et al. (2013) Change and Variability in East Antarctic Sea Ice Seasonality, 1979/80–2009/10. PLoS ONE 8(5): e64756. doi:10.1371/journal.pone.0064756
sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.