In AustralianAntarcticDivision/aceecostats: Ecosystem status and trends

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.

Tools

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.

Preliminaries

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

ISSUES

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.

Utilities.

dp <- "/rdsi/PRIVATE/raad/data_local/aad.gov.au/iceseason"

## 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 <- raster::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)
  }
}

Calculating ice season

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:

advance is the first day at or above the threshold after which the subsequent four days are also at or above the threshold
retreat is the day after which the threshold is not reached again.

library(raadtools)
files <- list.files(file.path("yearfiles"), full.names = TRUE, pattern = "^south.*grd$")

advance <- vector("list", length(files))
retreat <- vector("list", length(files))
daycount <- 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
  daycount[[i]] <- calc(readAll(brick(files[i])) > 15, mean, na.rm = TRUE) * 365

}


#advance <- brick(stack(advance), filename = file.path(dp, "south_advance.grd"))
#retreat <- brick(stack(retreat), filename = file.path(dp, "south_retreat.grd"))
#daycount <- brick(stack(daycount), filename = file.path(dp, "south_daycount.grd"))

advance <- readAll(brick(stack(advance)))
retreat <- readAll(brick(stack(retreat)))
daycount <- readAll(brick(stack(daycount)))

saveRDS(advance, file.path(dp, "ice_advance.rds"))
saveRDS(retreat, file.path(dp, "ice_retreat.rds"))
saveRDS(daycount, file.path(dp, "ice_daycount.rds"))

at 70%

library(raadtools)
files <- list.files(file.path("yearfiles"), full.names = TRUE, pattern = "^south.*grd$")

advance <- vector("list", length(files))
retreat <- vector("list", length(files))
daycount <- vector("list", length(files))

for (i in seq_along(files)) {
  ## obj contains 2 rasters "adv" and "ret"
  obj <- calc_ice_season(files[i], threshval = 70)
  advance[[i]] <- obj$adv
  retreat[[i]] <- obj$ret
  daycount[[i]] <- calc(readAll(brick(files[i])) > 70, mean, na.rm = TRUE) * 365

}


#advance <- brick(stack(advance), filename = file.path(dp, "south_advance.grd"))
#retreat <- brick(stack(retreat), filename = file.path(dp, "south_retreat.grd"))
#daycount <- brick(stack(daycount), filename = file.path(dp, "south_daycount.grd"))

advance <- readAll(brick(stack(advance)))
retreat <- readAll(brick(stack(retreat)))
daycount <- readAll(brick(stack(daycount)))

saveRDS(advance, file.path(dp, "ice_70_advance.rds"))
saveRDS(retreat, file.path(dp, "ice_70_retreat.rds"))
saveRDS(daycount, file.path(dp, "ice_70_daycount.rds"))

at 95%

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))
daycount <- vector("list", length(files))

library(future.apply)
plan(multiprocess)
dofun <- function(file) {
  obj <- calc_ice_season(file, threshval = 95)
  list(adv = obj$adv, ret = obj$ret, count = calc(readAll(brick(file)) > 95, mean, na.rm = TRUE) * 365)
}
out <- future_lapply(files, dofun)

# for (i in seq_along(files)) {
#   ## obj contains 2 rasters "adv" and "ret"
#   obj <- calc_ice_season(files[i], threshval = 95)
#   advance[[i]] <- obj$adv
#   retreat[[i]] <- obj$ret
#   daycount[[i]] <- calc(readAll(brick(files[i])) > 95, mean, na.rm = TRUE) * 365
#   print(i)
# }

advance <- readAll(brick(stack(lapply(out, function(aa) aa$adv))))
retreat <- readAll(brick(stack(lapply(out, function(aa) aa$ret))))
daycount <- readAll(brick(stack(lapply(out, function(aa) aa$count))))

saveRDS(advance, file.path(dp, "ice_95_advance.rds"))
saveRDS(retreat, file.path(dp, "ice_95_retreat.rds"))
saveRDS(daycount, file.path(dp, "ice_95_daycount.rds"))

Exploring ice season

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

References

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