##Writing a function for spatial data: The Location Quotient.

##In some cases it is necessary to conduct the same analysis multiple times 
##on either the same or different data. In such circumstances it is worth 
##writing a function to simplify the code. In this example the location quotient 
##provides a simple calculation easily written in to a function. 

##Background: The location quotient (LQ) is an index for comparing a region's 
##share of a particular activity with the share of that same activity found at 
##a more aggregate spatial level. In this example we take a shapefile of London 
##Boroughs that contains information on the population of each borough and the
##percentage of sports participation in each borough. In this case there is 
##little point in calculating the LQ as the percentage alone would be more meaningful. 
##The focus here is how to undertake the methods, not their appropriate use, or the 
##validity of the results. 

Data Requirements:

London Sport Participation Shapefile.

Install the following packages (if you haven't already done so):

maptools, RcolourBrewer.

Code (Comments are preceded by ##)

## Load the packages
library(maptools)
library(RColorBrewer)

## Set the working directory (the folder containing your shapefile).

setwd("/XX/XX")

## load the shapefile

sport<- readShapePoly("london_sport.shp")

##have a look at the attribute table headings

names(sport)

## Extract the attribute data of interest (in this case its the % participation and borough population.

partic_per<-sport$Partic_Per

bor_pop<- sport$Pop_2001

## with these two columns the LQ can be calculated for each borough.  
##This is a useful calculation and we may wish to apply it over a number of 
##datasets to see temporal change. We can therefore write a function.

"LQ"<- function(partic_per, bor_pop) {

	sdp<- partic_per/bor_pop
	dp<-  sum(partic_per)/ sum(bor_pop)
	sdp/dp
}

## apply the function above:

bor_LQ<- data.frame(sport$ons_label, LQ(partic_per, bor_pop))
colnames(bor_LQ)<- c("ons_label", "bor_LQ")


## Using the resulting LQ values we can plot a map.

## Specify the colour breaks in the data (how many colours and what values the colour change at)

brks<-c(0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6)

## Use the RColorBrewer package to define a colour palette.
colours <- brewer.pal(6, "Blues") 


## Plot map

brks<-c(0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6)
colours <- brewer.pal(6, "Blues") 
plot(sport, col=colours[findInterval(bor_LQ$bor_LQ, brks,
 all.inside=TRUE)], axes=F)
box()


## Add map annotations etc...

SpatialPolygonsRescale(layout.north.arrow(1), offset= c(505100,164800), scale = 6000, plot.grid=F)
title(paste ("LQ of London Sports Participation"))

legend(x=503800, y=164800, legend=leglabs(brks), fill=colours, bty="n") 

##Done!


##Disclaimer: The methods provided here may not be the best solutions, just the ones I happen to 
##know about! No support is provided with these worksheets. I have tried to make them as 
##self-explanatory as possible and will not be able to respond to specific requests for help. 
##I do however welcome feedback on the tutorials. 
##License: cc-by-nc-sa. Contact: james@spatialanalysis.co.uk