# This is the main R script to generate the interpolated maps for electorates and urban areas.
# 
# Required: Run make_seat_folders.R beforehand to generate subfolders in the images folder to store all
# the pictures.
# 
# Optional: If you have downloaded all the electoral boundary shapefiles from the AEC and ABS, then run
# the various "_shp" R files beforehand (do copy_seat_shp_back.R *last*) so that this script will
# use the boundaries defined in the shapefiles where possible.
#
# The main loop starts on line 252; everything before then is loading and preparing data, setting
# parameters, etc.  By default it will do a full run over all seats and urban areas; you can override
# that default by re-defining the seats and urban_areas arrays -- see the commented-out examples
# just before the start of the main loop.
#
# The whole thing takes about 6 hours to run on my computer.

library(sp)
library(spatstat)
library(rgdal)
library(maptools)
library(plyr)
library(ggplot2)
library(gstat)
library(rgeos)
library(scales)
library(fields)
library(EBImage)

# Auxiliary functions:
source("interpolate_seats_aux.R")

# Load the polygon holding the boundaries of Australia's land:
aus_land <- readOGR("shapefiles", "aust_simpl")

max_samples <- 10

# Do we want to use declustering weights?
decluster <- T

# Suburb list:
suburbs <- readOGR("shapefiles", "state_suburbs_points")
suburbs$SSC_NAME <- sub(" \\(.*", "", suburbs$SSC_NAME)

bad_suburbs <- c("Unclassified", "No usual address", "Migratory - Offshore - Shipping")
suburbs <- suburbs[-which(suburbs$SSC_NAME %in% bad_suburbs), ]

# Major town list:
major_towns.data <- read.csv("major_towns_sscnames.csv", header=T, sep=",")
# Convert factors to strings:
i <- sapply(major_towns.data, is.factor)
major_towns.data[i] <- lapply(major_towns.data[i], as.character)

indices_major_towns <- which(suburbs$SSC_NAME %in% major_towns.data$SSC_NAME)

# Make a second column which is mostly the same as SSC_NAME but which over-writes with major town
# names where necessary:
suburbs$Town <- suburbs$SSC_NAME
for (i in indices_major_towns) {
  suburbs$Town[i] <- major_towns.data$Town[which(major_towns.data$SSC_NAME == suburbs$SSC_NAME[i])]
}

# Booths to get rid of because they're a long way from the rest of the electorate:
outlier_booths <- c("Jervis Bay", "Wreck Bay", "Lord Howe Island", "Home Island", "West Island", "Cocos Islands PPVC", "Christmas Island")

# Following variable will decide the size of the grid cells:
num_grid_cells <- 100

# Threshold value (metres); if the grid mesh is greater than this, always select a major town to display
# in a given lattice cell of the map if possible.
threshold_large_cells <- 2000

# There are often many more suburbs than can reasonably be labelled on a plot.  The grid will be split into
# a num_places_y x num_places_x lattice, with one location randomly sampled from each lattice cell.
num_places_x <- 3
num_places_y <- 6

# Load the non-classic divisions:
nonclassic.data <- read.csv("booths_data/nonclassic_divs.csv", header=T, sep=",")
nonclassic_seats <- unique(nonclassic.data$Seat)

# Load the booth data:
years <- c(1993, 1996, 1998, 2001, 2004, 2007, 2010, 2013)
# years <- 2013

TCP_data_colnames <- c("State", "Seat", "Booth", "Lat", "Long", "Formal_votes", "TCP_1", "TCP_2")

# The data takes a bit of time to load, so only load it if the data isn't already loaded:
if (!exists("seatbooth.data")) {
  for (cal_year in years) {
    print(sprintf("Loading data, %d", cal_year))
    seatboothfile_str <- sprintf("booths_data/%d.csv", cal_year)
    boothfile_str <- sprintf("booths_summing_dups/%d.csv", cal_year)
    if (cal_year == years[1]) {
      seatbooth.data <- read.csv(file=seatboothfile_str, header=T, sep=",")
      
      # Extract the TCP columms:
      TCP_seatbooth.data <- seatbooth.data[TCP_data_colnames]
      # Remove the TCP columns from the main data frame:
      seatbooth.data <- seatbooth.data[, !(names(seatbooth.data) %in% c("TCP_1", "TCP_2"))]
      
      seatbooth.data$Year <- cal_year
      TCP_seatbooth.data$Year <- cal_year
      
      booth.data <- read.csv(file=boothfile_str, header=T, sep=",")
      booth.data$Year <- cal_year
    } else {
      temp.data <- read.csv(file=seatboothfile_str, header=T, sep=",")
      
      # Extract the TCP columms:
      temp_TCP_seatbooth.data <- temp.data[, TCP_data_colnames]
      # Remove the TCP columns from the main data frame:
      temp.data <- temp.data[, !(names(temp.data) %in% c("TCP_1", "TCP_2"))]
      
      temp.data$Year <- cal_year
      temp_TCP_seatbooth.data$Year <- cal_year
      
      seatbooth.data <- merge(seatbooth.data, temp.data, all=T)
      TCP_seatbooth.data <- merge(TCP_seatbooth.data, temp_TCP_seatbooth.data, all=T)
      rm(temp.data)
      rm(temp_TCP_seatbooth.data)
      
      temp.data <- read.csv(file=boothfile_str, header=T, sep=",")
      temp.data$Year <- cal_year
      booth.data <- merge(booth.data, temp.data, all=T)
    }
  }
  rm(temp.data)
  
  # Convert factors to strings:
  i <- sapply(seatbooth.data, is.factor)
  seatbooth.data[i] <- lapply(seatbooth.data[i], as.character)
  
  i <- sapply(booth.data, is.factor)
  booth.data[i] <- lapply(booth.data[i], as.character)
  
  i <- sapply(TCP_seatbooth.data, is.factor)
  TCP_seatbooth.data[i] <- lapply(TCP_seatbooth.data[i], as.character)
  
  # Remove points where longitude = 0 (= latitude):
  seatbooth.data <- seatbooth.data[-which(seatbooth.data$Long == 0), ]
  booth.data <- booth.data[-which(booth.data$Long == 0), ]
  TCP_seatbooth.data <- TCP_seatbooth.data[-which(TCP_seatbooth.data$Long == 0), ]
  
  # Remove the outlier booths:
  seatbooth.data <- seatbooth.data[-which(seatbooth.data$Booth %in% outlier_booths), ]
  booth.data <- booth.data[-which(booth.data$Booth %in% outlier_booths), ]
  TCP_seatbooth.data <- TCP_seatbooth.data[-which(TCP_seatbooth.data$Booth %in% outlier_booths), ]
  
  # For seatbooth's, remove points with no formal votes:
  seatbooth.data <- seatbooth.data[-which(seatbooth.data$Formal_votes == 0), ]
  
  # For seatbooth's, remove PPVC's:
  seatbooth.data <- seatbooth.data[-grep(pattern="PPVC", seatbooth.data$Booth), ]
  
  # For TCP_seatbooth's, remove rows from classic divisions:
  TCP_seatbooth.data <- TCP_seatbooth.data[-which(is.na(TCP_seatbooth.data$TCP_1)), ]
  
  # For booth tallies, remove points with less than a threshold number of votes
  # (Want to keep the small but real booths in the seatbooth.data frame because they'll be used
  # to estimate the electorate boundaries):
  threshold_votes <- 40
  booth.data <- booth.data[-which(booth.data$Formal_votes < threshold_votes), ]
  
  # Remove booths with less than a threshold number of formal votes unless the booth is sufficiently isolated:
  threshold_votes <- 100
  threshold_dist <- 3
  
  states <- unique(booth.data$State)
  
  small_nonisolated_booths <- numeric()
  
  for (cal_year in years) {
    for (state in states) {
      indices_year <- intersect(which(booth.data$Year == cal_year), which(booth.data$State == state))
      indices_small <- intersect(indices_year, which(booth.data$Formal_votes < threshold_votes))
      
      coords_small <- cbind(booth.data$Long[indices_small], booth.data$Lat[indices_small])
      coords_all <- cbind(booth.data$Long[indices_year], booth.data$Lat[indices_year])
      
      distances <- rdist.earth(coords_small, coords_all, miles=F)
      if (length(distances) > 0) {
        nearest_dist <- apply(distances, 1, function(x) sort(x)[2])
        small_nonisolated_booths <- c(small_nonisolated_booths, indices_small[which(nearest_dist < threshold_dist)])
      }
    }
  }
  
  booth.data <- booth.data[-small_nonisolated_booths, ]
  
  # Remove hyphens from "Kingsford-Smith":
  seatbooth.data$Seat[which(seatbooth.data$Seat == "Kingsford-Smith")] <- "Kingsford Smith"
  booth.data$Seat[which(booth.data$Seat == "Kingsford-Smith")] <- "Kingsford Smith"
  
  # Rename Prospect to McMahon:
  seatbooth.data$Seat[which(seatbooth.data$Seat == "Prospect")] <- "McMahon"
  booth.data$Seat[which(booth.data$Seat == "Prospect")] <- "McMahon"
  
  
  # Give booth.data and seatbooth.data some spatial properties:
  coordinates(seatbooth.data) <- ~Long+Lat
  proj4string(seatbooth.data) <- CRS("+proj=longlat +datum=WGS84")
  
  coordinates(booth.data) <- ~Long+Lat
  proj4string(booth.data) <- CRS("+proj=longlat +datum=WGS84")
  
  coordinates(TCP_seatbooth.data) <- ~Long+Lat
  proj4string(TCP_seatbooth.data) <- CRS("+proj=longlat +datum=WGS84")
}


