Mapping Scotland’s wind farms in R and Google Maps

Posted on Updated on

I thought I would put up some R code I have for mapping wind farms in Scotland. This is hopefully of use to some readers.

Data source

Scottish National Heritage publishes a regularly updated – i.e. every 12 months – data set of all wind farms in Scotland by status. It can be downloaded here, and it requires you to go through a short registration programme.

The data is split into four wind farm categories: Installed, Approved, Application, and Scoping. SNH provides a shape file which shows the geographic region covered by each wind farm. However, the data set provided does not go into details about the capacity of wind farms, number of turbines etc. All we can do is map it and potentially work out the area of each named wind farm.

Here I will use ggmap to plot the data. An introduction to that package, which is based on ggplot2, is given here. ggmap let’s you plot Google Maps maps, and has the full functionality of ggplot2.

Reading in and sorting the data

First we need to read the data in and convert it to a more usable format. I’m an old fashioned scientist who struggles to understand coordinate systems other than longitude and latitude in degrees.

The following code will do all of that and read in the necessary packages.

## This file reads in and maps wind farms in Scotland -------------------------------------------------------
## Data is available from Scottish National Heritage --------------------------------------------------------
## http://www.snh.gov.uk/publications-data-and-research/snhi-information-service/naturalspaces/ -------------

## Read in required packages --------------------------------------------------------------------------------
require(maptools);library(rgdal);require(geosphere);require(ggthemes);require(ggmap)
options(stringsAsFactors = FALSE)

## Read in raw data -----------------------------------------------------------------------------------------
## Change the file name --------------------------------------------------------------------------------------
raw = readShapePoly('~/Dropbox/Scottish_Wind_Farms/WINDFARM_SCOTLAND.shp')

## proj4 transformations for converting to long/lat coordinate system -----
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"

## specify coordinate system as WGS84 ------------------------------------------------------
raw@proj4string   # slot will be empty
raw@proj4string = CRS(bng)

raw_wgs84 = spTransform(raw, CRS(wgs84))
coordinates(raw_wgs84)

The raw data is now in longitude/latitude format. I need a data set that is more usable in ggmaps. So I will create a fortified data frame, which lets me identify the status of each farm.

This code does that:

## Fortify the data ----------------
tracker = 1
for(ii in unique(raw_wgs84$STATUS))
  {
  print(ii)
  op_wgs84 = subset( raw_wgs84, STATUS == ii )
  if( tracker == 1)
    all_farms = data.frame(fortify(op_wgs84), Status = ii ) else
      all_farms = rbind(all_farms,data.frame(fortify(op_wgs84), Status = ii))
    
  tracker = tracker + 1
}

Mapping the wind farms

We now have a data frame that can be mapped using ggmap. First, I will show all wind farms on a single map. To do this, I will use the colour blind scale from the package “ggthemes”.

map = get_map( zoom = 6, location = "Scotland",maptype = "hybrid")

ggmap(map)+
  geom_polygon(data = all_farms, aes(long, lat, group = group, fill = Status))+
  scale_fill_colorblind()+
  xlim(c(-8,0))+
  ylim(c(55,58.5))+
  theme(legend.position = c(0.85, 0.85))+
  xlab("Longitude")+
  ylab("Latitude")

The result is this:

Single_Map
This is perhaps unsatisfying. Even if the palette is for colour blind people, I find it hard to distinguish what is going on. It is probably better to use facetting to split the map out.

facet_wrap can be used to give four separate maps for the four different categories.

ggmap(map)+
  geom_polygon(data = all_farms, aes(long, lat, group = group), fill = "red")+
  scale_fill_colorblind()+
  xlim(c(-8,0))+
  ylim(c(55,58.5))+
  theme(legend.position = c(0.85, 0.85))+
  xlab("Longitude")+
  ylab("Latitude")+
  facet_wrap(~Status)

And we get this:

Panel_Map
Alternatively we can colour the wind farms yellow, which might aid some colour blind people:

Panel_MapAnd that’s that. I might write a follow up post looking at the areas of wind farms and their power density.The result is this: