inst/doc/gridded-event-detection.R

In RmarineHeatWaves: Detect Marine Heat Waves and Marine Cold Spells

## ----global_options, include = FALSE-------------------------------------
knitr::opts_chunk$set(fig.width = 4, fig.align = 'center',
                      echo = FALSE, warning = FALSE, 
                      message = FALSE, tidy = FALSE)

## ----load-pkg, echo = FALSE, eval = FALSE--------------------------------
#  library(RmarineHeatWaves)
#  library(ncdf4)
#  library(plyr)
#  library(dplyr)
#  library(tibble)
#  library(reshape2)
#  library(lubridate)
#  library(ggplot2)
#  library(doMC); doMC::registerDoMC(cores = 4)

## ----read-netcdf, echo = TRUE, eval = FALSE------------------------------
#  nc <- nc_open("/Users/ajsmit/spatial/OISSTv2/daily/OISST-V2-AVHRR_agg.nc")
#  # print(nc)
#  sst <- ncvar_get(nc, varid = "sst") # reads in the only var in netCDF file
#  dimnames(sst) <- list(lon = nc$dim$lon$vals,
#                        lat = nc$dim$lat$vals,
#                        t = nc$dim$time$vals)

## ----sst-df, echo = TRUE, eval = FALSE-----------------------------------
#  sst <- as_tibble(melt(sst, value.name = "temp"))
#  sst$t <- as.Date(sst$t, origin = '1978-01-01')

## ----remove-NA, echo = TRUE, eval = FALSE--------------------------------
#  sst1 <- sst %>%
#    na.omit()
#  rm(sst); rm(nc)

## ----detect-fun, echo = TRUE, eval = FALSE-------------------------------
#  OISST_detect <- function(dat) {
#    dat <- dat[,3:4]
#    start = "1983-01-01"
#    end = "2012-12-31"
#    whole <- make_whole(dat)
#    mhw <- detect(whole, climatology_start = start, climatology_end = end,
#                  min_duration = 5, max_gap = 2, cold_spells = FALSE)
#    events <- mhw$event
#    return(events)
#    }

## ----apply-detect-fun, echo = TRUE, eval = FALSE-------------------------
#  # it takes a long time...
#  system.time(OISST_events <- ddply(sst1, .(lon, lat), OISST_detect, .parallel = TRUE))
#   #   user  system elapsed
#   # 1862.3   205.9   704.8
#  save(OISST_events, file = "/Users/ajsmit/spatial/OISSTv2/MHWs/OISST_events.Rdata")

## ----load-SST, echo = TRUE, eval = FALSE---------------------------------
#  # in case it was calculated and saved earlier, we can load it here now:
#  load("/Users/ajsmit/spatial/OISSTv2/MHWs/OISST_events.Rdata")

## ----event-tally, echo = TRUE, eval = FALSE------------------------------
#  # summarise the number of unique longitude, latitude and year combination:
#  event_freq <- OISST_events %>%
#    mutate(year = year(date_start)) %>%
#    group_by(lon, lat, year) %>%
#    summarise(n = n())
#  head(event_freq)
#  
#  # create complete grid for merging with:
#  sst.grid <- sst1 %>%
#    select(lon, lat, t) %>%
#    mutate(t = lubridate::year(t)) %>%
#    unique() # slow...
#  colnames(sst.grid)[colnames(sst.grid) == 't'] <- 'year'
#  
#  # and merge:
#  OISST_n <- left_join(sst.grid, event_freq)
#  OISST_n[is.na(OISST_n)] <- 0

## ----trend-fun, echo = TRUE, eval = FALSE--------------------------------
#  lin_fun <- function(ev) {
#    mod1 <- glm(n ~ year, family = poisson(link = "log"), data = ev)
#    # mod1 <- lm(n ~ year, data = ev)
#    tr <- summary(mod1)$coefficients[2,c(1,4)] # extract slope coefficient and its p-value
#    return(tr)
#  }

## ----apply-trend-fun, echo = TRUE, eval = FALSE--------------------------
#  OISST_nTrend <- ddply(OISST_n, .(lon, lat), lin_fun)
#  OISST_nTrend$pval <- cut(OISST_nTrend[,4], breaks = c(0, 0.001, 0.01, 0.05, 1))
#  head(OISST_nTrend)

## ----load-geography, echo = TRUE, eval = FALSE---------------------------
#  ## Coastline of African continent and some borders:
#  load("/Users/ajsmit/Dropbox/repos/Trend_Analysis/graph/africa_coast.RData")
#  load("/Users/ajsmit/Dropbox/repos/Trend_Analysis/graph/africa_borders.Rdata")
#  load("/Users/ajsmit/Dropbox/repos/Trend_Analysis/graph/south_africa_coast.RData")
#  load("/Users/ajsmit/Dropbox/repos/Trend_Analysis/graph/sa_provinces_new.RData")

## ----the-figure, echo = TRUE, eval = FALSE-------------------------------
#  ggplot(OISST_nTrend, aes(x = lon, y = lat)) +
#    geom_rect(size = 0.2, fill = NA,
#         aes(xmin = lon - 0.1, xmax = lon + 0.1, ymin = lat - 0.1, ymax = lat + 0.1,
#             colour = OISST_nTrend[,5])) +
#    geom_raster(aes(fill = Estimate), interpolate = FALSE, alpha = 0.9) +
#    scale_fill_gradient2(name = "Count\n per\n year", high = "red", mid = "white",
#                         low = "darkblue", midpoint = 0) +
#    scale_colour_manual(breaks = c("(0,0.001]", "(0.001,0.01]", "(0.01,0.05]", "(0.05,1]"),
#                        values = c("firebrick1", "firebrick2", "firebrick3", "white"),
#                        name = "p-value", guide = FALSE) +
#    geom_polygon(data = south_africa_coast, aes(x = lon, y = lat, group = group),
#                 colour = NA, fill = "grey80") +
#    geom_path(data = africa_borders, aes(x = lon, y = lat, group = group),
#              size = 0.3, colour = "black") +
#    geom_polygon(data = south_africa_coast, aes(x = lon, y = lat, group = group),
#                 size = 0.5, colour = "black", fill = NA) +
#    coord_fixed(ratio = 1, xlim = c(12.12, 35.12), ylim = c(-26.88, -37.38), expand = TRUE) +
#    theme_bw()
#  ggsave(file = "README-grid-example1.png", width = 6, height = 3.2)
#  
#  ggplot(OISST_nTrend, aes(lon, lat)) +
#    geom_raster(aes(fill = pval), interpolate = FALSE) +
#    scale_fill_manual(breaks = c("(0,0.001]", "(0.001,0.01]", "(0.01,0.05]",
#                                 "(0.05,0.1]", "(0.1,0.5]", "(0.5,1]"),
#                      values = c("black", "grey20", "grey40",
#                                 "grey80", "grey90", "white"),
#                      name = "p-value") +
#    geom_polygon(data = south_africa_coast, aes(x = lon, y = lat, group = group),
#                 colour = NA, fill = "grey80") +
#    geom_path(data = africa_borders, aes(x = lon, y = lat, group = group),
#              size = 0.3, colour = "black") +
#    geom_polygon(data = south_africa_coast, aes(x = lon, y = lat, group = group),
#                 size = 0.5, colour = "black", fill = NA) +
#    coord_fixed(ratio = 1, xlim = c(12.12, 35.12), ylim = c(-26.88, -37.38), expand = TRUE) +
#    theme_bw()
#  ggsave(file = "README-grid-example2.png", width = 6, height = 3.2)