How to create contour with wind animation using gganimate?

Question

I have created a wind animation for each hour of a specific day (Click to play animation video). Instead of showing 19 points, I would like to create a contour plot interpolating/extrapolating using those 19 points at each hour over the entire area - just like this contour map produced using ArcGIS and spline interpolation tool.

The following code shows ggplot and gganimate I used to create the hourly wind animation. I have only created a small data frame as a subsample of the complete 24-hr data as I'm not familiar with attaching a csv into this forum. Is there a way to include a contour overlaying the area instead of geom_point?

library(ggplot2)
library(ggmap)
library(gganimate)

site <- c(1:18, 1:18)    
date <- data.frame(date=c(rep(as.POSIXct("2018-06-07 00:00:00"),18),rep(as.POSIXct("2018-06-07 01:00:00"),18)))    
long <- c(171.2496,171.1985, 171.2076, 171.2236,171.2165,171.2473,171.2448,171.2416,171.2243,171.2282,171.2344,171.2153,171.2532,171.2444,171.2443,171.2330,171.2356,171.2243)   
lati <- c(-44.40450,-44.38520,-44.38530,-44.38750,-44.39195,-44.41436,-44.38798,-44.38934,-44.37958,-44.37836,-44.37336,-44.37909,-44.40801, -44.40472,-44.39558,-44.40971,-44.39577,-44.39780)    
PM <- c(57,33,25,48,34,31,52,48,31,51,44,21,61,53,49,34,60,18,41,26,28,26,26,18,32,28,27,29,22,16,34,42,37,28,33,9)    
ws <- c(0.8, 0.1, 0.4, 0.4, 0.2, 0.1, 0.4, 0.2, 0.3, 0.3, 0.2, 0.7, NaN, 0.4, 0.3, 0.4, 0.3, 0.3, 0.8, 0.2, 0.4, 0.4, 0.1, 0.5, 0.5, 0.2, 0.3, 0.3, 0.3, 0.4, NaN, 0.5, 0.5, 0.4, 0.3, 0.2)    
wd <- c(243, 274, 227, 253, 199, 327, 257, 270, 209, 225, 230, 329, NaN, 219, 189, 272, 239, 237, 237, 273, 249, 261, 233, 306, 259, 273, 218, 242, 237, 348, NaN, 221, 198, 249, 236,252  )    
PMwind <- cbind(site,date,long,lati,PM, ws, wd)

tmlat <- c(-44.425, -44.365)                
tmlon <- c(171.175, 171.285)  

tim <- get_map(location = c(lon = mean(tmlon), lat = mean(tmlat)),
               zoom = 14,
               maptype = "terrain")

ggmap(tim) + 
    geom_point(aes(x=long, y = lati, colour=PM), data=PMwind,
               size=3,alpha = .8, position="dodge", na.rm = TRUE) +     
    geom_spoke(aes(x=long, y = lati, angle = ((270 -  wd) %% 360) * pi / 180), data=PMwind, 
               radius = -PMwind$ws * .01, colour="yellow", 
               arrow = arrow(ends = "first", length = unit(0.2, "cm"))) +
    transition_states(date, transition_length = 20, state_length = 60) +
    labs(title = "{closest_state}") +
    ease_aes('linear', interval = 0.1) +
    scale_color_gradient(low = "green", high = "red")+
    theme_minimal()+
    theme(axis.text.x=element_blank(), axis.title.x=element_blank()) +
    theme(axis.text.y=element_blank(), axis.title.y=element_blank()) +
    shadow_wake(wake_length = 0.01)

Answer 1

This can be done, but I'd say it's far from straightforward with current tools. To go from the dataset in the question to animated contours, we'll need to address the following obstacles:

1. We only have a couple of data points, spread out irregularly over the given area. Contour generation usually expects a regular grid of points.

2. geom_contour/stat_contour in ggplot2 don't deal well with open contours at the edges. See here for a GH discussion of what happens when we try to use the contour lines for filled polygons.

3. The polygons associated with contours don't necessarily persist over time: they appear, disappear, split into multiple smaller polygons, merge into bigger polygons, etc. This makes things difficult in gganimate, which needs to know which elements in frame T corresponds to which elements in frame T+1, in order to interpolate them properly.

The first two obstacles can be addressed with existing workarounds. The third requires some unorthodox hacking.

Part 1: Interpolate irregular points

Take the longitude / latitude / PM values of PMwind for each date value & use interp from the akima package to interpolate them into a regular grid. Bicubic spline interpolation with extrapolation set to TRUE will give you a regular grid of 40 x 40 (by default, change nx / ny parameter values if you prefer the grid to be coarser / finer) points with interpolated PM values.

library(dplyr)

PMwind2 <- PMwind %>%
  select(date, long, lati, PM) %>%
  tidyr::nest(-date) %>%
  mutate(data = purrr::map(data,
                           ~ akima::interp(x = .$long, y = .$lati, z = .$PM,
                                           linear = FALSE, # use spline interpolation
                                           extrap = TRUE) %>%
                             akima::interp2xyz(data.frame = TRUE))) %>%
  tidyr::unnest()

> str(PMwind2) # there are 2 x 40 x 40 observations, corresponding to 2 dates
'data.frame':   3200 obs. of  4 variables:
 $ date: POSIXct, format: "2018-06-07" "2018-06-07" "2018-06-07" ...
 $ x   : num  171 171 171 171 171 ...
 $ y   : num  -44.4 -44.4 -44.4 -44.4 -44.4 ...
 $ z   : num  31.8 31.4 31 30.6 30.3 ...