seats <- unique(seatbooth.data$Seat)

urban_areas <- c("Brisbane", "Sydney", "Melbourne", "Hobart", "Adelaide", "Perth", "Canberra", "Gold Coast", "Newcastle")
urban_states <- c("QLD", "NSW", "VIC", "TAS", "SA", "WA", "ACT", "QLD", "NSW")

# The script assumes that the folders for individual seat images have already been created, but that the list of
# wanted urban areas might change, so the latter are created here.  The R function dir.create will only generate
# a warning if the folder already exists.
for (urban_area in urban_areas) {
  sub_folder <- sprintf("%s_urban", urban_area)
  dir.create(file.path("./images", sub_folder))
}

# Get columns to interpolate
vote_cols <- grep("^ALP_2PP$", names(seatbooth.data))
vote_cols <- c(vote_cols, grep("_prim", names(seatbooth.data)))
# Tack on a flag for TCP at the end:
vote_cols <- c(vote_cols, 999)

# The following are booths that are used by many electorates.  For some reason WA was full of them.
capcity_booths <- c("Brisbane City", "Sydney", "Melbourne", "Adelaide", "Perth")
capcity_wa1 <- c("Bunbury Central", "Busselton", "Busselton West", "Dunsborough", "Falcon", "Lancelin")
capcity_wa2 <- c("Mandurah Central", "Margaret River", "Nannup", "Rottnest", "Yallingup", "Broome")
capcity_wa3 <- c("Carnarvon", "Exmouth", "Kalbarri", "Divisional Office (PREPOLL)")
capcity_booths <- c(capcity_booths, capcity_wa1, capcity_wa2, capcity_wa3)

# If we only want one seat, say:
# seats <- "Griffith"
# urban_areas <- numeric()

start_index <- 1
end_index <- length(seats) + length(urban_areas)

# If a run was interrupted and we want to resume from where we were up to:
# start_seat <- "Watson"
# start_index <- which(seats==start_seat)

# end_seat <- "Hume"
# end_index <- which(seats==end_seat)


