{"id":81,"date":"2017-04-19T11:59:05","date_gmt":"2017-04-19T10:59:05","guid":{"rendered":"http:\/\/research.shca.ed.ac.uk\/past-by-numbers\/?p=81"},"modified":"2017-04-28T10:12:58","modified_gmt":"2017-04-28T09:12:58","slug":"loading-and-visualizing-shapefiles-in-r","status":"publish","type":"post","link":"http:\/\/research.shca.ed.ac.uk\/past-by-numbers\/2017\/04\/19\/loading-and-visualizing-shapefiles-in-r\/","title":{"rendered":"Loading and visualizing Shapefiles in R"},"content":{"rendered":"

Geographical Information Systems such as QGIS<\/a> are common in the typical archaeologist’s toolbox. However, the complexities of our datasets mean that we need additional tools to identify patterns in time and space.<\/p>\n

Imagine that you are working with a set of points in a vector format such as a Shapefile. It can be a set of settlements, C14 dates or shipwrecks…it doesn’t matter; you will want to combine spatial analysis with other methods. This is particularly true when you are trying to get a feeling of the dataset by performing Exploratory Data Analysis<\/em>, so you need to compile summary statistics and create some meaningful visualisations.<\/p>\n

My personal approach is to combine QGIS with R and exploit the awesome ggplot<\/a> package to get some insight into the dataset beyond spatial patterns. R has some nice spatial functionality, but the cool thing is to integrate these methods with other plots to get a general picture of the case study you are studying.<\/p>\n

Example: Aircraft crashes in the Orkney islands<\/h2>\n

The Orkney islands were a key location for the British war effort during the Two World Wars because it hosted Scapa Flow<\/a>: the main naval base of the Royal Navy. For this reason it saw an intense aerial activity of all sorts: squadrons defending the base, German bombers attacking it, aircraft carrier operations…you name it. Inevitably this activity generated aircraft crashes due to accidents and combat and these events have been compiled by different initiatives such as the Project Adair-Whitaker<\/a> or the ARGOS group<\/a>.<\/p>\n

I created a nice Shapefile of aircraft crash sites based on the information provided by Canmore<\/a>. How can I explore the dataset beyond the spatial dynamics? How many German aircrafts were shot down? What types and squadrons sustained more losses? Let’s take a look!<\/p>\n

Download Data<\/h3>\n

The first thing to do is to load the Shapefile. This format requires of several files so I created a zip<\/em> file with all the set. First of all we need to download it to a temporary directory and extract its contents:<\/p>\n[generic linenumbers=”False”]\ntmp <- tempdir()
\nurl <- "https:\/\/github.com\/xrubio\/pastByNumbers\/raw\/master\/data\/aircraft_orkney.zip"
\nfile <- basename(url)
\ndownload.file(url, file)
\nunzip(file, exdir = tmp)
\n[\/generic]\n

Load the spatial data<\/h3>\n

We will use the rgdal<\/a> library to load the contents of the Shapefile into R:
\n[generic linenumbers=”False”]\nlibrary(rgdal)
\naircraftShp <- readOGR(dsn = tmp, layer="aircraft_orkney")
\nstr(aircraftShp)
\n[\/generic]\n

aircraftShp<\/em> is a variable of type SpatialPointsDataFrame<\/em> from package sp<\/em>. It is a container of some interesting information such as the bounding box delimiting our spatial entities or the Coordinate Reference System that it is being used. Ideally we would like to get the attributes of the spatial entities as a data frame so we can use it with most R packages:<\/p>\n[generic linenumbers=”False”]\naircraft <- as.data.frame(aircraftShp)
\nstr(aircraft)
\n[\/generic]\n

Maps and planes<\/h3>\n

First of all let’s plot the spatial coordinates of the aircraft crash sites against a nice background. I created a rectangle location based on the bounding box of aircraftShp<\/em>:<\/p>\n[generic linenumbers=”False”]\naircraftShp@bbox
\nlocation <- c(-3.85, 58.7, -1.95, 59.4)
\n[\/generic]\n

We can use ggmap<\/a> to download the background:
\n[generic linenumbers=”False”]\nlibrary(ggplot2)
\nlibrary(ggmap)
\nbgMap <- get_map(location=location, source="stamen", maptype="watercolor")
\nggmap(bgMap) + geom_point(data=aircraft, aes(x=coords.x1, y=coords.x2))
\n[\/generic]\n

\"\"
Aircraft crash sites around Scapa Flow<\/figcaption><\/figure>\n

Exploring other dimensions<\/h3>\n

This is ok but it is something we can already do with any GIS so…what’s the deal? The thing is that aircraft<\/em> is an R Data Frame so we can plot any variable using a graphic library, something that R does way better than any GIS. For example, let’s take a look at the temporal dimension:<\/p>\n[generic linenumbers=”False”]\nggplot(aircraft, aes(x=year, fill=force)) + geom_histogram(binwidth=1) + facet_wrap(~war, scales=”free_x”)
\n[\/generic]\n

\"\"
Losses by force and year<\/figcaption><\/figure>\n

We can also create a table of planes used by the different British squadrons that suffered losses:<\/p>\n[generic linenumbers=”False”]\nggplot(na.omit(aircraft), aes(y=type, col=force, x=as.factor(sqdn))) + geom_raster()
\n[\/generic]\n

\"\"
Aircraft types per squadron<\/figcaption><\/figure>\n

Combining everything<\/h3>\n

We can finally create an infographic combining space, time and other stuff with the library gridExtra<\/a>:<\/p>\n

First, we create nice version of the plots we just saw:<\/p>\n[generic linenumbers=”False”]\ng1 <- ggmap(bgMap, extent="panel", darken = c(.4,"white")) + geom_point(data=aircraft, aes(x=coords.x1, y=coords.x2, col=war), size=2) + geom_text(data=aircraft, aes(x=coords.x1, y=coords.x2,label=type), family = "Trebuchet MS", color="grey40", check_overlap=T, nudge_y=0.01) + theme(legend.position="bottom") + ggtitle("aircraft losses") + xlab("") + ylab("") + scale_color_manual(values=c("indianred2", "goldenrod1"))<\/p>\n

g2 <- ggplot(na.omit(aircraft), aes(y=type, col=force, x=as.factor(sqdn))) + geom_point(size=3) + theme_bw() + theme(legend.position="bottom") + xlab("squadron") + ylab("") + ggtitle("Britain – aircraft types per squadron")<\/p>\n

g3 <- ggplot(aircraft, aes(x=year, fill=force)) + geom_histogram(binwidth=1, col="black") + facet_wrap(~war, scales="free_x", ncol=1) + theme_bw() + theme(legend.position="bottom") + ggtitle("losses per year") + xlab("") + ylab("")
\n[\/generic]\n

And we can then combine the 3 plots into a nice visualization:
\n[generic linenumbers=”False”]\nlibrary(gridExtra)
\ngrid.arrange(arrangeGrob(arrangeGrob(g1,g2, heights=c(2\/3,1\/3), ncol=1), g3, widths=c(2\/3, 1\/3), ncol=2))
\n[\/generic]\n

<\/a>
Aicraft Losses in Orkney islands during the two World Wars<\/figcaption><\/figure>\n

Interpretation<\/h3>\n

These visualizations allow us to identify several interesting things going on:<\/p>\n