Part 2: Use an alternative package for generating contours with closed polygons at the edges.

Here, I used geom_contour_fill from the metR package, which is one of the fixes discussed in the GH thread. (The isoband package approach looks interesting too, but it's newer, & I haven't tested it out yet.)

library(ggplot2)
library(metR)

# define scale breaks to make sure the scale would be consistent across animated frames
scale.breaks = scales::fullseq(range(PMwind2$z), size = 10)

# define annotation layer & appropriate coord limits for map (metR's contour polygons
# don't go nicely with alpha < 1 in animation, as the order of layers could change, 
# but we can overlay the map as a semi-transparent annotation layer over the contour
# polygons, instead of having ggmap layer beneath semi-transparent contour polygons.)
map.annotation <- list(
  annotation_raster(tim %>% unlist() %>%
                      alpha(0.4) %>% # change alpha setting for map here
                      matrix(nrow = dim(tim)[1], 
                             byrow = TRUE),
                    xmin = attr(tim, "bb")$ll.lon,
                    xmax = attr(tim, "bb")$ur.lon,
                    ymin = attr(tim, "bb")$ll.lat,
                    ymax = attr(tim, "bb")$ur.lat),
  coord_quickmap(xlim = c(attr(tim, "bb")$ll.lon, attr(tim, "bb")$ur.lon),
                 ylim = c(attr(tim, "bb")$ll.lat, attr(tim, "bb")$ur.lat),
                 expand = FALSE))

p.base <- ggplot(PMwind2, aes(x = x, y = y, z = z))

# check static version of plot to verify that the geom layer works as expected
p.base + 
  geom_contour_fill(breaks = scale.breaks) +
  facet_wrap(~date) +
  map.annotation +
  scale_fill_gradient(low = "green", high = "red",
                      aesthetics = c("colour", "fill"),
                      limits = range(scale.breaks)) +
  theme_minimal()

Part 3: Instead of animating contour lines / polygons, animate the point values

After each frame of the animated plot has been generated (but before it has been printed / drawn to a graphics device), take its data, create a new plot (the plot we actually want), & send that to the graphics device instead. We can do so by inserting some code into plot_frame, a function in the ggproto object gganimate:::Scene, where the plotting happens.

Scene2 <- ggproto(
  "Scene2", gganimate:::Scene,
  plot_frame = function(self, plot, i, newpage = is.null(vp), vp = NULL, 
                        widths = NULL, heights = NULL, ...) {    
    plot <- self$get_frame(plot, i)

    # for each frame, use the plot data interpolated by gganimate to create a new plot
    new.plot <- ggplot(data = plot[["data"]][[1]],
                       aes(x = x, y = y, z = z)) + 
      geom_contour_fill(breaks = scale.breaks) +
      ggtitle(plot[["plot"]][["labels"]][["title"]]) +
      map.annotation +
      scale_fill_gradient(low = "green", high = "red",
                          limits = range(scale.breaks)) +
      theme_minimal()
    plot <- ggplotGrob(new.plot)

    # no change below
    if (!is.null(widths)) plot$widths <- widths
    if (!is.null(heights)) plot$heights <- heights
    if (newpage) grid::grid.newpage()
    grDevices::recordGraphics(
      requireNamespace("gganimate", quietly = TRUE),
      list(),
      getNamespace("gganimate")
    )
    if (is.null(vp)) {
      grid::grid.draw(plot)
    } else {
      if (is.character(vp)) seekViewport(vp)
      else pushViewport(vp)
      grid::grid.draw(plot)
      upViewport()
    }
    invisible(NULL)
  })

We also need to define a series of intermediary functions in order for the animation to use this Scene2 instead of the original gganimate:::Scene. I've used the same approach to answer another question before here, with some discussion on the pros & cons of doing so.

library(magrittr)

create_scene2 <- function(transition, view, shadow, ease, transmuters, nframes) {
  if (is.null(nframes)) nframes <- 100
  ggproto(NULL, Scene2, transition = transition, 
          view = view, shadow = shadow, ease = ease, 
          transmuters = transmuters, nframes = nframes)
}

ggplot_build2 <- gganimate:::ggplot_build.gganim
body(ggplot_build2) <- body(ggplot_build2) %>%
  as.list() %>%
  inset2(4,
         quote(scene <- create_scene2(plot$transition, plot$view, plot$shadow, 
                                      plot$ease, plot$transmuters, plot$nframes))) %>%
  as.call()

prerender2 <- gganimate:::prerender
body(prerender2) <- body(prerender2) %>%
  as.list() %>%
  inset2(3,
         quote(ggplot_build2(plot))) %>%
  as.call()

animate2 <- gganimate:::animate.gganim
body(animate2) <- body(animate2) %>%
  as.list() %>%
  inset2(7,
         quote(plot <- prerender2(plot, nframes_total))) %>%
  as.call()

Finally, here's the result:

library(gganimate)

animate2(p.base + 
           geom_point(aes(color = z)) + # this layer will be replaced by geom_contour_fill in 
                                        # the final plot; it's here as the placeholder in 
                                        # order for gganimate to interpolate the relevant data
           transition_time(date) +
           ggtitle("{frame_time}"),
         nframes = 30, fps = 10)        # you can increase nframes for smoother transition
                                        # (which would also be much bigger in file size)