I found an interesting London Fire Brigade data-set which I wanted to see visualised. Here it is for all of London:

London FireStation call-outs

And here for only two areas (i.e. all the fire stations serving these two areas/OAs – personal to me):

London FireStation call-outs (select)

The exercise of plotting just the fire-stations serving my regions of interest was basically:

  1. Generate the polygons which define my areas of interest
  2. Intersect the fire data (incidence) with those polygons and keep just those fires
  3. Retrieve all the fire-stations which served those fires
  4. Finally, keep just those fire-stations

Such a task is trivial when working with spatial-data. I mostly made use of the gIntersection command in the rgeos package.

This project could have been potentially more interesting if I decided to generate a route linking the origin to the destination (rather than a straight-line). Such a task could be accomplished using the Google Directions API, GraphHopper API or ArcGIS. However, I decided I wanted to see what a minimalist map would look like.

R-script:

library(ggplot2)
library(rgdal)
library(foreach)
library(plyr)
####################################################################################################
setwd(".../8. FireMap")

###################
# [A] Import Fire Data
###################
input_csv <- list()
counter = 0
foreach (part_file=c("LFB incident data 1 Jan 2012 to 30 Sep 2015.csv", 
                     "LFB incident data 1 Jan 2009 to 31 Dec 2011.csv")) %do% {
                       
  print(part_file)
  counter = counter + 1
  input_csv[[counter]] <- read.csv(part_file, header = TRUE, as.is = TRUE)
  
  # Reshape from wide to long (so no difference if first pump or second pump)
  input_csv_pump1 <- input_csv[[counter]][,c("Easting_rounded", "Northing_rounded", "FirstPumpArriving_DeployedFromStation", "DateOfCall")]
  names(input_csv_pump1) <- c("Easting", "Northing", "Station", "Date")
  input_csv_pump2 <- input_csv[[counter]][,c("Easting_rounded", "Northing_rounded", "SecondPumpArriving_DeployedFromStation", "DateOfCall")]
  names(input_csv_pump2) <- c("Easting", "Northing", "Station", "Date")
  
  input_csv[[counter]] <- rbind(input_csv_pump1,input_csv_pump2)
  rm(input_csv_pump1)
  rm(input_csv_pump2)
  
  # Remove if no station information
  input_csv[[counter]] <- input_csv[[counter]][input_csv[[counter]]$Station!="",]
}

# Combine the two data-sets
final.data <- do.call(rbind, input_csv)
rm(input_csv)

# Cut to one year (30 Sep 2014 to 30 Sep 2015)
# I guess I could have not imported the second data-set then
final.data$rdate <- as.Date(final.data$Date, "%d-%b-%y")
final.data <- final.data[final.data$rdate >= as.Date("2014/09/30"),]
final.data <- final.data[,c("Easting", "Northing", "Station")]

# Collapse count
final.data$c <- 1
final.data  <- ddply(final.data , c("Easting", "Northing", "Station"), summarize, count=sum(c)) 
# Sort in descending order
final.data <- final.data[order(-final.data$count),]
final.data <- final.data[final.data$Station!="Esher",]

####################
# [B] Import Fire Station Locations
###################
# Produced from: http://www.london-fire.gov.uk/A-ZFireStations.asp
# Added 10  fire stations closed in late 2014 manually
# Not sure why 'Esher' Firestation included in data (maybe only as secondary)

stations_geo <- read.csv("stations.csv")
# Merge on the station latitude/longitude
merged_geo <- merge(final.data, stations_geo, by="Station", all.x=TRUE)
names(merged_geo) <- c("Station", "inc_Easting", "inc_Northing", "Count", "Station_Postcode", "Station_Lat", "Station_Lon")
assertthat::noNA(merged_geo)

# Transform the Easting and Northing into Latitude/Longitudes
# Key proj4 transformations
wgs84 = '+proj=longlat +datum=WGS84'
bng = '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs'

# Create spatial object where eastings is x-axis and northings is y-axis
coordinates(merged_geo) <- c("inc_Easting", "inc_Northing")
# Set projection of data
proj4string(merged_geo) <- CRS(bng)
# BNG to WGS84 conversion
merged_geo_wgs84 = spTransform(merged_geo, CRS(wgs84))
# Data-frame
converted_coordinates <- as.data.frame(merged_geo_wgs84)

# Segments to plot
segments_df <- converted_coordinates[,c("Station_Lat", "Station_Lon", "inc_Easting", "inc_Northing", "Count", "Station")]
names(segments_df) <- c("from_lat", "from_lon", "to_lon", "to_lat", "count", "station")
segments_df$id <- row.names(segments_df)

##############################################
# [C] Keep only fires in regions_of_interest 
##############################################

# Regions_of_Interest
london_shape_loc <- 'wards'
london <- readOGR(london_shape_loc, layer="London_Ward_CityMerged") 
london <- spTransform(london, CRS("+proj=longlat +datum=WGS84"))
regions_of_interest <- london[london$NAME %in% c("Coldharbour","Greenford Broadway"),]
plot(regions_of_interest)

library(rgeos)
# Intersect the fire incidents with the regions of interest to drop those outside
regions_of_interest_wg84 <- gIntersection(regions_of_interest,merged_geo_wgs84, byid=TRUE)

# Convert SpatialPoints back into a data-frame
regions_of_interest_coordinates <- as.data.frame(regions_of_interest_wg84)
regions_of_interest_coordinates$id <- row.names(regions_of_interest_coordinates)
regions_of_interest_coordinates$id <- gsub( ".* ","",regions_of_interest_coordinates$id)

#Trim to those stations which serve regions_of_interest lane
segment_regions_of_interest <- merge(regions_of_interest_coordinates,
                             segments_df, 
                             by="id", all.x=TRUE)
stations <- as.data.frame(unique(segment_regions_of_interest$station))
names(stations) <- c("station")
# All of these stations served fires in our areas of interest
print(stations)
# Keep just those stations
segments_df_trim <- merge(stations,segments_df,
                           by='station',
                           all.x=TRUE)
assertthat::noNA(segments_df_trim)

###########
# [D] Plot
###########
#Remove axis
xquiet<- scale_x_continuous("", breaks=NULL)
yquiet<-scale_y_continuous("", breaks=NULL)
quiet<-list(xquiet, yquiet)

# Base-map (Greater London)
london_shape_loc <- 'london_sport'
london <- readOGR(london_shape_loc, layer="london_sport") 
london <- spTransform(london, CRS("+proj=longlat +datum=WGS84"))

# ADD Water (clipped from OS)
rivers_loc <- 'os_vector_map_merged(IK)'
rivers <- readOGR(rivers_loc, layer="TQ_TidalWater_Merge_Clip") 
rivers <- spTransform(rivers, CRS("+proj=longlat +datum=WGS84"))

# ALL
ggplot() +
  # River as background
  geom_polygon(data=rivers, aes(long, lat, group=group),
               color='white', fill='white', size=1, alpha=1)+
  # Plot segments
  geom_segment(data=segments_df, aes(x=from_lon, xend = to_lon, y=from_lat, yend=to_lat, alpha=count, 
                                          color=station), size = 1) +
  scale_alpha_continuous(range = c(0.05, 0.85)) +
  # Base
  geom_polygon(data=london, aes(long, lat, group=group),
               color='white', fill=NA, size=0.8, alpha=0.8)+
  theme(panel.background = element_rect(fill='#2C3539',colour='#2C3539')) +
  #http://geepeeex.com/LongitudesAndLatitudes.htm
  quiet + coord_equal(ratio=122/78)
  
ggsave(file='firestations_all_V2.png', width=40,height=20,dpi=360)