Spatial Modelling 1
## Shapefile type: Polygon, (5), # of Shapes: 10## Shapefile type: Polygon, (5), # of Shapes: 155In this session we explore the details of spatial modelling, mapping and more advanced GIS work as they support social and cultural analytics. We agree with https://pakillo.github.io/R-GIS-tutorial/ that R is also great for spatial analysis. Check out http://spatialanalysis.co.uk/wp-content/uploads/2012/02/bike_ggplot.png for an example visualisation.
Before we proceed, you need to register for a Google Maps API key. First, you need to go to https://developers.google.com/maps/documentation/maps-static/get-api-key and follow the instructions under Quick Guide. Once you have a key, follow the instructions at https://github.com/dkahle/ggmap to register the Google Map key in R. If you have not done this yet, please leave swirl at the next question and return once you have a registered Google API key.
We use a number of libraries, starting with maptools, which provides us with a functions that help with reading and processing spatial data. You can load the package with library(maptools). Please note that this also loads the standard sp package, which provides classes for storing and manipulating different types of spatial data. It is kind of a standard in R. So new spatial packages expect data to be an sp object.
library(maptools)In your environment you should see spatial datasets related to Haiti. In particular, we will investigate a dataset related to the earthquake that took place there in 2010. The dataset was crowd-sourced using the Ushahidi platform (https://wiki.ushahidi.com/display/WIKI/Datasets), a non-profit software company that enables crowdsourcing of information related to natural disasters and geopolitical events via web form, text messages or social media. We have some data collected during the 2010 Haiti earthquake crisis and aftermath, which we borrowed from Wes Mckinney’s Python for Data Analysis. Let us take a look at the data with head(haiti_quake).
head(haiti_quake)## # A tibble: 6 x 55
## Serial INCIDENT_TITLE INCIDENT_DATE LOCATION DESCRIPTION CATEGORY
## <int> <chr> <dttm> <chr> <chr> <chr>
## 1 4052 * URGENT * Ty… 2010-07-05 17:26:00 Jacmel,… Birthing C… 1. Urge…
## 2 4042 Citi Soleil s… 2010-05-18 16:26:00 Citi So… "We are wo… 1. Urge…
## 3 4041 Radio Commerc… 2010-04-26 13:14:00 Radio C… "i'm Louin… 5e. Com…
## 4 4040 Contaminated … 2010-04-26 14:19:00 Marc ne… "How do we… 4. Mena…
## 5 4039 Violence at &… 2010-04-26 14:27:00 unable … "Goodnight… 4. Mena…
## 6 4038 No electricit… 2010-03-15 10:58:00 Pernier "why the p… 2. Urge…
## # ... with 49 more variables: LATITUDE <dbl>, LONGITUDE <dbl>,
## # APPROVED <chr>, VERIFIED <chr>, category_1 <int>, category_1a <int>,
## # category_1b <int>, category_1c <int>, category_1d <int>,
## # category_2 <int>, category_2a <int>, category_2b <int>,
## # category_2c <int>, category_2d <int>, category_2e <int>,
## # category_2f <int>, category_3 <int>, category_3a <int>,
## # category_3b <int>, category_3c <int>, category_3d <int>,
## # category_3e <int>, category_4 <int>, category_4a <int>,
## # category_4c <int>, category_4e <int>, category_5 <int>,
## # category_5a <int>, category_5b <int>, category_5c <int>,
## # category_5d <int>, category_5e <int>, category_6 <int>,
## # category_6a <int>, category_6b <int>, category_6c <int>,
## # category_7 <int>, category_7a <int>, category_7b <int>,
## # category_7c <int>, category_7d <int>, category_7g <int>,
## # category_7h <int>, category_8 <int>, category_8a <int>,
## # category_8c <int>, category_8d <int>, category_8e <int>,
## # category_8f <int>Each observation is a report by a mobile phone indicating an emergency or some other problem after the quake. Each has an associated timestamp and a location as latitude and longitude, which we will use to map the incidents. The many category fields might be confusing, but they relate to specific events such as occurrence of violence or food shortage.
haiti_quake is a normal data frame. But if you check your environment, you can see that it also contains other spatial data such as SpatialPointsDateFrames, SpatialPolygonsDataFrames and SpatialLinesDataFrames. The names should be quite telling but we will discuss them in detail, too. There are two spatial datasets for Haiti; one containing the departmental boundaries of Haiti and another one with all the roads.
As a first step, we transform haiti_roads into a data frame and take a look at it in a format we understand. Run head(data.frame(haiti_roads)) to print the first rows of the roads data. You can see that haiti_roads contains road names, but also additional information such as whether these are one-way, etc.
head(data.frame(haiti_roads))## name type osm_id ref oneway bridge
## 0 Rue Magny residential 23860949 <NA> 0 0
## 1 Boulevard Henri Truman primary 23861001 <NA> 1 0
## 2 Boulevard Henri Truman primary 23861018 <NA> 1 0
## 3 Boulevard Henri Truman primary 23861022 <NA> 1 1
## 4 Nationale No 2 primary 23861052 2 1 0
## 5 Boulevard des Industries secondary 23861200 <NA> 1 0Let us also take a look at the departments of Haiti with head(data.frame(haiti_departments)). Again, we can see spatial information but also additional information per department such as the population density.
head(data.frame(haiti_departments))## ID_ADM1 ADM1 Shape_Leng Shape_Area Capital Area Population
## 1 1 Ouest 7.216540 0.4267732 Port-au-Prince 4827 4029705
## 2 2 Sud-Est 4.042320 0.1740848 Jacmel 2023 632601
## 3 3 Nord 3.834594 0.1810655 Cap-Haïtien 2106 1067177
## 4 4 Nord-Est 2.478687 0.1400675 Fort-Liberté 1805 393967
## 5 5 Artibonite 4.993639 0.4176159 Gonaives 4984 1727524
## 6 6 Centre 3.601807 0.2969891 Hinche 3675 746236
## Population_density
## 1 834.8260
## 2 312.7044
## 3 506.7317
## 4 218.2643
## 5 346.6140
## 6 203.0574To proceed with the spatial analysis of the Ushahidi data, we would like to create a spatial dataset of points with additional information. This SpatialPointsDataFrame will contain the information in haiti_quake and is quite easy to produce in R. We first need to extract the longitude and latitude information from haiti_quake. Hopefully, you remember cbind and can run without problems coords_tmp <- cbind(haiti_quake\(LONGITUDE, haiti_quake\)LATITUDE).
coords_tmp <- cbind(haiti_quake$LONGITUDE, haiti_quake$LATITUDE)Next we simply create the SpatialPointsDataFrame with these coordinates and the other information from haiti_quake. Please type haiti_quake_sdf <- SpatialPointsDataFrame(coords_tmp, data = data.frame(haiti_quake), proj4string=crswgs84). In case you wonder about proj4string, you should read through Part III of https://cran.r-project.org/doc/contrib/intro-spatial-rl.pdf. It explains that CRS stands for Coordinate Reference System of spatial objects. It sets where coordinates are placed on the Earth’s surface, which is difficult because we need to project from an almot spherical earth surface to a two-dimensional map. CRS information should always be provided, because the projections from the earth sphere to the map can be very different. Every spatial object should have a coordinate reference system (CRS) associated with it. In this case, we use the ‘British National Grid - WGS84’, a commonly used CRS.
haiti_quake_sdf <- SpatialPointsDataFrame(coords_tmp, data = data.frame(haiti_quake), proj4string=crswgs84)I have created a little demo for you to demonstrate how projections can influence the way the world looks like on a map. Run demo_projections() in order to compare four common projections for mapping the world onto a two-dimensional surface. You can clearly see how different continents look. You might want to check online the details of the four projections.
demo_projections()## OGR: Corrupt data## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3## OGR: Corrupt data## OGR: Corrupt data
## OGR: Corrupt data
## OGR: Corrupt data
## OGR: Corrupt data
## OGR: Corrupt data
## OGR: Corrupt dataBack to Haiti. Check the new spatial dataset haiti_quake_sdf with head(data.frame(haiti_quake_sdf)).
head(data.frame(haiti_quake_sdf))## Serial INCIDENT_TITLE
## 1 4052 * URGENT * Type O blood donations needed in #Jacmel #Haiti
## 2 4042 Citi Soleil school
## 3 4041 Radio Commerce in Sarthe
## 4 4040 Contaminated water in Baraderes.
## 5 4039 Violence at "arcahaie bas Saint-Ard"
## 6 4038 No electricity in pernier
## INCIDENT_DATE
## 1 2010-07-05 17:26:00
## 2 2010-05-18 16:26:00
## 3 2010-04-26 13:14:00
## 4 2010-04-26 14:19:00
## 5 2010-04-26 14:27:00
## 6 2010-03-15 10:58:00
## LOCATION
## 1 Jacmel, Haiti
## 2 Citi Soleil, Haiti
## 3 Radio Commerce Shelter, Sarthe
## 4 Marc near Baraderes
## 5 unable to find "arcahaie bas Saint-Ard" on map
## 6 Pernier
## DESCRIPTION
## 1 Birthing Clinic in Jacmel #Haiti urgently needs a O- blood transfusion 4 woman who just gave birth. please see Via @coreyrateau (Twitter)
## 2 We are working with Haitian (NGO) -The Christian Foundation of Succor to Haiti. \r \rPasteur Dieusel operates a School in Citi Soleil. We have secured one permanent school building and two large tents to temporily replace/relocate the damaged/dangerous buildings that are along the putrid river. \r \rHe needs funding for the construction of 5 more permanent buildings in the new location to provide space for 600 more children.\r \rPasteur Dieusel is a good man. He would appreciate any assistance you can provide with his efforts to help the Citi Soleil children.\r \rCharles G. Willard\rcgw@wrtb.org
## 3 i'm Louinel from Sarthe. I'd to know what can the government do for us in Radio Commerce Shelter\r\rIDUshahidi: \t13715162
## 4 How do we treat water in areas without Pipe?\t\rIDUshahidi: \t13688395
## 5 Goodnight at (arcahaie bas Saint-Ard) 2 young boys tiga and pastor they are brothers they put with big doudou they give avinel 3 bullets Tuesday April 20 .They father's named is tiyonyon, hurry come get these thugs because they have big guns in the hands. \rIDUshahidi: \t13710901
## 6 why the people who lives in pernier doesn't find in Electricity ( EDH) \rIDUshahidi: \t12176330
## CATEGORY
## 1 1. Urgences | Emergency, 3. Public Health,
## 2 1. Urgences | Emergency,
## 3 5e. Communication lines down,
## 4 4. Menaces | Security Threats, 4e. Assainissement eau et hygiene | Water sanitation and hygiene promotion,
## 5 4. Menaces | Security Threats,
## 6 2. Urgences logistiques | Vital Lines, 2f. Sans courant | Power Outage,
## LATITUDE LONGITUDE APPROVED VERIFIED category_1 category_1a category_1b
## 1 18.23333 -72.53333 YES NO 1 0 0
## 2 18.57108 -72.33467 YES NO 1 0 0
## 3 18.59371 -72.31008 YES NO 0 0 0
## 4 18.48280 -73.63880 YES NO 0 0 0
## 5 18.41500 -73.19500 YES NO 0 0 0
## 6 18.51744 -72.23684 YES NO 0 0 0
## category_1c category_1d category_2 category_2a category_2b category_2c
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 1 0 0 0
## category_2d category_2e category_2f category_3 category_3a category_3b
## 1 0 0 0 1 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 1 0 0 0
## category_3c category_3d category_3e category_4 category_4a category_4c
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 1 0 0
## 5 0 0 0 1 0 0
## 6 0 0 0 0 0 0
## category_4e category_5 category_5a category_5b category_5c category_5d
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 1 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0
## category_5e category_6 category_6a category_6b category_6c category_7
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 1 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0
## category_7a category_7b category_7c category_7d category_7g category_7h
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0
## category_8 category_8a category_8c category_8d category_8e category_8f
## 1 0 0 0 0 0 0
## 2 0 0 0 0 0 0
## 3 0 0 0 0 0 0
## 4 0 0 0 0 0 0
## 5 0 0 0 0 0 0
## 6 0 0 0 0 0 0
## coords.x1 coords.x2 optional
## 1 -72.53333 18.23333 TRUE
## 2 -72.33467 18.57108 TRUE
## 3 -72.31008 18.59371 TRUE
## 4 -73.63880 18.48280 TRUE
## 5 -73.19500 18.41500 TRUE
## 6 -72.23684 18.51744 TRUEThis way you can create any spatial data frame in R. Of course, creating lines and polygons requires more complicated definitions than simple points. There are many advantages of using the spatial data frames. The first one is that it is very easy to plot them by simply typing in plot(haiti_departments).
plot(haiti_departments)
You should see the departments of Haiti. We can add the roads by adding another plot command. Type in plot(haiti_roads, col = “lightgrey”, add = TRUE). Do you remember the add attribute?
plot(haiti_departments)
plot(haiti_roads, col = 'lightgrey', add = TRUE)
It will look a bit crowded but we can add the Ushahidi messages’ locations to the map, too. Type in plot(haiti_quake_sdf, col = “red”, pch = 1, add = TRUE).
plot(haiti_departments)
plot(haiti_roads, col = 'lightgrey', add = TRUE)
plot(haiti_quake_sdf, col = 'red', pch = 1, add = TRUE)
I hope you can see that some Haiti departments have seen more emergency notifications than others. We suspect there is some kind of grouping of emergencies going on here. We will come back later to how to formally investigate this. Let us finally add a title with title(“Ushahidi Messages Haiti Earthquake”), which completes our plot.
plot(haiti_departments)
plot(haiti_roads, col = 'lightgrey', add = TRUE)
plot(haiti_quake_sdf, col = 'red', pch = 1, add = TRUE)
title("Ushahidi Messages Haiti Earthquake")
We can subset spatial data like any other data we have met. For instance, we can identify the northern departments in Haiti. Start a new plot with plot(haiti_departments).
plot(haiti_departments)
Entering plot(subset(haiti_departments, ADM1 %in% c(‘Nord’, ‘Nord-Est’)), col = “lightblue”, add = TRUE), will colour the 2 northern departments in blue. Please, also note the %in% operator that checks that ADM1 is in the vector of northern departments.
plot(haiti_departments)
plot(subset(haiti_departments, ADM1 %in% c("Nord", "Nord-Est")), col = "lightblue", add = TRUE)
In the haiti_quake_sdf dataset, the many category columns describe the different kinds of emergencies. Let us compare 2 types (category_1 and category_2) to see whether we can identify spatial pattern. Let us start a new plot with plot(haiti_departments).
plot(haiti_departments)
Now, we only add category 1 emergencies. These are those observations which have a value larger than 0 in the category_1 column. Run plot(subset(haiti_quake_sdf, category_1 > 0), col = “red”, pch = 1, alpha = 0.3, add = TRUE).
plot(haiti_departments)
plot(subset(haiti_quake_sdf, category_1 > 0), col = "red", pch = 1, add = TRUE)
Category 1 emergencies are obviously clustered in Port au Prince, the capital of Haiti. In order to add category_2, execute plot(subset(haiti_quake_sdf, category_2 > 0), col = “blue”, pch = 1, add = TRUE).
plot(haiti_departments)
plot(subset(haiti_quake_sdf, category_1 > 0), col = "red", pch = 1, add = TRUE)
plot(subset(haiti_quake_sdf, category_2 > 0), col = "blue", pch = 1, add = TRUE)
Port-au-Prince is again a focus point of these emergencies, but there are also a few other distinct report locations compared to category_2. Let us see next whether we can identify a relationship between the time of the incident and its location. To this end, we will use INCIDENT_DATE, which a good person has prepared for us as a timestamp feature. We start a new plot with plot(haiti_departments).
plot(haiti_departments)
We know this outline of Haiti by now quite well. We would like to group five different groups of incidents’ times and colour them accordingly on the map. With palette(c(“blue”,“red”,“green”,“black”,“yellow”)), we set the five colour palette for a col argument with a numeric index. Try it.
palette(c("blue","red","green","black","yellow"))The plot is now easy. Just run plot(haiti_quake_sdf, col = cut(haiti_quake_sdf\(INCIDENT_DATE, breaks = quantile(haiti_quake_sdf\)INCIDENT_DATE)), pch = 1, add = TRUE). The quantile function is used to split up the dataset according to 5 different time slots. The cut function then cuts this quantile range into the 5 different time slots.
plot(haiti_departments)
plot(haiti_quake_sdf, col = cut(haiti_quake_sdf$INCIDENT_DATE, breaks = quantile(haiti_quake_sdf$INCIDENT_DATE)), pch = 1, add = TRUE)
The resulting plot does not really reveal distinct times for particular location reports. There is are some incidents further away from the capital of Haiti that obviously relate to particular events but, otherwise, there is no real trend here to observe. We could continue now and try and investigate INCIDENT_DATE further to see whether this correct. But there is much more to show in terms of spatial analysis, and we should move on. Before we do anything else, however, we first set back the colour palette to its default with palette(“default”).
palette("default")GIS work is about much more than simply visualising data on maps. We are, for instance, often interested how one spatial object relates to another. In our case, we are would like to know which incidents (points) occurred in a particular department (area) in Haiti (in Nord). There is a library in R called GISTools that offers great tools to deal with these kinds of questions. Load it with library(GISTools).
library(GISTools)Next we extract the spatial information about Nord from haiti_departments. This works in the same way as for any other data structure in R. Define haiti_department_nord <- subset(haiti_departments, ADM1==“Nord”).
haiti_department_nord <- subset(haiti_departments, ADM1=="Nord")Plot Nord with plot(haiti_department_nord).
plot(haiti_department_nord)
In order to find the incidents in Nord, we simply need to intersect haiti_department_nord with haiti_quake_sdf to clip all the relevant points. We can do this with haiti_department_nord_quake <- gIntersection(haiti_department_nord, haiti_quake_sdf). gIntersection is part of GISTools and intersects two geographies.
haiti_department_nord_quake <- gIntersection(haiti_department_nord, haiti_quake_sdf)Plot the result with plot(haiti_department_nord_quake, col = “red”, pch = 1, add = TRUE).
plot(haiti_department_nord)
plot(haiti_department_nord_quake, col = "red", pch = 1, add = TRUE)
Haiti is part of a larger island, which also includes the Dominican Republic. Its shape is loaded for you in dominican_republic. To recreate the island, we need to union the two geographies with gUnion of GISTools. Run island <- gUnion(haiti_departments, dominican_republic).
island <- gUnion(haiti_departments, dominican_republic)The result is the shape of the complete island. Run plot(island).
plot(island)
We talked earlier about incidents that seemed to be close to each and constitute location cluster. In the GIS world, one way to further investigation this claim is to rasterize a shape and see whether there are any clusters in the raster. This means we need to transfer values associated with spatial data objects (points, lines, polygons) to raster cells in a grid. To that end, we need to first load library(raster), which provides us with all we need for rasters - a data matrix, the ability to map information on a grid and add CRS information.
library(raster)According to the description of raster, we can create a spatial data structure that divides a region into rectangles called ‘cells’ that can store one or more values for each of these cells. Such a data structure is also referred to as a ‘grid’ and is often contrasted with ‘vector’ data that is used to represent points, lines, and polygons. If we find cells with the same or similar values, we know that the data is clustered in particular locations. We define a raster r with the extend of haiti_departments by running r <- raster(nrow = 130, ncol = 360, ext = extent(haiti_departments)).
r <- raster(nrow = 130, ncol = 360, ext = extent(haiti_departments))Running r <- rasterize(as(haiti_quake_sdf, “SpatialPoints”), r, fun = sum) will update the raster z with the spatial points of the Ushahidi incidents. With the attribute fun = sum, we tell the function to add up all the points per raster cell.
r <- rasterize(as(haiti_quake_sdf, "SpatialPoints"), r, fun = sum)We add haiti_departments with plot(haiti_departments, add = TRUE).
plot(island)
plot(haiti_departments, add = TRUE)
Finally, we realise the full raster image with plot(r, breaks=c(0, 20000, 40000, 60000, 90000), col=rev(brewer.pal(5, “Oranges”)), legend = FALSE, add = TRUE). If you look closely you can see a concentration of points in Port-au-Prince. But overall this is not yet a very good visualisation. We will later see how to develop better ones that show point patterns using Google Maps. Also check out the function kde.points in GISTools.
plot(island)
plot(haiti_departments, add = TRUE)
plot(r, breaks=c(0, 20000, 40000, 60000, 90000), col=rev(brewer.pal(5, "Oranges")), legend = FALSE, add = TRUE)
We have just done our first Point Pattern Analysis, where points are references for particular events, features, etc. In R, the standard package for Point Pattern Analysis is spatstat (https://cran.r-project.org/web/packages/spatstat/spatstat.pdf), which unfortunately needs its own data format (i.e. ppp). Fortunately, there are ways to transform a data frame into a ppp object. First we load library(spatstat).
library(spatstat)We start a new plot with plot(haiti_departments).
plot(haiti_departments)
In spatstat, the function owin is used to set the observation window. The package maptools provides a way to transform a spatial polygons into an object of class owin, using the function as.owin. Please, enter window <- as.owin(haiti_departments).
window <- as.owin(haiti_departments)The function ppp of spatstat creates the point pattern object. Define haiti_ppp <- ppp(x=haiti_quake_sdf@coords[,1],y = haiti_quake_sdf@coords[,2],window=window). Please, also note how you can access latitude and longitude in a spatial data frame using the @ sign. We have not used this yet, as we knew the coordinates already when we created haiti_quake_sdf. But in case you do not know the coordinates, this is the way to get them. To increase its correctness, we would have to clean the data further first, but we want to ignore this right now in order to concentrate on discussing the method.
haiti_ppp <- ppp(x=haiti_quake_sdf@coords[,1],y=haiti_quake_sdf@coords[,2],window=window)As said our assumption is that the Ushahidi incidents are unevenly distributed across Haiti. This spatial variation of the intensity can be identified using a method called quadrat counting, where the area is divided into rectangles (5x5 in this case) and the number of events in each of them is counted. To perform quadrat, please run plot(quadratcount(haiti_ppp, nx = 5, ny = 5), add=TRUE, col=“red”, alpha = 0.4) to visualise the counts for each cell and identify where they are clustered. We clearly see most reports come from Port au Prince, which we will investigate further next.
plot(haiti_departments)
plot(quadratcount(haiti_ppp, nx = 5, ny = 5), add=TRUE, col="red", alpha = 0.4)
I think we are all impressed now by the power of these GIS tools in R. At least, we should be. Next, let us introduce another popular method, which builds a choropleth map shading regions on the map according to their attributes. A choropleth map describes a map where parts are coloured according to some variable. You can simply type spplot(haiti_departments, z = “Population_density”) to map the population density per department. z defines the feature to prodcue the choropleth. Try it.
spplot(haiti_departments, z = "Population_density")
Choropleths are useful to visualise the distribution of attributes according to features in a spatial dataset, as just presented for Haiti’s population density. But we are of course interested in the 2010 earthquake and not general statistics of Haiti. To advance our understanding of the earthquake, a basic statistics counts the number of reported incidences per department in Haiti. GISTools contains poly.counts, which given a set of points, and a set of polygons, computes the number of points in each polygon. Define haiti_departments$hdq <- poly.counts(haiti_quake_sdf, haiti_departments).
haiti_departments$hdq <- poly.counts(haiti_quake_sdf, haiti_departments)With prop.table(haiti_departments$hdq/poly.areas(haiti_departments)) we can calculate the relative distribution of incidents. We can clearly see that the area around the capital dominates the distribution of incidents.
prop.table(haiti_departments$hdq/poly.areas(haiti_departments))## 0 1 2 3 4 5
## 0.730641402 0.085380427 0.032373017 0.002241893 0.037345657 0.010573328
## 6 7 8 9
## 0.037442909 0.008304957 0.034237125 0.021459286Next, we can plot the choropleth with spplot(haiti_departments, z=“hdq”, col.regions=brewer.pal(3,“Oranges”),colorkey=list(space=“bottom”), scales = list(draw = TRUE, col.regions=bpy.colors(3)), main=list(label=“Ushahidi Incidents per Department”,cex=2, fontfamily=“mono”), at = c(0,50,250,3000)). Take care with all the attributes, which create a rich visual experience.
spplot(haiti_departments, z="hdq", col.regions=brewer.pal(3,"Oranges"),colorkey=list(space="bottom"), scales = list(draw = TRUE, col.regions=bpy.colors(3)), main=list(label="Ushahidi Incidents per Department",cex=2, fontfamily="mono"), at = c(0,50,250,3000))
Despite the improvements, spplot is still quite limited. We can improve our choropleth visualisations easily by employing the tmap package, which uses the more advanced ggpplot visualisations. Load the library with library(tmap).
library(tmap)The tmap package is designed to map thematic data. Run tm_shape(haiti_departments) + tm_fill(“Population_density”, style=“quantile”, palette=“Reds”, title=“Population density (per square km)”) + tm_borders(“white”) + tm_legend(outside = TRUE, text.size = .8) + tm_layout(frame = FALSE) and simply admire the output. We will be going through the art of tmap productions in more detail later.
tm_shape(haiti_departments) + tm_fill("Population_density", style="quantile", palette="Reds", title="Population density \n(per square km)") + tm_borders("white") + tm_legend(outside = TRUE, text.size = .8) + tm_layout(frame = FALSE)
Much prettier! There is a lot going on here, but if you are used to ggplot the tmap specification is almost self-explanatory. The fill.style parameter defines the binning scheme. Here, we choose a quantile scheme, which means that the population density is evenly split up. Check out the help pages for the many different options. If you want to manually define your breaks, you can pass the option fixed to the fill.style parameter and then pass the desired breaks to the breaks parameter. For example, to break the data into 3 bins, run tm_shape(haiti_departments) + tm_fill(“Population_density”, style=“fixed”, breaks=c(0,230,350,1000 ), labels=c(“under 230”, “230 to 350”, “above 350”), palette=“Blues”, title=“Population density (per square km)”) + tm_borders(“grey”) + tm_legend(outside = TRUE, text.size = .8) + tm_layout(frame = FALSE). That is it for tmap right now. We will return to its abilities later.
tm_shape(haiti_departments) + tm_fill("Population_density", style="fixed", breaks=c(0,230,350,1000 ), labels=c("under 230", "230 to 350", "above 350"), palette="Blues", title="Population density \n(per square km)") + tm_borders("grey") + tm_legend(outside = TRUE, text.size = .8) + tm_layout(frame = FALSE)
We will come back later to tmap but for now we move on to study online sources of map data. While not too long ago it was difficult to get access to map data, nowadays we are of course used to have map data easily available on the web. One of the projects, which has made free online maps available is OpenStreetMap (https://www.openstreetmap.org/), an open source community resource. The R library OpenStreetMap provides an easy interface. Load it with library(OpenStreetMap). One of the projects, which has made this possible is OpenStreetMap (https://www.openstreetmap.org/), an open source community resource. The R library OpenStreetMap provides an easy interface. Load it with library(OpenStreetMap). This part is currently not active, as OpenStreetMap depends on Java on your system. The Java link to R can cause issues on Windows machines.
If you used OpenStreetMap before, you will know that it uses high-resolution map tiles (http://wiki.openstreetmap.org/wiki/Tiles). Tiles are typically 256×256 pixel images of a particular area that can be tiled together to get a full image of the geography you are interested in. To download the right tile from the OpenStreetMap repository, we need to define the bounding box (bbox) of spatial data for Haiti. Take a look at bbox(haiti_departments) and you will see four locations/points, the upper left and right corners as well as the lower left and right corners. Together they make up the tile square we are interested in.
OpenStreetMap needs the upper left corner and lower right corner. Let us define the upper left one (ul) first with ul <- as.vector(cbind(bbox(haiti_departments)[2,2],bbox(haiti_departments)[1,1])).
Next run lr <- as.vector(cbind(bbox(haiti_departments)[2,1],bbox(haiti_departments)[1,2])).
OpenStreetMaps’s function openmap downloads the right tiles for us. Type in open_haiti <- openmap(ul, lr, zoom = 9, type = “osm”) to retrieve the right tiles of Haiti. You need to be connected to the Internet. We found the correct zoom level of 9 by trial and error. You might want to leave swirl at this moment to experiment with different zoom levels. You might remember that you can do this with play() and then resume the lesson with nxt(). You can also go to https://www.openstreetmap.org/ and export the right parameters once you found the map you are looking for.
Type in plot(open_haiti) to plot open_haiti.
Now, we add the incidents with plot(spTransform(haiti_quake_sdf, osm()), add = TRUE, alpha = 0.1, col = “red”, pch=1). While most of the parameters are pretty standard, you have not learned about spTransform yet. It delivers conversions between projections. OpenStreet maps use the mercator projection (https://en.wikipedia.org/wiki/Mercator_projection), which is the standard for most map tiling systems. It can be easily derived form another projection and back. We will come back to the topic of transformation later.
The OpenStreetMap community was widely praised for its reaction to help with the crisis after the Haitian earthquake from 2010. Another big online map provider also helped. Google Maps were widely used to map incidents of the Haiti earthquake. The ggmap package makes it easy to retrieve raster map tiles from popular online mapping services like Google Maps but also others like OpenStreetmap. Load the library with library(ggmap).
library(ggmap)In particular, ggmap enables the integration of additional attributes and contextual information of geographies in the ggplot2 plotting framework, whose power we have already experienced. Type in library(ggplot2) to load the relevant library.
library(ggplot2)Just like for OpenStreetMap, we can download the tiles we are interested in. We noticed earlier a concentration of incidents in the capital of Haiti, Port-au-Prince. Let us load the Port-au-Prince tiles with haiti_gmap <- get_map(location = “Port-au-Prince”, zoom = 12).
haiti_gmap <- get_map(location = "Port-au-Prince", zoom = 12)## Source : https://maps.googleapis.com/maps/api/staticmap?center=Port-au-Prince&zoom=12&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-3VMyfPew## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Port-au-Prince&key=xxx-3VMyfPewUsing ggplot, we can simply map the incidents’ locations in Port-au-Prince with ggmap(haiti_gmap) + geom_point(aes(x = LONGITUDE, y = LATITUDE), colour = “red”, alpha = 0.1, size = 2, data = haiti_quake). You should be able to understand this expression. Otherwise, check online for the details of the ggplot.
ggmap(haiti_gmap) + geom_point(aes(x = LONGITUDE, y = LATITUDE), colour = "red", alpha = 0.1, size = 2, data = haiti_quake)
In Port-au-Prince locations seem to be concentrated around the Palais National. Earlier we tried to use raster images to identify these clusters. With ggplot, we can use geom_density2d to retrieve contours from a 2d density estimate and finally, with stat_density2d, calculate a hot spot layer and polygons with a fill colour based on relative frequency of points. ..level.. identifies the levels of the contours. Run the whole map plot with ggmap(haiti_gmap, extent = “device”) + geom_density2d(data = haiti_quake, aes(x = LONGITUDE, y = LATITUDE), size = 0.3) + stat_density2d(data = haiti_quake, aes(x = LONGITUDE, y = LATITUDE, fill = ..level.., alpha = ..level..), size = 0.01, bins = 16, geom = “polygon”) + scale_fill_gradient(low = “green”, high = “red”) + scale_alpha(range = c(0, 0.3), guide = FALSE). These are a lot of arguments. Try to ask your friend the Internet about the details and remember that you often do not have to understand the whole expression to reproduce it. In this case you just need to know how to address the right map (haiti_gmap) and feed the data (haiti_quake).
ggmap(haiti_gmap, extent = "device") + geom_density2d(data = haiti_quake, aes(x = LONGITUDE, y = LATITUDE), size = 0.3) + stat_density2d(data = haiti_quake, aes(x = LONGITUDE, y = LATITUDE, fill = ..level.., alpha = ..level..), size = 0.01, bins = 16, geom = "polygon") + scale_fill_gradient(low = "green", high = "red") + scale_alpha(range = c(0, 0.3), guide = FALSE)
This map clearly indicates a few hot spots across the city but one particular strong one around the Palais National. This concludes our introduction to spatial analysis with R and our work on the Haiti data.
Let’s check what we have learned next.
What does CRS stand for?
- Coordinate Reference Station
- Coordinate Reference System
- Canadian Research System
- Don’t know
Coordinate Reference System
Next we would like to plot in lightblue the two departments Sud and Sud-Est. Start the new plot by entering which command?
plot(haiti_departments)
Finish the plot by using subset to select the correct departments. Remember the colour is lightblue.
plot(haiti_departments)
plot(subset(haiti_departments, ADM1 %in% c("Sud", "Sud-Est")), col = "lightblue", add = TRUE)
Create a new plot of quadrat counting, where the area is divided into rectangles again, but 4x4 this time. Each quadrant should again contain a count of the number of earthquake events. Use add=TRUE, col=“red” and alpha = 0.4.
plot(haiti_departments)
plot(subset(haiti_departments, ADM1 %in% c("Sud", "Sud-Est")), col = "lightblue", add = TRUE)
plot(quadratcount(haiti_ppp, nx = 4, ny = 4), add=TRUE, col="red", alpha = 0.4)
Finally, let’s get a quick Google Map for London. Define lond_map using get_map with location = “London” and zoom = 16.
lond_map <- get_map(location = "London", zoom = 16)## Source : https://maps.googleapis.com/maps/api/staticmap?center=London&zoom=16&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-3VMyfPew## Source : https://maps.googleapis.com/maps/api/geocode/json?address=London&key=xxx-3VMyfPewNow plot lond_map with ggmap.
ggmap(lond_map)
Let’s move on to the second part of this session.