Typographic conventions<\/h2>\n
To ensure reproducibility and allow automatic syntax highlighting, this document has been written in RMarkdown. Be aware of the following typographic conventions: R code (e.g. are compiled in-line: the\u00a0 For this tutorial you need to install R, the latest version of which can be downloaded from\u00a0http:\/\/cran.r-project.org\/<\/a>. A number of R editors such as\u00a0RStudio<\/a>\u00a0can be used to make R more user friendly, but these are not needed to complete the tutorial.<\/p>\n R has a huge and growing number of spatial data packages. These can be installed in one go with the\u00a0 The packages we will be using are\u00a0 All of the data used for the tutorial can be downloaded from here:<\/p>\n http:\/\/spatial.ly\/wp-content\/uploads\/2013\/12\/spatialggplot.zip<\/a><\/p>\n Save this to a new folder, then in R specify the path of that folder as you working directory. Use the\u00a0 If you are working in RStudio, you can create a project that will automatically set your working directory.<\/p>\n One of the most important steps in handling spatial data with R is the ability to read in shapefiles. There are a number of ways to do this. The most simple is\u00a0 This method works OK, but it is no longer considered best practice since it doesn\u2019t load in the spatial referencing information etc associated with the shapefile. A more powerful way to read in geographical data is to use the\u00a0 In the code above\u00a0 The file contains the borough population and the percentage of the population engaging in sporting activities and was taken from the\u00a0active people survey<\/a>. The boundary data is from the\u00a0Ordnance Survey<\/a>.<\/p>\n All shapefiles have an attribute table. This is loaded with\u00a0 R hides the geometry of spatial data unless you print the object (using the\u00a0 Type\u00a0 In the above code\u00a0 You will see you get a warning. This is simply saying that you are changing the coordinate reference system, not reprojecting the data. Epsg:27700 is the code for British National Grid. If we wanted to reproject the data into something like WGS84 for latitude and longitude we would use the following code:<\/p>\n The different epsg codes are a bit of hassle to remember but you can find them all at\u00a0spatialreference.org<\/a>.<\/p>\n This next section of the practical introduces a slightly different method of creating plots in R using the ggplot2 package. The package is an implementation of the Grammar of Graphics (Wilkinson 2005) – a general scheme for data visualization that breaks up graphs into semantic components such as scales and layers. ggplot2 can serve as a replacement for the base graphics in R (the functions you have been plotting with today) and contains a number of default options that match good visualisation practice.<\/p>\n The maps we produce will not be that meaningful – the focus here is on sound visualisation with R and not sound analysis (obviously the value of the former diminished in the absence of the latter!) Whilst the instructions are step by step you are encouraged to deviate from them (trying different colours for example) to get a better understanding of what we are doing.<\/p>\n Good examples of graphs can also be found on the website\u00a0cookbook-r.com<\/a>.<\/p>\n Load the package:<\/p>\n It is worth noting that the basic\u00a0 As a first attempt with ggplot2 we can create a scatter plot with the attribute data in the sport object created above. Type:<\/p>\n What you have just done is set up a ggplot object where you say where you want the input data to come from.\u00a0 Within the brackets you can alter the nature of the points. Try something like\u00a0 If you want to scale the points by borough population and colour them by sports participation this is also fairly easy by adding another\u00a0 The real power of ggplot2 lies in its ability to add layers to a plot. In this case we can add text to the plot.<\/p>\n This idea of layers (or geoms) is quite different from the standard plot functions in R, but you will find that each of the functions does a lot of clever stuff to make plotting much easier (see the documentation for a full list).<\/p>\n The following steps will create a map to show the percentage of the population in each London Borough who regularly participate in sports activities.<\/p>\n To get the shapefiles into a format that can be plotted we have to use the\u00a0 This step has lost the attribute information associated with the sport object. We can add it back using the merge function (this performs a data join). To find out how this function works look at the output of typing\u00a0 Take a look at the\u00a0 It is now straightforward to produce a map using all the built in tools (such as setting the breaks in the data) that ggplot2 has to offer.\u00a0 Now, just typing\u00a0 The default colours are really nice but we may wish to produce the map in black and white, which should produce a map like that shown below:<\/p>\n Saving plot images is also easy. You just need to use\u00a0 ggmap is a package that uses the ggplot2 syntax as a template to create maps with image tiles from the likes of Google and OpenStreetMap:<\/p>\n The sport object is in British National Grid but the ggmap image tiles are in WGS84. We therefore need to use the sport.wgs84 object created in the reprojection operation earlier.<\/p>\n The first job is to calculate the bounding box (bb for short) of the sport.wgs84 object to identify the geographic extent of the image tiles that we need.<\/p>\n This is then fed into the\u00a0 In much the same way as we did above we can then layer the plot with different geoms.<\/p>\n First fortify the sport.wgs84 object and then merge with the required attribute data (we already did this step to create the sport.f object).<\/p>\n We can now overlay this on our base map.<\/p>\n The code above contains a lot of parameters. Use the ggplot2 help pages to find out what they are. The resulting map looks okay, but it would be improved with a simpler base map in black and white. A design firm called stamen provide the tiles we need and they can be brought into the plot with the\u00a0 We can then produce the plot as before.<\/p>\n Finally, if we want to increase the detail of the base map, get_map has a zoom parameter.<\/p>\n This section builds on the previous information on plotting and highlights some of R’s more advanced spatial functions from the\u00a0 To reaffirm our starting point, let’s re-plot the only spatial dataset in our workspace, and count the number of polygons:<\/p>\n Because we are using borough-level data, and boroughs are official administrative zones, there is much data available at this level. We will use the example of crime data to illustrate this data availability, and join this with the current spatial dataset. As before, we can download and import the data from within R:<\/p>\n Initially, the\u00a0 Now that we have crime data at the borough level, the challenge is to join it by name. This is not always straightforward. Let us see which names in the crime data match the spatial data:<\/p>\n The first line of code above shows that all but one of the borough names matches; the second tells us that it is City of London that is named differently in the crime data. Look at the results (not shown here) on your computer.<\/p>\n The above code block first identified the row with the faulty name and then renamed the level to match the\u00a0 We are now ready to join the datasets. It is recommended to use the\u00a0 The documentation for join will be displayed if the plyr package is loaded (if not, load or install and load it!). It requires all joining variables to have the same name, so we will rename the variable to make the join work:<\/p>\n In addition to joining by zone name, it is also possible to do\u00a0spatial joins<\/a>\u00a0in R. There are three main varieties: many-to-one – where the values of many intersecting objects contribute to a new variable in the main table – one-to-many, or one-to-one. Because boroughs in London are quite large, we will conduct a many-to-one spatial join. We will be using Tube Stations as the spatial data to join, with the aim of finding out which and how many stations are found in each London borough.<\/p>\n The above code loads the data correctly, but also shows that there are problems with it: the Coordinate Reference System (CRS) differs from that of our shapefile. Although OSGB 1936 (or EPSG 27700) is the ‘correct’ CRS for the UK, we will convert the stations dataset into lat-long coordinates, as this is a more common CRS and enables easy base map creation:<\/p>\n Now we can clearly see that the stations overlay the boroughs. The problem is that the stations dataset is far more extensive than London borough dataset; we need<\/p>\n There are a number of functions that we can use to clip the points so that only those falling within London boroughs are retained:<\/p>\n We can write off the first one straight away as it is depreciated by the second. It seems that\u00a0 In the above code, only the first line actually ‘does’ anything in our workspace, by creating the object\u00a0 The first line instructs R to look at each column ( Now that we know how\u00a0plot(x, y)<\/code>) is written in a\u00a0monospace<\/code>\u00a0font while prose is not. Blocks of code such as,<\/p>\nc(1:3, 5)^2\n<\/code><\/pre>\n## [1] 1 4 9 25\n<\/code><\/pre>\n##<\/code>\u00a0indicates this is output from R. Some of the output from the code below is quite long; we only show the output that is useful. A single hash (#<\/code>) is a comment for humans to read that R will ignore. All images in this document are small and low-quality to save space; they should display better on your computer screen and can be saved at any resolution. The code presented here is not the only way to do things: we encourage you to play with it and try things out to gain a deeper understanding of R. Don’t worry, you cannot ‘break’ anything using R and all the input data can be re-loaded if things do go wrong.<\/p>\nPrerequisites and packages<\/h2>\n
ctv<\/code>\u00a0package and the command\u00a0install.views(\"Spatial\")<\/code>. We do NOT recommend running this command for this tutorial: partly because downloading and compiling all spatial packages takes a long time (hundreds of megabytes) and also because we will add new packages when they are needed to see what each does. We do recommend taking a quick browse at the range of spatial packages on offer though:\u00a0http:\/\/cran.r-project.org\/web\/views\/Spatial.html<\/a>.<\/p>\nggplot2<\/code>,\u00a0rgdal<\/code>,\u00a0rgeos<\/code>,\u00a0maptools<\/code>\u00a0and\u00a0ggmap<\/code>. To test whether ggplot2 is installed, for example, enter\u00a0library(ggpot2)<\/code>. If you get an error message, it needs to be installed:\u00a0install.packages(\"ggplot2\")<\/code>.<\/p>\nsetwd<\/code>\u00a0command to do this. If your username is \u201cusername\u201d and you saved the files into a folder called \u201crmapping\u201d on your Desktop, for example, you would type the following:<\/p>\nsetwd(\"C:\/Users\/username\/Desktop\/rmapping\/R\")\n<\/code><\/pre>\nLoading spatial data<\/h2>\n
readShapePoly()<\/code>\u00a0in the\u00a0maptools<\/code>\u00a0package:<\/p>\nlibrary(maptools) # load the package\nsport <- readShapePoly(\"london_sport.shp\") # read in the shapefile\n<\/code><\/pre>\nrgdal<\/code>\u00a0function\u00a0readOGR<\/code>, which automatically extracts this information. This is R\u2019s interface to the \u201cGeospatial Abstraction Library (GDAL)\u201d which is used by other open source GIS packages such as QGIS and enables R to handle a broader range of spatial data formats.<\/p>\nlibrary(rgdal)\nsport <- readOGR(dsn = \".\", \"london_sport\")\n<\/code><\/pre>\n## OGR data source with driver: ESRI Shapefile \n## Source: \".\", layer: \"london_sport\"\n## with 33 features and 4 fields\n## Feature type: wkbPolygon with 2 dimensions\n<\/code><\/pre>\ndsn<\/code>\u00a0is an\u00a0argument<\/em>\u00a0of the\u00a0function<\/em>\u00a0readOGR<\/code>. R functions have a default order of functions, so\u00a0dsn =<\/code>\u00a0does not actually need to be typed:\u00a0readOGR(\".\", \"london_sport\")<\/code>\u00a0works the same, but it is good to remember the meaning of each argument when beginning to use R, so we sometimes include argument names when it is relevant. Here,\u00a0dsn<\/code>\u00a0stands for \u201cdata source name\u201d which is the folder containing the spatial data \u2013 this was pre-specified when you set your working directory. The next argument is acharacter string<\/em>, text identifying the file required. There is no need to add a file extension.<\/p>\nreadOGR<\/code>\u00a0and can be treated in a similar way to an R\u00a0data frame<\/a>.<\/p>\nprint()<\/code>). Let’s take a look at the headings of sport, using the following command:\u00a0names(sport)<\/code>\u00a0The data contained in spatial data are kept in a ‘slot’ that can be accessed using the @ symbol:\u00a0sport@data<\/code>. This is useful if you do not wish to work with the spatial components of the data at all times.<\/p>\nsummary(sport)<\/code>\u00a0to get some additional information about the data object. Spatial objects in R contain a variety of additional information:<\/p>\nObject of class SpatialPolygonsDataFrame\nCoordinates:\n min max\nx 503571.2 561941.1\ny 155850.8 200932.5\nIs projected: TRUE \nproj4string :\n[+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 ....]\n<\/code><\/pre>\nproj4string<\/code>\u00a0represents the coordinate reference system used in the data. In this file it has been incorrectly specified so we can change it with the following:<\/p>\nproj4string(sport) <- CRS(\"+init=epsg:27700\")\n<\/code><\/pre>\n## Warning: A new CRS was assigned to an object with an existing CRS:\n## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs\n## without reprojecting.\n## For reprojection, use function spTransform in package rgdal\n<\/code><\/pre>\nsport.wgs84 <- spTransform(sport, CRS(\"+init=epsg:4326\"))\n<\/code><\/pre>\nggplot2<\/h1>\n
ggplot2<\/code>\u00a0is one of the best documented packages in R. The full documentation for it can be found online and it is recommended you test out the examples on your own machines and play with them:\u00a0http:\/\/docs.ggplot2.org\/current\/<\/a>\u00a0.<\/p>\nlibrary(ggplot2)\n<\/code><\/pre>\nplot()<\/code>\u00a0function requires no data preparation but additional effort in colour selection\/adding the map key etc.\u00a0qplot()<\/code>\u00a0and\u00a0ggplot()<\/code>\u00a0(from the ggplot2 package) require some additional steps to format the spatial data but select colours and add keys etc automatically. More on this later.<\/p>\np <- ggplot(sport@data, aes(Partic_Per, Pop_2001))\n<\/code><\/pre>\nsport@data<\/code>\u00a0is actually a data frame contained within the wider spatial object\u00a0sport<\/code>\u00a0(the\u00a0@<\/code>\u00a0enables you to access the attribute table of the sport shapefile). The characters inside the\u00a0aes<\/code>\u00a0argument refer to the parts of that data frame you wish to use (the variables\u00a0Partic_Per<\/code>\u00a0and\u00a0Pop_2001<\/code>). This has to happen within the brackets of\u00a0aes()<\/code>, which means, roughly speaking ‘aesthetics that vary’.
\nIf you just type p and hit enter you get the error\u00a0No layers in plot<\/code>. This is because you have not told ggplot what you want to do with the data. We do this by adding so-called \u201cgeoms\u201d, in this case\u00a0geom_point()<\/code>.<\/p>\np + geom_point()\n<\/code><\/pre>\n<\/p>\n
p + geom_point(colour = \"red\", size=2)<\/code>\u00a0and experiment.<\/p>\naes()<\/code>\u00a0argument.<\/p>\np + geom_point(aes(colour = Partic_Per, size = Pop_2001))\n<\/code><\/pre>\np + geom_point(aes(colour = Partic_Per, size = Pop_2001)) + geom_text(size = 2, \n aes(label = name))\n<\/code><\/pre>\n<\/p>\n
fortify()<\/code>\u00a0function. Spatial objects in R have a number of slots containing the various items of data (polygon geometry, projection, attribute information) associated with a shapefile. Slots can be thought of as shelves within the data object that contain the different attributes. The \u201cpolygons\u201d slot contains the geometry of the polygons in the form of the XY coordinates used to draw the polygon outline. The generic plot function can work out what to do with these, ggplot2 cannot. We therefore need to extract them as a data frame. The fortify function was written specifically for this purpose. For this to work, either\u00a0gpclib<\/code>\u00a0or\u00a0rgeos<\/code>\u00a0packages must be installed.<\/p>\n# library(gpclib); gpclibPermit() # uncomment if rgeos not installed\nsport.f <- fortify(sport, region = \"ons_label\")\n<\/code><\/pre>\n## Loading required package: rgeos\n## rgeos version: 0.2-19, (SVN revision 394)\n## GEOS runtime version: 3.3.8-CAPI-1.7.8 \n## Polygon checking: TRUE\n<\/code><\/pre>\n?merge<\/code>.<\/p>\nsport.f <- merge(sport.f, sport@data, by.x = \"id\", by.y = \"ons_label\")\n<\/code><\/pre>\nsport.f<\/code>\u00a0object to see its contents. You should see a large data frame containing the latitude and longitude (they are actually Easting and Northing as the data are in British National Grid format) coordinates alongside the attribute information associated with each London Borough. If you type\u00a0print(sport.f)<\/code>\u00a0you will just how many coordinate pairs are required! To keep the output to a minimum, take a peak at the object just using the\u00a0head<\/code>\u00a0command:<\/p>\nhead(sport.f[, 1:8])\n<\/code><\/pre>\n## id long lat order hole piece group name\n## 1 00AA 531027 181611 1 FALSE 1 00AA.1 City of London\n## 2 00AA 531555 181659 2 FALSE 1 00AA.1 City of London\n## 3 00AA 532136 182198 3 FALSE 1 00AA.1 City of London\n## 4 00AA 532946 181895 4 FALSE 1 00AA.1 City of London\n## 5 00AA 533411 182038 5 FALSE 1 00AA.1 City of London\n## 6 00AA 533843 180794 6 FALSE 1 00AA.1 City of London\n<\/code><\/pre>\ncoord_equal()<\/code>\u00a0is the equivalent of asp=T in regular plots with R:<\/p>\nMap <- ggplot(sport.f, aes(long, lat, group = group, fill = Partic_Per)) + geom_polygon() + \n coord_equal() + labs(x = \"Easting (m)\", y = \"Northing (m)\", fill = \"% Sport Partic.\") + \n ggtitle(\"London Sports Participation\")\n<\/code><\/pre>\nMap<\/code>\u00a0should result in your first ggplot-made map of London! There is a lot going on in the code above, so think about it line by line: what has each of the elements of code above has been designed to do. Also note how the\u00a0aes()<\/code>\u00a0components can be combined into one set of brackets after\u00a0ggplot<\/code>, that has relevance for all layers, rather than being broken into separate parts as we did above. The different plot functions still know what to do with these. The\u00a0group=group<\/code>\u00a0points ggplot to the group column added byfortify()<\/code>\u00a0and it identifies the groups of coordinates that pertain to individual polygons (in this case London Boroughs).<\/p>\nMap + scale_fill_gradient(low = \"white\", high = \"black\")\n<\/code><\/pre>\n<\/p>\n
ggsave<\/code>\u00a0after each plot, e.g.\u00a0ggsave(\"my_map.pdf\")<\/code>\u00a0will save the map as a pdf, with default settings. For a larger map, you could try the following:<\/p>\nggsave(\"my_large_plot.png\", scale = 3, dpi = 400)\n<\/code><\/pre>\nAdding base maps to ggplot2 with ggmap<\/h1>\n
library(ggmap) # you may have to use install.packages to install it first\n<\/code><\/pre>\nb <- bbox(sport.wgs84)\nb[1, ] <- (b[1, ] - mean(b[1, ])) * 1.05 + mean(b[1, ])\nb[2, ] <- (b[2, ] - mean(b[2, ])) * 1.05 + mean(b[2, ])\n# scale longitude and latitude (increase bb by 5% for plot) replace 1.05\n# with 1.xx for an xx% increase in the plot size\n<\/code><\/pre>\nget_map<\/code>\u00a0function as the location parameter. The syntax below contains 2 functions.\u00a0ggmap<\/code>\u00a0is required to produce the plot and provides the base map data.<\/p>\nlnd.b1 <- ggmap(get_map(location = b))\n<\/code><\/pre>\n## Warning: bounding box given to google - spatial extent only approximate.\n<\/code><\/pre>\nsport.wgs84.f <- fortify(sport.wgs84, region = \"ons_label\")\nsport.wgs84.f <- merge(sport.wgs84.f, sport.wgs84@data, by.x = \"id\", by.y = \"ons_label\")\n<\/code><\/pre>\nlnd.b1 + geom_polygon(data = sport.wgs84.f, aes(x = long, y = lat, group = group, \n fill = Partic_Per), alpha = 0.5)\n<\/code><\/pre>\nget_map<\/code>\u00a0function:<\/p>\nlnd.b2 <- ggmap(get_map(location = b, source = \"stamen\", maptype = \"toner\", \n crop = T))\n<\/code><\/pre>\nlnd.b2 + geom_polygon(data = sport.wgs84.f, aes(x = long, y = lat, group = group, \n fill = Partic_Per), alpha = 0.5)\n<\/code><\/pre>\nlnd.b3 <- ggmap(get_map(location = b, source = \"stamen\", maptype = \"toner\", \n crop = T, zoom = 11))\n\nlnd.b3 + geom_polygon(data = sport.wgs84.f, aes(x = long, y = lat, group = group, \n fill = Partic_Per), alpha = 0.5)\n<\/code><\/pre>\n<\/p>\n
Joining and clipping<\/h1>\n
rgeos<\/code>\u00a0package. We look at joining new datasets to our data – an attribute join – spatial joins, whereby data is added to the target layer depending on the location of the origins and clipping.<\/p>\nlibrary(rgdal)\nlnd <- readOGR(dsn = \".\", \"london_sport\")\n<\/code><\/pre>\n## OGR data source with driver: ESRI Shapefile \n## Source: \".\", layer: \"london_sport\"\n## with 33 features and 4 fields\n## Feature type: wkbPolygon with 2 dimensions\n<\/code><\/pre>\nplot(lnd)\n<\/code><\/pre>\n<\/p>\n
nrow(lnd)\n<\/code><\/pre>\n## [1] 33\n<\/code><\/pre>\nDownloading additional data<\/h2>\n
# download.file('http:\/\/data.london.gov.uk\/datafiles\/crime-community-safety\/mps-\n# recordedcrime-borough.csv', destfile = 'mps-recordedcrime-borough.csv')\n# uncomment and join the above code to download the data\n\ncrimeDat <- read.csv(\"mps-recordedcrime-borough.csv\") # flags an error\n<\/code><\/pre>\nread.csv<\/code>\u00a0command flags an error: open the raw .csv file in a text editor such as Notepad, Notepad++ or GVIM, find the problem and correct it. Alternatively, you can work out what the file encoding is and use the correct argument (this is not recommended – simpler just to edit the text file in most cases).<\/p>\ncrimeDat <- read.csv(\"mps-recordedcrime-borough.csv\", fileEncoding = \"UCS-2LE\")\nhead(crimeDat)\nsummary(crimeDat$MajorText)\ncrimeTheft <- crimeDat[which(crimeDat$MajorText == \"Theft & Handling\"), ]\nhead(crimeTheft, 2) # change 2 for more rows\ncrimeAg <- aggregate(CrimeCount ~ Spatial_DistrictName, FUN = \"sum\", data = crimeTheft)\nhead(crimeAg, 2) # show the aggregated crime data\n<\/code><\/pre>\nlnd$name %in% crimeAg$Spatial_DistrictName\n<\/code><\/pre>\n## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\n## [12] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE\n## [23] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE\n<\/code><\/pre>\nlnd$name[which(!lnd$name %in% crimeAg$Spatial_DistrictName)]\n<\/code><\/pre>\n## [1] City of London\n## 33 Levels: Barking and Dagenham Barnet Bexley Brent Bromley ... Westminster\n<\/code><\/pre>\nlevels(crimeAg$Spatial_DistrictName)\nlevels(crimeAg$Spatial_DistrictName)[25] <- as.character(lnd$name[which(!lnd$name %in% \n crimeAg$Spatial_DistrictName)])\nlnd$name %in% crimeAg$Spatial_DistrictName # now all columns match\n<\/code><\/pre>\nlnd<\/code>\u00a0dataset. Note that we could not rename the variable directly, as it is stored as a factor.<\/p>\njoin<\/code>\u00a0function in the\u00a0plyr<\/code>\u00a0package but the\u00a0merge<\/code>\u00a0function could equally be used.<\/p>\nhelp(join)\nlibrary(plyr)\nhelp(join) # now help should appear\n<\/code><\/pre>\nhead(lnd$name)\nhead(crimeAg$Spatial_DistrictName) # the variables to join\ncrimeAg <- rename(crimeAg, replace = c(Spatial_DistrictName = \"name\"))\nhead(join(lnd@data, crimeAg)) # test it works\n<\/code><\/pre>\n## Joining by: name\n<\/code><\/pre>\nlnd@data <- join(lnd@data, crimeAg)\n<\/code><\/pre>\n## Joining by: name\n<\/code><\/pre>\nAdding point data for clipping and spatial join<\/h2>\n
download.file(\"http:\/\/www.personal.leeds.ac.uk\/~georl\/egs\/lnd-stns.zip\", \"lnd-stns.zip\")\nunzip(\"lnd-stns.zip\")\nlibrary(rgdal)\nstations <- readOGR(dsn = \".\", layer = \"lnd-stns\", p4s = \"+init=epsg:27700\")\nproj4string(stations) # this is the full geographical detail.\nproj4string(lnd)\nbbox(stations)\nbbox(lnd)\n<\/code><\/pre>\nstationsWGS <- spTransform(stations, CRSobj = CRS(proj4string(lnd)))\nstations <- stationsWGS\nrm(stationsWGS)\nplot(lnd)\npoints(stations[sample(1:nrow(stations), size = 500), ])\n<\/code><\/pre>\n<\/p>\n
Clipping<\/h2>\n
?overlay\n?sp::over\nlibrary(rgeos)\n?rgeos::gIntersects\n<\/code><\/pre>\ngIntersects<\/code>\u00a0can produce the same output as\u00a0over<\/code>, based on\u00a0discussion<\/a>\u00a0in the community, so either can be used. (See this\u00a0discussion<\/a>\u00a0for further alternatives.) In this tutorial we will use\u00a0gIntersects<\/code>, for clipping although we could equally use\u00a0gContains<\/code>,gWithin<\/code>\u00a0and other\u00a0g...<\/code>\u00a0functions – see rgeos help pages by typing\u00a0?gOverlaps<\/code>\u00a0or other functions for more.\u00a0gIntersects<\/code>\u00a0will output information for each point, telling us which polygon it interacts with (i.e. the polygon it is in):<\/p>\nint <- gIntersects(stations, lnd, byid = T) # find which stations intersect \nclass(int) # it's outputed a matrix\ndim(int) # with 33 rows (one for each zone) and 2532 cols (the points)\nsummary(int[, c(200, 500)]) # not the output of this\nplot(lnd)\npoints(stations[200, ], col = \"red\") # note point id 200 is outside the zones\npoints(stations[500, ], col = \"green\") # note point 500 is inside\nwhich(int[, 500] == T) # this tells us that point 500 intersects with zone 32\npoints(coordinates(lnd[32, ]), col = \"black\") # test the previous statement\n<\/code><\/pre>\n<\/p>\n
int<\/code>. The proceeding lines are dedicated to exploring this object and what it means. Note that it is a matrix with columns corresponding to the points and rows corresponding to boroughs. The borough in which a particular point can be extracted from\u00a0int<\/code>\u00a0as we shall see below. For the purposes of clipping, we are only interested in whether the point intersects with\u00a0any<\/em>\u00a0of the boroughs. This is where the function\u00a0apply<\/code>, which is unique to R, comes into play:<\/p>\nclipped <- apply(int == F, MARGIN = 2, all)\nplot(stations[which(clipped), ]) # shows all stations we DO NOT want\nstations.cl <- stations[which(!clipped), ] # use ! to select the invers\npoints(stations.cl, col = \"green\") # check that it's worked\n<\/code><\/pre>\n<\/p>\n
stations <- stations.cl\nrm(stations.cl) # tidy up: we're only interested in clipped ones\n<\/code><\/pre>\nMARGIN = 2<\/code>, we would use\u00a0MARGIN = 1<\/code>\u00a0for row-by-row analysis) and report back whether\u00a0all<\/code>\u00a0of the values are false. This creates the inverse selection that we want, hence the use of\u00a0!<\/code>\u00a0to invert it. We test that the function works on a new object (often a good idea, to avoid overwriting useful data) with plots and, once content that the clip has worked, save the sample of points to our main\u00a0stations<\/code>\u00a0object and remove the now duplicated\u00a0stations.cl<\/code>\u00a0object.<\/p>\nAggregating the data to complete the spatial join<\/h2>\n
gIntersects<\/code>\u00a0works in general terms and for clipping, let’s use it to allocate a borough to each of our station points, which we will then aggregate up. Data from these points (e.g. counts, averages in each area etc.) can then be transferred to the main polygons table: the essence of a spatial join. Again,\u00a0apply<\/code>\u00a0is our friend in this instance, meaning we can avoid\u00a0for<\/code>\u00a0loops:<\/p>\n