Add ice edge to maps

[nmbool]

Hi Mike,

It would be great if you could cover how to add the mean ice extent in the Antarctic for September (a line) to my maps. cheers,

Nat

[mdsumner]

This is a good example, there's a lot to talk about. Can you answer

• mean of many Septembers, or just one specific one?
• mean latitude extent of all days in September, or the contour of the mean monthly conc ?...
• what are your maps like? do they have a custom map projection?

Sea ice extent line on a map

(This is raadtools specific, so not available to everyone without setting that up. )

Plot in longlat, note that we use concentration 15% to contour at, and this is of mean monthly conc.

library(raadtools)

## a particular extent in longlat
ex <- extent(100, 160, -70, -50)

## a specific September, we project it to longlat so the line "wraps" properly
ice <- projectRaster(readice("2016-09-15", time.resolution = "monthly"), crs = "+proj=longlat +ellps=WGS84")
line <- rasterToContour(ice, level = 15)

## plot your map (or whatever map)
plot(ex, asp = 1.2)
plot(line, add = TRUE)

We might:

• show this in leaflet
• show different ways of extracting the contours, problems with the dateline etc.

[mdsumner]

Here is the actual study map we are using. I'm working with generic pieces rather than our model output, and this one is pretty involved.

## a custom projection
prj <- "+proj=omerc +lonc=165 +lat_0=-22 +alpha=23 +k=0.99984 +x_0=0 +y_0=0 +no_uoff +gamma=23 +ellps=WGS84 +units=m +no_defs"
## the longlat basis
llprj <-  "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
library(raster)
## a background grid for the region
master <- raster(extent(-16800000, 12500000, -20300000, 19200000), res = 150000, crs = prj)

library(raadtools)
## get north and south hemi ice 
sice <- readice("2016-09-15", time.resolution = "monthly")
nice <- readice("2016-09-15", time.resolution = "monthly", hemisphere = "north")
sline <- spTransform(rasterToContour(sice, level = 15), prj)
nline <- spTransform(rasterToContour(nice, level = 15), prj)

## clean up the lines
sline <- disaggregate(sline); sline <- sline[which.max(rgeos::gLength(sline, byid = TRUE)), ]
nline <- disaggregate(nline); nline <- nline[which.max(rgeos::gLength(nline, byid = TRUE)), ]

library(graticule)
g1 <- graticule(seq(85,250 , by = 15), seq(-85, 0, by = 15), ylim = c(-85, 0), proj = prj)
g2 <- graticule(seq(100, 250, by = 15), seq(85, 25, by = -15), ylim = c(85, 0), proj = prj)

library(maptools)
data(wrld_simpl)
coords <- coordinates(spTransform(as(as(wrld_simpl, "SpatialLinesDataFrame"), "SpatialPointsDataFrame"), prj))
op <- par(mar = rep(0, 4))
plot(NA, xlim = c(xmin(master), xmax(master)), ylim = c(ymin(master), ymax(master)), axes = FALSE, xlab = "", ylab = "", asp = 1)

points(coords, pch = ".", col = "black")
plot(g1, add = TRUE, lty = 2, col = "darkgrey"); plot(g2, add = TRUE, lty = 2, col = "darkgrey")

plot(sline, add = TRUE, col = "dodgerblue")
plot(nline, add = TRUE, col = "dodgerblue")