for(i in start_index:end_index) {
  
  if (i <= length(seats)) {
    seat <- seats[i]
    in_seat <- T
  } else {
    urban_i <- i - length(seats)
    # seat <- urban_areas[urban_i]
    seat <- urban_areas[urban_i]
    urban_area <- urban_areas[urban_i]
    this_state <- urban_states[urban_i]
    in_seat <- F
    if (i == length(seats) + 1) {
      # Double the resolution of the grid when we start the urban areas:
      num_grid_cells <- 2*num_grid_cells
    }
  }
  
  print(seat)
  
  if (in_seat) {
   indices_seat_seatbooth <- which(seatbooth.data$Seat == seat)
    
    this_state <- seatbooth.data$State[floor(median(indices_seat_seatbooth))]
    
    # Get just the rows of the data frame from the current state:
    indices_state_seatbooth <- which(seatbooth.data$State == this_state)
    
    state_seatbooth.data = seatbooth.data[indices_state_seatbooth, ]
        
    indices_capcity <- which(state_seatbooth.data$Booth %in% capcity_booths)
    
    x <- state_seatbooth.data$Long[indices_state_seatbooth]
    y <- state_seatbooth.data$Lat[indices_state_seatbooth]
    
    # Re-find the indices of the current seat, in the new data frame restricted
    # to the state
    indices_seat_seatbooth <- which(state_seatbooth.data$Seat == seat)
    
    # Need to remove the capital city booth if it's outside a bounding box defined by the remaining
    # booths in the seat.  This is repeated for each year.  Only once this cleaning has been done is the
    # common grid for the interpolations created.
    for (cal_year in years) {
      indices_year <- which(state_seatbooth.data$Year == cal_year)
      indices_seat_year <- intersect(indices_year, indices_seat_seatbooth)
      
      capcity_indices <- intersect(indices_seat_year, indices_capcity)
      
      if (length(capcity_indices) > 0) {
        indices_excl_capcity <- indices_seat_year[which(!(indices_seat_year %in% capcity_indices))]
        min_x <- min(state_seatbooth.data$Long[indices_excl_capcity])
        max_x <- max(state_seatbooth.data$Long[indices_excl_capcity])
        min_y <- min(state_seatbooth.data$Lat[indices_excl_capcity])
        max_y <- max(state_seatbooth.data$Lat[indices_excl_capcity])
        
        for (i in 1:length(capcity_indices)) {
          capcity_x <- state_seatbooth.data$Long[capcity_indices[i]]
          capcity_y <- state_seatbooth.data$Lat[capcity_indices[i]]
          
          if ((capcity_x > max_x) || (capcity_x < min_x) || (capcity_y > max_y) || (capcity_y < min_y)) {
            indices_seat_seatbooth <- indices_seat_seatbooth[which(!indices_seat_seatbooth==capcity_indices[i])]
            indices_seat_year <- indices_seat_year[which(!indices_seat_year==capcity_indices[i])]
          }
        }
      }
    }
    
    # Now create the grid that is common to all years for this electorate:
    median_long <- median(state_seatbooth.data$Long[indices_seat_seatbooth])
  } else {
    # With the urban areas, need only load the polygon and get the bbox.
    poly_folder <- sprintf("shapefiles/urban_areas/%s", urban_area)
    urban_poly <- readOGR(poly_folder, urban_area)
    median_long <- mean(urban_poly@bbox[1, ])
  }
  
  # ********************************************************************
  # Following is common for both seats and urban areas: 
  zone <- floor((186+median_long)/6)
  CRS_string <- sprintf("+proj=utm +south +zone=%d +datum=WGS84", zone)
  
  suburbs.proj <- spTransform(suburbs,  CRS(CRS_string))
  coords_suburbs <- as.data.frame(coordinates(suburbs.proj))
  
  indices_state_booth <- which(booth.data$State == this_state)
  state_booth.data = booth.data[indices_state_booth, ]
  state_booth.data.proj <- spTransform(state_booth.data, CRS(CRS_string))
  coords_booth <- as.data.frame(coordinates(state_booth.data.proj))
  easting_booth <- coords_booth$Long
  northing_booth <- coords_booth$Lat
  # ********************************************************************
  
  if (in_seat) {
    # While we're here, transform the TCP data frame:
    TCP_seatbooth.data.proj <- spTransform(TCP_seatbooth.data,  CRS(CRS_string))
    
    state_seatbooth.data.proj <- spTransform(state_seatbooth.data,  CRS(CRS_string))
    coords <- as.data.frame(coordinates(state_seatbooth.data.proj))
    easting <- coords$Long
    northing <- coords$Lat
    
    easting_seat <- easting[indices_seat_seatbooth]
    northing_seat <- northing[indices_seat_seatbooth]
  } else {
    urban_poly.proj <- spTransform(urban_poly, CRS(CRS_string))
    easting_seat <- urban_poly.proj@bbox[1, ]
    northing_seat <- urban_poly.proj@bbox[2, ]
  }
  
  range_x <- max(easting_seat) - min(easting_seat)
  range_y <- max(northing_seat) - min(northing_seat)
  
  # Round the mesh size up to the nearest 50 metres:
  grid_mesh <- round_any(max(c(1.2*range_x, 1.2*range_y)) / num_grid_cells, 50, f = ceiling)
  
  min_x <- round_any(min(easting_seat) - 0.1*range_x, grid_mesh, f = floor)
  max_x <- round_any(max(easting_seat) + 0.1*range_x, grid_mesh, f = ceiling)
  min_y <- round_any(min(northing_seat) - 0.1*range_y, grid_mesh, f = floor)
  max_y <- round_any(max(northing_seat) + 0.1*range_y, grid_mesh, f = ceiling)
  
  # Make a subet of the data frame that includes only booths within the grid area;
  # used later for setting the colour scale.
  if (in_seat) {
    grid_seatbooth.data <- state_seatbooth.data[which((easting > min_x) & (easting < max_x) & (northing > min_y) & (northing < max_y)), ]
  } else {
    indices_booth_urbanbbox <- which((easting_booth > min(easting_seat)) & (easting_booth < max(easting_seat)) & (northing_booth > min(northing_seat)) & (northing_booth < max(northing_seat)))
  }
  grid_booth.data <- state_booth.data[which((easting_booth > min_x) & (easting_booth < max_x) & (northing_booth > min_y) & (northing_booth < max_y)), ]
  
  
  
  
  # Make the grid itself and give it spatial properties:
  east_grid <- seq(min_x, max_x, grid_mesh)
  north_grid <- seq(min_y, max_y, grid_mesh)
  aspect <- length(north_grid) / length(east_grid)
  png_height <- ceiling(min(c(500, aspect*360 + 80)))
  
  grid <- expand.grid(x=east_grid, y=north_grid)
  coordinates(grid) <- ~x+y
  proj4string(grid) <- CRS(CRS_string)
  
  # Sample from the suburbs
  suburb_mesh_x <- (max_x - min_x) / num_places_x
  suburb_mesh_y <- (max_y - min_y) / num_places_y
  east_suburb_lattice <- seq(min_x + suburb_mesh_x/2, max_x - suburb_mesh_x/2, suburb_mesh_x)
  north_suburb_lattice <- seq(min_y + suburb_mesh_y/2, max_y - suburb_mesh_y/2, suburb_mesh_y)
  suburb_sample_grid <- expand.grid(x=east_suburb_lattice, y=north_suburb_lattice)
  
  suburb_sample_grid$done <- 0
  
  indices_sampled_suburbs <- numeric()
  
  # First run through, pick out any major towns:
  if (grid_mesh > threshold_large_cells) {
    for (ct_lattice in sample(1:dim(suburb_sample_grid)[1])) {
      x_low <- suburb_sample_grid$x[ct_lattice] - suburb_mesh_x/2
      x_high <- suburb_sample_grid$x[ct_lattice] + suburb_mesh_x/2
      y_low <- suburb_sample_grid$y[ct_lattice] - suburb_mesh_y/2
      y_high <- suburb_sample_grid$y[ct_lattice] + suburb_mesh_y/2
      
      indices_lattice_cell_suburbs <- which((coords_suburbs$coords.x1 > x_low) & (coords_suburbs$coords.x1 < x_high) & (coords_suburbs$coords.x2 > y_low) & (coords_suburbs$coords.x2 < y_high))
      indices_lattice_cell_majortowns <- intersect(indices_lattice_cell_suburbs, indices_major_towns)
      
      if (length(indices_lattice_cell_majortowns) > 0) {
        test_sample <- indices_lattice_cell_majortowns[1]
        indices_sampled_suburbs <- c(indices_sampled_suburbs, test_sample)
        suburb_sample_grid$done[ct_lattice] <- 1
      }
    }
  }
  
  for (ct_lattice in sample(which(suburb_sample_grid$done == 0))) {
    x_low <- suburb_sample_grid$x[ct_lattice] - suburb_mesh_x/2
    x_high <- suburb_sample_grid$x[ct_lattice] + suburb_mesh_x/2
    y_low <- suburb_sample_grid$y[ct_lattice] - suburb_mesh_y/2
    y_high <- suburb_sample_grid$y[ct_lattice] + suburb_mesh_y/2
    
    indices_lattice_cell_suburbs <- which((coords_suburbs$coords.x1 > x_low) & (coords_suburbs$coords.x1 < x_high) & (coords_suburbs$coords.x2 > y_low) & (coords_suburbs$coords.x2 < y_high))
    
    # Want to try to ensure that the place names don't get written over the top of one another.
    # To that end, as we tour randomly the lattice grid and sample one placename from each
    # lattice cell, check that the latest sample isn't within 3 grid cells of a previously-
    # sampled placename.  (Should make the threshold of 3 dependent on the num_grid_cells).
    got_sample <- 0
    
    # Urban area plots have twice the resolution, so need to double the number of grid cells 
    # required to separate the place names.
    if (in_seat) {
      mult <- 1
    } else {
      mult <- 2
    }
    
    while (got_sample == 0) {
      if (length(indices_lattice_cell_suburbs) > 0) {
        test_sample <- sample(indices_lattice_cell_suburbs, 1)
        if (sum(which(abs(coords_suburbs$coords.x2[indices_sampled_suburbs] - coords_suburbs$coords.x2[test_sample]) < 3*grid_mesh*mult)) > 0) {
          indices_lattice_cell_suburbs <- indices_lattice_cell_suburbs[-which(indices_lattice_cell_suburbs==test_sample)]
        } else {
          got_sample = 1
          indices_sampled_suburbs <- c(indices_sampled_suburbs, test_sample)
        }
      } else {
        # Don't return a sample, but say we do anyway to exit the while loop.
        got_sample <- 1
      }
    }
  }
  
  
  state_poly <- aus_land[aus_land$STATE_NAME == long_state(this_state), ]
  state_poly <- spTransform(state_poly, CRS(CRS_string))
  
  national_poly <- spTransform(aus_land, CRS(CRS_string))
  
  # Want to know if the grid crosses a state boundary:
  grid_corners <- expand.grid(x=c(min_x, max_x), y=c(min_y, max_y))
  coordinates(grid_corners) <- ~x+y
  proj4string(grid_corners) <- CRS(CRS_string)
  
  grid_corner.data <- over(SpatialPoints(coordinates(grid_corners), proj4string=CRS(CRS_string)), national_poly)
  corner_states <- unique(grid_corner.data$STATE_NAME[which(!is.na(grid_corner.data$STATE_NAME))])
  if (length(corner_states) == 0) {
    need_national_poly <- 1
  } else {
    if ((length(corner_states) > 1) | (corner_states[1] != long_state(this_state))) {
      need_national_poly <- 1
    } else {
      need_national_poly <- 0
    }
  }
  
  if (need_national_poly == 0) {
    onland.data <- over(SpatialPoints(coordinates(grid), proj4string=CRS(CRS_string)), state_poly)
    land_cells <- which(!is.na(onland.data$STATE_CODE))
    state_cells <- land_cells
    water_cells <- which(is.na(onland.data$STATE_CODE))
  } else {
    onland.data <- over(SpatialPoints(coordinates(grid), proj4string=CRS(CRS_string)), national_poly)
    land_cells <- which(!is.na(onland.data$STATE_NAME))
    state_cells <- which(onland.data$STATE_NAME == long_state(this_state))
    water_cells <- which(is.na(onland.data$STATE_NAME))
  }
  
  grid_ppp <- ppp(grid$x, grid$y, c(min_x, max_x), c(min_y, max_y))
  
  if (!in_seat) {
    # When working with seats, the seat area will change each year with redistributions.
    # But for urban areas, we only need to calculate the "seat_cells" once.
    in_urban.data <- over(SpatialPoints(coordinates(grid), proj4string=CRS(CRS_string)), urban_poly.proj)
    seat_cells <- intersect(land_cells, which(!is.na(in_urban.data[[names(in_urban.data)[1]]])))
    nonseat_cells <- setdiff(land_cells, seat_cells)
  }
  
  # Plot the basic map, showing land, water, and the sampled place names.
  # We do three times, with different numbers of digits in the colourbar scale labels.
  # 'digits' refers to the number of digits in the upper limit of the scale.
  digits <- c(1, 2, 3)
  land_val <- c(4, 40, -1)
  water_val <- c(6, 60, 101)
  end_px <- as.complex(c(0, 0, 0)) # to be filled in later
  
  plot_suburbs <- suburbs$SSC_NAME[indices_sampled_suburbs]
  if (grid_mesh > threshold_large_cells) {
    plot_suburbs <- suburbs$Town[indices_sampled_suburbs]
  }
  
  colour_vec <- c("#666666", "#87ceeb")
  
  if (in_seat) {
    sub_folder <- seat
  } else {
    sub_folder <- sprintf("%s_urban", seat)
  }
  
  # Following two variables are named identically to how they appear in the plot_map function
  # because of a scoping bug with aes in ggplot2:
  suburb_names <- plot_suburbs
  suburb_xy <- cbind(coords_suburbs$coords.x1[indices_sampled_suburbs], coords_suburbs$coords.x2[indices_sampled_suburbs])
  
  for (ct in 1:length(digits)) {
    grid.df <- basic_map(grid, land_cells, land_val[ct], water_cells, water_val[ct])
    
    colour_val <- rescale(c(land_val[ct], water_val[ct]))
    title_str <- sprintf("%s\n ", seat)
    output_map <- sprintf("images/%s/map%d_%s.png", sub_folder, digits[ct], seat)
    
    end_px[ct] <- plot_map(c(min_x, max_x), c(min_y, max_y), grid.df, colour_vec, colour_val, suburb_xy, suburb_names, title_str, output_map, png_height)
  }
  
  # With the grid set up, interpolate each election's results:
  for (cal_year in years) {
    print(sprintf("Starting year %d", cal_year))
    
    # Make a dataframe for this state in this year:
    indices_boothdf_year <- which(state_booth.data$Year == cal_year)
    state_year.data <- spTransform(state_booth.data[indices_boothdf_year, ],  CRS(CRS_string))
    coords <- as.data.frame(coordinates(state_year.data))
    
    if (in_seat) {
      indices_boothdf_seat <- which(state_year.data$Seat == seat)
      
      # Used to set the size scale for the booth basemap plot:
      min_formalvotes <- min(state_year.data$Formal_votes[indices_boothdf_seat])
      max_formalvotes <- max(state_year.data$Formal_votes[indices_boothdf_seat])
      state_year.data$Votes_toplot <- state_year.data$Formal_votes
      state_year.data$Votes_toplot[which(state_year.data$Votes_toplot > max_formalvotes)] <- max_formalvotes
      state_year.data$Votes_toplot[which(state_year.data$Votes_toplot < min_formalvotes)] <- min_formalvotes
      
      
      # Need to calculate the seat_cells each year when working with seats rather than urban areas.
      if (length(indices_boothdf_seat) > 0) {
        do_interp <- T
        shp_file <- paste("shapefiles/seats/", seat, "/", cal_year, ".shp", sep="")
        if (file.exists(shp_file)) {
          year_str <- sprintf("%d", cal_year)
          seat_folder <- sprintf("shapefiles/seats/%s", seat)
          seat_poly <- readOGR(seat_folder, year_str)
          seat_poly <- spTransform(seat_poly, CRS(CRS_string))
          inseat.data <- over(SpatialPoints(coordinates(grid), proj4string=CRS(CRS_string)), seat_poly)
          seat_cells <- intersect(land_cells, which(!is.na(inseat.data[[names(inseat.data)[1]]])))
          nonseat_cells <- setdiff(land_cells, seat_cells)
        } else {
          booths_ppp <- ppp(coords$Long, coords$Lat, c(min(coords$Long), max(coords$Long)), c(min(coords$Lat), max(coords$Lat)))
          nearest <- nncross(grid_ppp[state_cells, ], booths_ppp, what="which", k=1)
          
          # Assign the nearest booth index to the nearest_booths vector.
          # In cases where the nearest booth is shared, those indices will be over-written in the loop below.
          nearest_booths <- state_year.data$Booth[nearest]
          
          nearest_ismultiple <- which(state_year.data$Seat[nearest] == "MULTIPLE")
          
          kth_nbr = 2
          
          while (length(nearest_ismultiple) > 0) {
            # Calculate kth-nearest neighbours for all grid cells even though we only need a subset,
            # to keep the indexing consistent.
            k_nearest <- nncross(grid_ppp[state_cells, ], booths_ppp, what="which", k=kth_nbr)
            
            k_nearest_booths <- state_year.data$Booth[k_nearest]
            # Update list of nearest neighbours:
            nearest[nearest_ismultiple] <- k_nearest[nearest_ismultiple]
            
            # nearest_booths[nearest_ismultiple] <- k_nearest_booths[nearest_ismultiple]
            nearest_ismultiple <- intersect(nearest_ismultiple, which(state_year.data$Seat[nearest] == "MULTIPLE"))
            
            kth_nbr = kth_nbr + 1
          }
          
          # Need to get rid of it or R spits on the next loop.
          rm(booths_ppp)
          
          seat_cells <- state_cells[which(nearest %in% indices_boothdf_seat)]
          nonseat_cells <- setdiff(land_cells, seat_cells)
        }
      } else {
        do_interp <- F
      }
    } else {
      # Always do the interpolations when we're working with an urban area:
      do_interp <- T
    }
    
    if (do_interp) {
      # Reset the grid data frame to remove any stray interpolated variables from it:
      
      grid.df <- as.data.frame(grid)
      indices_year_booth <- which(state_booth.data$Year == cal_year)
      
      if (in_seat) {
        indices_seatbooth_year <- which(state_seatbooth.data$Year == cal_year)
        indices_seatbooth_seat <- which(state_seatbooth.data$Seat == seat)
        indices_seatbooth_seat_year <- intersect(indices_seatbooth_year, indices_seatbooth_seat)
        
        indices_tcp_year <- which(TCP_seatbooth.data$Year == cal_year)
        indices_tcp_seat <- which(TCP_seatbooth.data$Seat == seat)
        indices_tcp_seat_year <- intersect(indices_tcp_year, indices_tcp_seat)
        
        easting_seatbooth_year <- easting[indices_seatbooth_seat_year]
        northing_seatbooth_year <- northing[indices_seatbooth_seat_year]
        
        # Want to find declustering weights; choose as a window the median nearest-neighbour spacing:
        thisseat_booths_ppp <- ppp(easting_seatbooth_year, northing_seatbooth_year, c(min_x, max_x), c(min_y, max_y))
        median_nearest <- median(nndist(thisseat_booths_ppp))
        
        # Set it up for both the booth dataframe and the seatbooth dataframe:
        state_booth.data.proj$wts <- NA
        state_seatbooth.data.proj$wts <- NA
        TCP_seatbooth.data.proj$wts <- NA
        
        easting_booth_year <- easting_booth[indices_year_booth]
        northing_booth_year <- northing_booth[indices_year_booth]
        state_booths_year_ppp <- ppp(easting_booth_year, northing_booth_year, c(min(easting_booth), max(easting_booth)), c(min(northing_booth), max(northing_booth)))
        
        
        state_booth.data.proj$wts[indices_year_booth] <- decl_wts(state_booths_year_ppp, median_nearest/2)
        state_seatbooth.data.proj$wts[indices_seatbooth_seat_year] <- decl_wts(thisseat_booths_ppp, median_nearest/2)
        
        if (seat %in% nonclassic_seats) {
          TCP_seatbooth.data.proj$wts[indices_tcp_seat_year] <- decl_wts(thisseat_booths_ppp, median_nearest/2)
        }
      } else {
        # indices_booth_urbanbbox
        indices_booth_year_urbanbbox <- intersect(indices_year_booth, indices_booth_urbanbbox)
        urban_booths_ppp <- ppp(easting_booth[indices_booth_year_urbanbbox], northing_booth[indices_booth_year_urbanbbox], c(min_x, max_x), c(min_y, max_y))
        median_nearest <- median(nndist(urban_booths_ppp))
        
        easting_booth_year <- easting_booth[indices_year_booth]
        northing_booth_year <- northing_booth[indices_year_booth]
        state_booths_year_ppp <- ppp(easting_booth_year, northing_booth_year, c(min(easting_booth), max(easting_booth)), c(min(northing_booth), max(northing_booth)))
        
        state_booth.data.proj$wts <- NA
        state_booth.data.proj$wts[indices_year_booth] <- decl_wts(state_booths_year_ppp, median_nearest/2)
      }
      
      
      for (vote_col in vote_cols) {
        if (vote_col < 999) {
          col_name <- names(state_booth.data)[vote_col]
          indices_votedefined <- which(!is.na(state_booth.data[[col_name]]))
          indices_inputdata <- intersect(indices_year_booth, indices_votedefined)
          
          if (in_seat) {
            # Above indices are from the booth data.  Need to check from the seatbooth.data if the party
            # actually contested the seat this year.
            
            indices_seatbooth_votedefined <- which(!is.na(state_seatbooth.data[[col_name]]))
            indices_party_contested <- intersect(indices_seatbooth_seat_year, indices_seatbooth_votedefined)
          } else {
            # For the urban areas, only need that the party contested somewhere in the state this year, and that
            # the current col_name isn't seat-specific.
            
            indices_party_contested <- indices_inputdata
            if ((col_name == "IND_prim") | (col_name == "TCP_1")) {
              indices_party_contested <- numeric()
            }
          }
          
          
          # Get min and max party vote percentages across all years for this grid area, to set colour scale:
          
          if (col_name == "IND_prim") {
            if (in_seat) {
              indices_seatbooth_seat_ind <- intersect(indices_seatbooth_seat, indices_seatbooth_votedefined)
              min_vote <- min(state_seatbooth.data[[col_name]][indices_seatbooth_seat_ind])
              max_vote <- max(state_seatbooth.data[[col_name]][indices_seatbooth_seat_ind])
            }
          } else {
            if (in_seat) {
              indices_grid_seatbooth_party <- which(!is.na(grid_seatbooth.data[[col_name]]))
              min_vote <- min(grid_seatbooth.data[[col_name]][indices_grid_seatbooth_party])
              max_vote <- max(grid_seatbooth.data[[col_name]][indices_grid_seatbooth_party])
            } else {
              indices_grid_booth_party <- which(!is.na(grid_booth.data[[col_name]]))
              
              if (length(indices_grid_booth_party) == 0) {
                # Don't interpolate.
                indices_party_contested = numeric()
              }
              
              min_vote <- min(grid_booth.data[[col_name]][indices_grid_booth_party])
              max_vote <- max(grid_booth.data[[col_name]][indices_grid_booth_party])
            }
          }
        } else {
          # TCP
          
          indices_party_contested <- numeric()
          
          # Over-write indices_party_contested if we do actually want to plot the TCP:
          if ((seat %in% nonclassic_seats) & (in_seat)) {
            indices_tcp_year <- which(TCP_seatbooth.data$Year == cal_year)
            indices_tcp_seat <- which(TCP_seatbooth.data$Seat == seat)
            indices_tcp_seat_year <- intersect(indices_tcp_year, indices_tcp_seat)
            if (length(indices_tcp_seat_year) > 0) {
              indices_party_contested <- indices_tcp_seat_year
              col_name <- "TCP_1"
              
              indices_nonclassic_year <- which(nonclassic.data$Year == cal_year)
              indices_nonclassic_seat <- which(nonclassic.data$Seat == seat)
              index_nonclassic_seat_year <- intersect(indices_nonclassic_year, indices_nonclassic_seat)
              
              TCP_party1 <- toString(nonclassic.data$TCP_1[index_nonclassic_seat_year])
              TCP_party2 <- toString(nonclassic.data$TCP_2[index_nonclassic_seat_year])
              
              min_vote <- min(TCP_seatbooth.data[[col_name]][indices_tcp_seat_year])
              max_vote <- max(TCP_seatbooth.data[[col_name]][indices_tcp_seat_year])
            }
          }
        }
        
        # Bit of an ugly way to try to avoid grid cells outside min-max range:
        min_vote <- max(c(min_vote - 1, 0))
        max_vote <- min(c(max_vote + 1, 100))
        
        if (length(indices_party_contested) > 0) {
          
          if (decluster) {
            if (col_name == "TCP_1") {
              TCP_seatbooth.data.proj$var_wts <- TCP_seatbooth.data.proj$TCP_1 * TCP_seatbooth.data.proj$wts
              interp_idw_numer.df <- as.data.frame(krige(var_wts ~ 1, TCP_seatbooth.data.proj[indices_party_contested, ], grid[seat_cells, ], nmax=max_samples))
              interp_idw_denom.df <- as.data.frame(krige(wts ~ 1, TCP_seatbooth.data.proj[indices_party_contested, ], grid[seat_cells, ], nmax=max_samples))
            } else if (col_name == "IND_prim") {
              state_seatbooth.data.proj$var_wts <- state_seatbooth.data.proj$IND_prim * state_seatbooth.data.proj$wts
              interp_idw_numer.df <- as.data.frame(krige(var_wts ~ 1, state_seatbooth.data.proj[indices_party_contested, ], grid[seat_cells, ], nmax=max_samples))
              interp_idw_denom.df <- as.data.frame(krige(wts ~ 1, state_seatbooth.data.proj[indices_party_contested, ], grid[seat_cells, ], nmax=max_samples))
            } else {
              state_booth.data.proj$var_wts <- state_booth.data.proj[[col_name]] * state_booth.data.proj$wts
              interp_idw_numer.df <- as.data.frame(krige(var_wts ~ 1, state_booth.data.proj[indices_inputdata, ], grid[seat_cells, ], nmax=max_samples))
              interp_idw_denom.df <- as.data.frame(krige(wts ~ 1, state_booth.data.proj[indices_inputdata, ], grid[seat_cells, ], nmax=max_samples))
            }
            interp_idw.df <- interp_idw_numer.df
            interp_idw.df$var1.pred <- interp_idw_numer.df$var1.pred / interp_idw_denom.df$var1.pred
          } else {
            interp_formula <- as.formula(sprintf("%s ~ 1", col_name))
            if (col_name == "TCP_1") {
              # Input data comes from the TCP data frame:
              interp_idw <- krige(interp_formula, TCP_seatbooth.data.proj[indices_party_contested, ], grid[seat_cells, ])
            } else if (col_name == "IND_prim") {
              interp_idw <- krige(interp_formula, state_seatbooth.data.proj[indices_party_contested, ], grid[seat_cells, ])
            } else {
              interp_idw <- krige(interp_formula, state_booth.data.proj[indices_inputdata, ], grid[seat_cells, ], nmax=10)
            }
            interp_idw.df <- as.data.frame(interp_idw)
          }
          
          # Reduce any interpolated cells above the max_vote to that threshold (ugly):
          interp_idw.df$var1.pred[which(interp_idw.df$var1.pred > 0.999*max_vote)] <- 0.999*max_vote
          
          # Following is an ugly hack to make min_vote positive
          if (min_vote == 0.0) {
            min_vote <- 0.01
            bad_cells <- which(interp_idw.df$var1.pred < min_vote + 0.1)
            interp_idw.df$var1.pred[bad_cells] <- min_vote + 0.1
          }
          
          # Resulting interpolation is only defined on some cells.
          # Try to re-write the following silly thing with a merge, all=T
          interp_idw.df$xy <- complex(real=interp_idw.df$x, imaginary=interp_idw.df$y)
          
          grid.df$xy <- complex(real=grid.df$x, imaginary=grid.df$y)
          grid.df$interp <- NA
          
          # Get indices where coords of the cells estimated into interp_idw.df match the grid.df:
          grid.df$interp[match(interp_idw.df$xy, grid.df$xy)] <- interp_idw.df$var1.pred
          
          if (length(nonseat_cells) > 0) {
            if ((col_name == "ALP_2PP") | (col_name == "TCP_1")) {
              grid.df$interp[nonseat_cells] <- -1
            } else {
              grid.df$interp[nonseat_cells] <- min_vote - 0.001
            }
          }
          
          if (length(water_cells) > 0) {
            if ((col_name == "ALP_2PP") | (col_name == "TCP_1")) {
              grid.df$interp[water_cells] <- 101
            } else {
              grid.df$interp[water_cells] <- max_vote + 0.01
            }
            
          }
          
          interp_raster <- ggplot() + scale_x_continuous(limits=c(min_x, max_x), label=function(x) format(x/1000)) + scale_y_continuous(limits=c(min_y, max_y), label=function(x) format(x/1000)) + coord_equal()
          interp_raster <- interp_raster + theme(axis.title.y = element_text(vjust=0.4))
          interp_raster_basic <- interp_raster
          
          # Plot a raster of values at 101 to trick the colour scale into going that high even when there's
          # no water.  Then plot the proper raster over the top of it:
          if ((col_name == "ALP_2PP") | (col_name == "TCP_1")) {
            grid.df$fakewater <- 101
            grid.df$fakenonseat <- -1
          } else {
            grid.df$fakewater <- max_vote + 0.01
            grid.df$fakenonseat <- min_vote - 0.001
          }
          
          water_layer <- geom_raster(data=grid.df, aes(x=x, y=y, fill=fakewater))
          nonseat_layer <- geom_raster(data=grid.df, aes(x=x, y=y, fill=fakenonseat))
          main_layer <- geom_raster(data=grid.df, aes(x=x, y=y, fill=interp))
          
          interp_raster <- interp_raster + water_layer + nonseat_layer + main_layer
          
          if (col_name == "ALP_2PP") {
            colour_vec <- c("#666666", "#0000FF", "#0000FF", "#FFFFFF", "#FF0000", "#FF0000", "#87ceeb")
            colour_val <- rescale(c(-1, 0, 20, 50, 80, 100, 101))
          } else if (col_name == "TCP_1") {
            # Need to append "_prim" to the two-candidate parties to pass the string to party_colours:
            p1 <- sprintf("%s_prim", TCP_party1)
            p2 <- sprintf("%s_prim", TCP_party2)
            
            colours_touse <- c(party_colours(p2), party_colours(p1))[c(2, 4)]
            if ((colours_touse[1] == "#444444") | (colours_touse[2] == "#444444")) {
              nonseat_colour <- "#4A664A"
            } else {
              nonseat_colour <- "#666666"
            }
            
            mid_colour <- "#FFFFFF"
            if ("#FFFF00" %in% colours_touse) {
              # Don't want to have a scale go yellow-white
              mid_colour <- "#888888"
              nonseat_colour <- "#4A664A"
            }
            
            colour_vec <- c(nonseat_colour, colours_touse[1], colours_touse[1], mid_colour, colours_touse[2], colours_touse[2], "#87ceeb")
            colour_val <- rescale(c(-1, 0, 20, 50, 80, 100, 101))
          } else {
            colours_touse <- party_colours(col_name)
            if (colours_touse[2] == "#444444") {
              # Make the non-seat land cells slightly greener rather than grey if the colour scale goes
              # to grey.
              nonseat_colour <- "#4A664A"
            } else {
              nonseat_colour <- "#666666"
            }
            colour_vec <- c(nonseat_colour, colours_touse, "#87ceeb")
            colour_val <- rescale(c(min_vote - 0.001, min_vote, max_vote, max_vote + 0.01))
          }
          
          interp_raster <- interp_raster + scale_fill_gradientn(colours=colour_vec, space="rgb", values=colour_val)
          interp_raster <- interp_raster + guides(fill = guide_colorbar(title=""))
          
          if (in_seat) {
            # Plot booth locations:
            # interp_raster <- interp_raster + geom_point(aes(coords$Long[indices_boothdf_seat], coords$Lat[indices_boothdf_seat]), colour="#000000")
            booth_clr <- booth_colour(colour_vec)
            interp_raster <- interp_raster + geom_point(aes(coords$Long, coords$Lat, size=state_year.data$Votes_toplot), shape="+", colour=booth_clr, alpha=0.35, range=c(3, 6)) + guides(size=F)
          }
          
          if (col_name == "TCP_1") {
            title_str <- sprintf("%s, %s\nTCP: %s v %s", seat, cal_year, TCP_party1, TCP_party2)
          } else {
            title_str <- sprintf("%s, %s\n%s", seat, cal_year, col_name)
          }
          
          interp_raster <- interp_raster + ggtitle(title_str)
          
          if (in_seat) {
            sub_folder <- seat
          } else {
            sub_folder <- sprintf("%s_urban", seat)
          }
          
          if ((col_name != "ALP_2PP") & (col_name != "TCP_1")) {
            output_map_colour <- sprintf("images/%s/%s_%s_colour.png", sub_folder, seat, col_name)
            
            if (!file.exists(output_map_colour)) {
              # During the first year we have this party in this seat, plot a dummy raster plot to file with a colour-scale
              # that's guaranteed not to be broken.  We'll later use EBImage to paste the proper
              # colour scale over ggplot2's default, which occasionally strays into the 
              # fakewater or fakenonseat colours.
              grid.df$temp_var <- NA
              grid.df$temp_var[1] <- min_vote
              grid.df$temp_var[2] <- max_vote
              
              colour_vec_temp <- colours_touse
              colour_val_temp <- rescale(c(min_vote, max_vote))
              
              colourscale_layer <- geom_raster(data=grid.df, aes(x=x, y=y, fill=temp_var))
              
              interp_raster_basic <- interp_raster_basic + colourscale_layer
              interp_raster_basic <- interp_raster_basic + scale_fill_gradientn(colours=colour_vec_temp, space="rgb", values=colour_val_temp)
              interp_raster_basic <- interp_raster_basic + guides(fill = guide_colorbar(title=""))
              interp_raster_basic <- interp_raster_basic + ggtitle(title_str)
              
              output_map_colour <- sprintf("images/%s/%s_%s_colour.png", sub_folder, seat, col_name)
              
              png_print(output_map_colour, 500, png_height, interp_raster_basic)
            }
          }
          
          
          if ((col_name == "ALP_2PP") & (seat %in% nonclassic_seats) & (in_seat)) {
            indices_tcp_year <- which(TCP_seatbooth.data$Year == cal_year)
            indices_tcp_seat <- which(TCP_seatbooth.data$Seat == seat)
            indices_tcp_seat_year <- intersect(indices_tcp_year, indices_tcp_seat)
            if (length(indices_tcp_seat_year) == 0) {
              # Plot the TPP raster as the TCP raster:
              title_str <- sprintf("%s, %s\nTCP: ALP v LNP", seat, cal_year)
              tcp_interp_raster <- interp_raster + ggtitle(title_str)
              output_map <- sprintf("images/%s/%s_TCP_1_%d.png", sub_folder, seat, cal_year)
              png_print(output_map, 500, png_height, tcp_interp_raster)
            }
          }
          
          output_map <- sprintf("images/%s/%s_%s_%d.png", sub_folder, seat, col_name, cal_year)
          png_print(output_map, 500, png_height, interp_raster)
          
          # Paste the proper colour scale over the possibly-broken one that ggplot2 just made:
          if ((col_name != "ALP_2PP") & (col_name != "TCP_1")) {
            output_map_colour <- sprintf("images/%s/%s_%s_colour.png", sub_folder, seat, col_name)
            colourscale_img <- readImage(output_map_colour)
            main_img <- readImage(output_map)
            
            # Need to do two rows because one might be one of the white grid lines (precise row names assume 500x500px image):            
            real_end <- raster_corner(main_img)
            real_end_r <- Re(real_end)
            
            main_img[real_end_r:500, , ] <- colourscale_img[real_end_r:500, , ]
            writeImage(main_img, output_map)
            
            # Check that the dimensions match one of the already-existing maps; if not, make a new one.
            map_digits <- which(end_px == real_end)
            if (length(map_digits) == 0) {
              new_digit <- length(digits) + 1
              digits <- c(digits, new_digit)
              
              colour_vec <- c("#666666", "#87ceeb")
              colour_val <- rescale(c(min_vote - 0.001, max_vote + 0.01))
              title_str <- sprintf("%s\n ", seat)
              output_map <- sprintf("images/%s/map%d_%s.png", sub_folder, new_digit, seat)
              
              grid.df.temp <- basic_map(grid, land_cells, min_vote - 0.001, water_cells, max_vote + 0.01)
              end_px[new_digit] <- plot_map(c(min_x, max_x), c(min_y, max_y), grid.df.temp, colour_vec, colour_val, suburb_xy, suburb_names, title_str, output_map, png_height)
            }
          }
          
          # For the TPP, make another raster with a moving colour-scale.
          if (col_name == "ALP_2PP") {
            if (in_seat) {
              ALP_2PP_votes <- state_seatbooth.data$ALP_2PP[indices_seatbooth_seat_year]
              formal_votes <- state_seatbooth.data$Formal_votes[indices_seatbooth_seat_year]
            } else {
              indices_grid_booth_year <- intersect(which(grid_booth.data$Year == cal_year), indices_grid_booth_party)
              
              ALP_2PP_votes <- grid_booth.data$ALP_2PP[indices_grid_booth_year]
              formal_votes <- grid_booth.data$Formal_votes[indices_grid_booth_year]
            }
            mean_2PP <- sum(ALP_2PP_votes * formal_votes) / sum(formal_votes)
            
            scale_low <- max(0, mean_2PP - 30)
            scale_high <- min(100, mean_2PP + 30)
            
            colour_vec <- c("#666666", "#0000FF", "#0000FF", "#FFFFFF", "#FF0000", "#FF0000", "#87ceeb")
            colour_val <- rescale(c(-1, 0, scale_low, mean_2PP, scale_high, 100, 101))
            
            title_str <- sprintf("%s, %s\n%s; white = %2.1f", seat, cal_year, col_name, mean_2PP)
            
            # Overwrite old fill gradient and title:
            interp_raster <- interp_raster + scale_fill_gradientn(colours=colour_vec, space="rgb", values=colour_val) + ggtitle(title_str)
            
            output_map <- sprintf("images/%s/%s_%smoving_%d.png", sub_folder, seat, col_name, cal_year)
            png_print(output_map, 500, png_height, interp_raster)
            
            if ((seat %in% nonclassic_seats) & (in_seat)) {
              indices_tcp_year <- which(TCP_seatbooth.data$Year == cal_year)
              indices_tcp_seat <- which(TCP_seatbooth.data$Seat == seat)
              indices_tcp_seat_year <- intersect(indices_tcp_year, indices_tcp_seat)
              if (length(indices_tcp_seat_year) == 0) {
                # Plot the TPPmoving raster as the TCPmoving raster:
                title_str <- sprintf("%s, %s\nTCP: ALP v LNP; white = %2.1f", seat, cal_year, mean_2PP)
                tcp_interp_raster <- interp_raster + ggtitle(title_str)
                output_map <- sprintf("images/%s/%s_TCP_1moving_%d.png", sub_folder, seat, cal_year)
                png_print(output_map, 500, png_height, tcp_interp_raster)
              }
            }
          }
          
          # Make a moving colour scale for the TCP as well
          if (col_name == "TCP_1") {
            TCP_1_votes <- TCP_seatbooth.data.proj$TCP_1[indices_party_contested]
            formal_votes <- TCP_seatbooth.data.proj$Formal_votes[indices_party_contested]
            mean_TCP <- sum(TCP_1_votes * formal_votes) / sum(formal_votes)
            scale_low <- max(0, mean_TCP - 30)
            scale_high <- min(100, mean_TCP + 30)
            colour_val <- rescale(c(-1, 0, scale_low, mean_TCP, scale_high, 100, 101))
            title_str <- sprintf("%s, %s\nTCP: %s v %s; white = %2.1f", seat, cal_year, TCP_party1, TCP_party2, mean_TCP)
            # Overwrite old fill gradient and title:
            interp_raster <- interp_raster + scale_fill_gradientn(colours=colour_vec, space="rgb", values=colour_val) + ggtitle(title_str)
            output_map <- sprintf("images/%s/%s_%smoving_%d.png", sub_folder, seat, col_name, cal_year)
            png_print(output_map, 500, png_height, interp_raster)
          }
        }
      }
    }
  }
  
  # Delete the images that were only used to have a definitely-working colour bar:
  if (grepl(" ", seat)) {
    shell_str <- sprintf("del \"images\\%s\\*_colour.png\"", sub_folder)
  } else {
    shell_str <- sprintf("del images\\%s\\*_colour.png", sub_folder)
  }
  shell(shell_str)
  
  # Convert the frames into gifs.
  
  # A vector containing the names of various intermediate files created during the gif-making, mostly copies of
  # first and last frames.  They will be deleted at the end.
  intermediate_files <- character()
  
  # Number of copies to make of first and last frames:
  num_repeats <- 3
  gif_delay <- 75
  
  for (digit_num in digits) {
    mapfile_str <- sprintf("images/%s/map%d_%s.png", sub_folder, digit_num, seat)
    # Make new copies of the maps:
    for (ct in 1:num_repeats) {
      file_str_out <- sprintf("images/%s/map%d_%s_%d.png", sub_folder, digit_num, seat, ct)
      file.copy(mapfile_str, file_str_out, overwrite=T)
      intermediate_files <- c(intermediate_files, file_str_out)
    }
  }
 
  # Remove the last elemnt of vote_cols which is 999, replace the name manually with "TCP_1":
  col_names <- names(state_booth.data)[vote_cols[1:(length(vote_cols)-1)]]
  col_names <- c("ALP_2PPmoving", col_names, "TCP_1", "TCP_1moving")
  
  for (col_name in col_names) {
   got_image <- logical()
    for (cal_year in years) {
      file_str <- sprintf("images/%s/%s_%s_%d.png", sub_folder, seat, col_name, cal_year)
      got_image <- c(got_image, file.exists(file_str))
    }
    
    indices_year <- which(got_image)
    
    if (length(indices_year) > 0) {
      first_year <- years[indices_year[1]]
      last_year <- years[indices_year[length(indices_year)]]
      
      copy_frame_years <- unique(c(first_year, last_year))
      
      # Make copies of the first and last frames, so that the gif pauses:
      for (copy_year in copy_frame_years) {
        file_str_in <- sprintf("images/%s/%s_%s_%d.png", sub_folder, seat, col_name, copy_year)
        for (ct in 1:num_repeats) {
          file_str_out <- sprintf("images/%s/%s_%s_%d_%d.png", sub_folder, seat, col_name, copy_year, ct)
          file.copy(file_str_in, file_str_out, overwrite=T)
          intermediate_files <- c(intermediate_files, file_str_out)
        }
      }
      
      # The map will be the first frame, but we have to work out which map to use -- hopefully matches the digits
      # in the colourbar scale.
      main_img <- readImage(file_str_in)
      
      # Start by getting the right-hand edge of the plotting area of the main image.
      # Need to do two rows because one might be one of the white grid lines (precise row names assume 500x500px image)
      end_main <- raster_corner(main_img)
      real_end_r <- Re(end_main)
      
      num_digits <- digits[which(end_px == end_main)]
      
      if (length(num_digits) == 0) {
        print("Couldn't find the matching map.")
        asdflguh <- asdkfjh(asdfg)
      }
      
      if (grepl(" ", seat)) {
        partial_path <- sprintf("\"images\\%s\\%s_%s", sub_folder, seat, col_name)
        shell_str <- sprintf("convert -loop 0 -delay %d -layers Optimize \"images\\%s\\map%d*.png\" %s_*.png\" %s.gif\"", gif_delay, sub_folder, num_digits[1], partial_path, partial_path)
      } else {
        partial_path <- sprintf("images\\%s\\%s_%s", sub_folder, seat, col_name)
        shell_str <- sprintf("convert -loop 0 -delay %d -layers Optimize images\\%s\\map%d*.png %s_*.png %s.gif", gif_delay, sub_folder, num_digits[1], partial_path, partial_path)
      }
      shell(shell_str)
    }
  }
  file.remove(intermediate_files)
}
