Mass. Power Outages
One of the power companies in Mass. was kind enough to use a mapping service that makes the power outage data easily available. We’ll need some pacakes for this.
library(maptools) # geo stuff
library(rgeos) # moar geo stuff
library(rgdal) # even moar geo stuff
library(ggplot2) # best plotting evar
library(ggthemes) # good themes
library(ggalt) # better projections
library(viridis) # better colors
library(scales) # formattingThis is an “empty” shapefile that contains the data about the power service status in each town.
outage_url <- "http://mema.mapsonline.net/power_outage_public.geojson"
ma_outages <- suppressWarnings(readOGR(outage_url, "OGRGeoJSON",
verbose=FALSE, stringsAsFactors=FALSE))You can see what it hold by looking at the data slot:
dplyr::glimpse(ma_outages@data)
## Observations: 32
## Variables: 10
## $ town (chr) "BERLIN", "BOSTON", "BOURNE", "BREWSTER", "...
## $ county (chr) "WORCESTER", "SUFFOLK", "BARNSTABLE", "BARN...
## $ memaregion (chr) "4", "1", "2", "2", "2", "2", "1", "4", "2"...
## $ utility (chr) "National Grid", "Eversource - Eastern MA",...
## $ total_cust (int) 1348, 301361, 11312, 8845, 5595, 8355, 1551...
## $ total_cust_utility (int) 1348, 301361, 11312, 8845, 5595, 8355, 1551...
## $ no_power (int) 7, 1, 63, 94, 1, 65, 15, 1, 1, 167, 30, 3, ...
## $ pct_nopow (dbl) 5.192878e-03, 3.318279e-06, 5.569307e-03, 1...
## $ notes (chr) "", "", "", "", "", "", "", "", "", "", "",...
## $ last_update (chr) "2016/01/23 05:00 pm", "2016/01/23 05:00 pm...Now we need some county and town maps. Mass. gov provides those (and they are shockingly decent). We conver them to lat/lon so it speeds up coord_proj and/or coord_map. And, we download them in such a way as to cache them and not re-download them every time.
cnty_url <- "http://wsgw.mass.gov/data/gispub/shape/state/counties.zip"
cnty_fil <- basename(cnty_url)
if (!file.exists(cnty_fil)) download.file(cnty_url, cnty_fil)
unzip(cnty_fil)
counties <- readOGR("COUNTIES_POLY.shp", "COUNTIES_POLY",
verbose=FALSE, stringsAsFactors=FALSE)
counties <- SpatialPolygonsDataFrame(spTransform(counties, CRS("+proj=longlat")),
counties@data)
twn_url <- "http://wsgw.mass.gov/data/gispub/shape/state/towns.zip"
twn_fil <- basename(twn_url)
if (!file.exists(twn_fil)) download.file(twn_url, twn_fil)
unzip(twn_fil)
towns <- readOGR("TOWNS_POLY.shp", "TOWNS_POLY",
verbose=FALSE, stringsAsFactors=FALSE)
towns <- SpatialPolygonsDataFrame(spTransform(towns, CRS("+proj=longlat")),
towns@data)We need those polygons in a format we can use with ggplot2, so we “fortify” them (turn the complex SpatialPolygonsDataFrames into a plain old data.frame) using the county and town names as a way to index the polygons.
counties_map <- fortify(counties, region="COUNTY")
towns_map <- fortify(towns, region="TOWN")This is the projection most folks are used to when they see a Mass. map:
ma_proj <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"Now, we plot the map. First we lay down the towns, then fill in the towns with outages, then lay down the county borders.
gg <- ggplot()
gg <- gg + geom_map(map=towns_map, data=towns_map,
aes(x=long, y=lat, map_id=id),
color="#7f7f7f", fill=NA, size=0.15)
gg <- gg + geom_map(map=towns_map, data=ma_outages@data,
aes(fill=pct_nopow, map_id=town),
color="white", size=0.15)
gg <- gg + geom_map(map=counties_map, data=counties_map,
aes(x=long, y=lat, map_id=id),
color="#2b2b2b", fill=NA, size=0.5)
gg <- gg + scale_fill_viridis(name="% w/o\nPower", labels=percent)
gg <- gg + coord_proj(ma_proj)
gg <- gg + labs(title="Mass. Power Outages",
x=sprintf("%s total w/o power",
comma(sum(ma_outages@data$no_power))))
gg <- gg + theme_map()
gg <- gg + theme(legend.position="right")
gg <- gg + theme(axis.title=element_text())
gg <- gg + theme(axis.title.y=element_blank())
gg <- gg + theme(plot.title=element_text(hjust=0, size=16, face="bold"))
gg
Every time you re-run this you’ll get new outage data. You could even write a script that grabs the outage data regularly and tracks it over time.