Skip to content

[Code] Animating a spatio-temporal dataset over a geographic map in R

December 7, 2012

[This is a code post. I wrote this code, and it did something useful for me. Feel free to use the code for your own purposes. If it’s useful, don’t hesitate to tell me you liked it!]

I wrote this code to generate an animation of seabird recovery locations after the 2011 Rena oil spill. It generates a .gif animation where each frame is a map of seabird recovery locations on that day of the response effort, with datapoints colour-coded according to how heavily oiled the bird was. Colour-coding is based on a ramp where zero scores are grey and all positive-number scores follow a colour-ramp to easily discriminate oiled from non-oiled birds. The code should generalise to any situation where you need an animated map showing the locations of datapoints through time, where the colour of the datapoints depends on some extra variable. The base map is ‘nzHires’ available in the ‘maps’ package, but any base map should work provided  you set xlim and ylim in the plot.default call to appropriate values.

The input is a .csv dataset where each bird was one row, and there were columns for the bird’s latitude (called googleLat), longitude (called googleLong), degree of external oiling (called percentOil), the day the bird was found coded as numeric days since the start of the response effort (called cumulativeFoundDate), the day the bird was processed coded as numeric days since the start of the response effort (called cumulativeProcessDate), and the day the bird was found coded as dd/mm/yy (called foundDate).

 


read.df <- read.csv(file.choose()) # read in "Birds only with georefs"
library(maps)
library(mapproj)
library(animation)
library(mapdata)
library(fBasics)

ani.record(reset = TRUE) # clear history before recording

for(i in 1:max(read.df$cumulativeProcessDate)) {

subset <- subset(read.df, cumulativeFoundDate==i)
attach(subset)
par(mar=c(5.1,4.1,4.1,1), oma=c(2,1,0,0), bg="white")
date <- as.Date(subset$foundDate[1], "%d/%m")
palette(c("gray",rampPalette(100, name="blue2red")))
plot.default(jitter(googleLat, amount=0.01)~jitter(googleLong, amount=0.01),
             xlim=c(175,179.5), ylim=c(-38.5,-37),
             pch=16,
             col=percentOil+1,
             bg=map("nzHires", xlim=c(175,179.5), ylim=c(-38.5,-37), col="black", lwd=2),
             xlab="longitude",  ylab="latitude")
legend(x=178.2, y=-37.0,
       legend=c("Not Oiled","1% oiled","25% Oiled", "50% Oiled", "75% Oiled", "100% Oiled"),
       pch=16, col=(c("gray",2,25,50,75,100)),
)
mtext(format(date, "%b %d"), SOUTH<-1, cex=1, outer=TRUE)
ani.record() # record the current frame
detach(subset)
}

oopts = ani.options(interval = 0.5)
ani.replay() # replays the animation in the R graphics window

saveGIF({
ani.replay()},
movie.name = "percentOilAsColours.gif", interval = 0.2)

Advertisements

From → Code

One Comment

Trackbacks & Pingbacks

  1. [Ashmore Trip 1] Day 11: Filling the gaps |

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: