---
title:  "Mapping Environmental Action in Scotland"
abstract:    
# thanks: "Replication files are available on the author's Github account (https://github.com/kidwellj/mapping_environmental_action). **Current version**: `r format(Sys.time(), '%B %d, %Y')`
style:  jeremy1
author: "[Jeremy H. Kidwell](http://jeremykidwell.info)"
affiliation: University of Birmingham
institute: University of Birmingham
e-mail: "[j.kidwell@bham.ac.uk](mailto:j.kidwell@bham.ac.uk)"
date: "`r Sys.Date()`"
bibliography: biblio.bib
linkcolor: black
geometry: margin=1in
# fontfamily: mathpazo
fontsize: 11pt
output:
  html_document:
    theme: readable
    keep_md: true
    code_folding: hide
    self_contained: true
    toc: true
    toc_depth: 2
    number_sections: true
    fig_caption: true
    fig_retina: 1
  pdf_document:
    toc: false
    keep_tex: true
    number_sections: true
    fig_caption: true
    citation_package: natbib
    latex_engine: xelatex
  always_allow_html: yes
  
---

```{r setup, include=FALSE}
require(knitr)
require(kableExtra)
knitr::opts_chunk$set(fig.path='figures/', warning=FALSE, echo=FALSE, message=FALSE, dpi=300, fig.width=7)
# TODO: consider implementing knitcitations - https://github.com/cboettig/knitcitations
# TODO: fix simultaneous output towards PDF, see here: https://stackoverflow.com/questions/23621012/display-and-save-the-plot-simultaneously-in-r-rstudio
```

```{r load_packages, message=FALSE, warning=FALSE, include=FALSE}
## Default repo
# setwd("/Users/jeremy/gits/mapping_environmental_action")
# setwd("/Users/kidwellj/OneDrive\ -\ bham.ac.uk/writing/201708_mapping_environmental_action")

# Set repository to be new standard, e.g. cloud server. 
# This will avoid a dialogue box if packages are to be installed for below on first run.
local({r <- getOption("repos")
       r["CRAN"] <- "https://cloud.r-project.org" 
       options(repos=r)
})
# TODO: remove sp etc. once sf is fully implemented
# TODO: automatically test for packages below on given execution environment and run install.packages() as needed.
require(RCurl) # used for fetching reproducible datasets
require(sf) # new simplefeature data class, supercedes sp in many ways
require(sp) # needed for proj4string, deprecated by sf()
require(rgdal) # deprecated by sf()
require(GISTools) # deprecated by sf()
require(rgeos) # deprecated by sf()
require(maptools)
require(ggplot2)
require(tmap) # using as an alternative to base r graphics and ggplot for geospatial plots
require(tmaptools) # for get_asp_ratio below
require(grid) # using for inset maps on tmap
require(broom) # required for tidying SPDF to data.frame for ggplot2
require(tidyr)  # using for grouped bar plot
require(plyr)
require(dplyr)
require(reshape2)  # using for grouped bar plot
require(scales)
# require(sqldf) # using sqldf to filter before loading very large data sets
require(plotly) # allows for export of plots to dynamic web pages
require(gtable) # more powerful package for multi-plot layouts, not necessary for knitr
require(showtext) # for loading in fonts
require(extrafont) # font support

# Set up local workspace:
if (dir.exists("data") == FALSE) {
  dir.create("data") 
}
if (dir.exists("figures") == FALSE) {
  dir.create("figures") 
}
if (dir.exists("derivedData") == FALSE) {
  dir.create("derivedData")
}

# Define Coordinate Reference Systems (CRS) for later use
# Note: I've used British National Grid (27000) in this paper, but have found that 
# it is falling out of use in many cases, so will be defaulting to WGS84 in future 
# data-sets and papers.

# Working with EPSG codes for spatialfeature CRS given the usage of this approach with sf()
bng <- CRS("+init=epsg:27700")
wgs84 <- CRS("+init=epsg:4326")

# Configure fonts for plots below

## Loading Google fonts (http://www.google.com/fonts)
# font_add_google("Merriweather", "merriweather")
# The following will load in system fonts (uncomment and run as needed on first execution)
# font_import(pattern="[A/a]rial", prompt=FALSE)
```

# Introduction[^15541312]

Until recently, environmentalism has been treated by governments and environmental charities as a largely secular concern. In spite of the well-developed tradition of "eco-theology" which began in earnest in the UK in the mid-twentieth century (and which has many precursors in previous centuries), third-sector groups and governments, particularly in Britain and Europe, have largely ignored religious groups as they have gone about their business crafting agendas for behaviour change, developing funding programmes, and developing platforms to mitigate ecological harm, motivate consumers and create regulation regimes. That this has changed is evidenced by the fact that several prominent non-religious environmental groups have commissioned studies and crafted outreach programmes to persons with a particular faith tradition or to "spiritual communities" including RSPB (2013) and the Sierra Club USA (2008).[^158261118] Further, since 2008, the Scottish Government has provided a significant portion of funding for the ecumenical charity, Eco-Congregation Scotland, which works to promote literacy on environmental issues in religious communities in Scotland and helps to certify congregations under their award programme. What is not well known, however, even by these religious environmental groups themselves, is whether or how their membership might be different from other environmental groups. This study represents an attempt to illuminate this new interest with some more concrete data about religious groups in Scotland and how they may differ from non-religious counterparts.

# Eco-Congregation Scotland: The Basics

```{r load_ecs_data, message=FALSE, warning=FALSE}
# read in Eco-Congregation Scotland data and-------------------
# ...turn it into a SpatialPointsDataFrame---------------------
# TODO: upload ECS-GIS-Locations_3.0.csv to zenodo repository, i.e.
ecs <- read.csv("data/ECS-GIS-Locations_3.0.csv", comment.char="#")
# unnecessary with advent of sf (above)
coordinates(ecs) <- c("X", "Y")
# Modified to use EPSG code directly 27 Feb 2019
proj4string(ecs) <- bng
ecs_sf <- st_as_sf(ecs, coords = c("X", "Y"), crs=paste0("+init=epsg:",27700))
```

There are `r length(ecs)` eco-congregations in Scotland. By some measurements, particularly in terms of individual sites and possibly also with regards to volunteers, this makes Eco-Congregation Scotland one of the largest environmental third-sector groups in Scotland.[^159141043] 

In seeking to conduct GIS and statistical analysis of ECS, it is important to note that there some ways in which these sites are statistically opaque. Our research conducted through interviews at a sampling of sites and analysis of a variety of documents suggests that there is a high level of diversity both in terms of the number of those participating in environmental action and the types of action underway at specific sites. Work at a particular site can also ebb and flow over the course of time. Of course, as research into other forms of activism and secular environmental NGOs has shown, this is no different from any other third sector volunteer group. Variability is a regular feature of groups involved in activism and/or environmental concern. 

For the sake of this analysis, we took each Eco-Congregation Scotland site to represent a point of analysis as if each specific site represented a community group which had "opted-in" on environmental concern. On this basis, in this section, in the tradition of human geography, we "map" environmental action  among religious communities in Scotland a variety of ways. This is the first major geographical analysis of this kind conducted to date in Europe. We measure the frequency and location of ECS sites against a variety of standard geo-referenced statistical data sets, seeking to provide a statistical and geographically based assessment of the participation of religious groups in relation to the following:

- Location within Scotland
- Religious affiliation
- Relation to the Scottish Index of Multiple Deprivation (SIMD)
- Relation to the 8-Fold Scottish Government Urban-Rural Scale 
- Proximity to "wilderness" (based on several different designations)

For the sake of comparison, we also measured the geographical footprint of two other forms of community group in Scotland, (1) Transition Towns (taking into account their recent merge with Scotland Communities Climate Action Network) and (2) member groups of the Development Trust Association Scotland ("DTAS"). These two groups provide a helpful basis for comparison as they are not centralised and thus have a significant geographical dispersion across Scotland. They also provide a useful comparison as transition is a (mostly) non-religious environmental movement, and community development trusts are not explicitly linked to environmental conservation (though this is often part of their remit), so we have a non-religious point of comparison in Transition and a non-environmental point of comparison with DTAS

# Technical Background

Analysis was conducted using QGIS 2.8 and R `r getRversion()`, and data-sets were generated in CSV format.[^15541313] To begin with, I assembled a data set consisting of x and y coordinates for each congregation in Scotland and collated this against a variety of other specific data. Coordinates were checked by matching UK postcodes of individual congregations against geo-referencing data in the Office for National Statistics postcode database. In certain instances a single "congregation" is actually a series of sites which have joined together under one administrative unit. In these cases, each site was treated as a separate data point if worship was held at that site at least once a month, but all joined sites shared a single unique identifier. As noted above, two other datasets were generated for the sake of comparative analysis.[^177171536] These included one similar Environmental Non-Governmental Organisation (ENGO) in Scotland (1) Transition Scotland (which includes Scotland Communities Climate Action Network);[^15541342] and another community-based NGO, Scottish Community Development Trusts.[^158261232] As this report will detail, these three overlap in certain instances both literally and in terms of their aims, but each also has a separate identity and footprint in Scotland. Finally, in order to normalise data, we utilised the PointX POI dataset which maintains a complete database of Places of Worship in Scotland.[^15541614]

# Background and History of Eco-Congregation Scotland

Eco-Congregation Scotland began a year before the official launch of Eco-Congregation England and Wales, in 1999, as part of an effort by Kippen Environment Centre (later renamed to Forth Environment Link, or "FEL") a charity devoted to environmental education in central Scotland[^158261210] to broaden the scope of its environmental outreach to churches in central Scotland.[^15826124] Initial funding was provided, through Kippen Environment Centre by way of a "sustainable action grant" (with funds drawn from a government landfill tax) through a government programme called Keep Scotland Beautiful (the Scottish cousin of Keep Britain Tidy). After this initial pilot project concluded, the Church of Scotland provided additional funding for the project in the form of staff time and office space. Additional funding a few years later from the Scottish Government helped subsidise the position of a business manager, and in 2011 the United Reformed Church contributed additional funding which subsidised the position of a full-time environmental chaplain for a 5-year term, bringing the total staff to five.

```{r calculate_ecs_by_year, message=FALSE, warning=FALSE}
# Tidy up date fields and convert to date data type
ecs$registration <- as.Date(ecs$registration, "%Y-%m-%d")
# TODO: Fix issues here with R complaining that "character string is not in a std...
# ecs$award1 <- as.Date(ecs$award1)
# ecs$award2 <- as.Date(ecs$award2)
# ecs$award3 <- as.Date(ecs$award3)
# ecs$award4 <- as.Date(ecs$award4)
# TODO: add "R" to command in paragraph below once this is resolved, do search for all non 'r instances
ecs_complete_cases <- ecs[complete.cases(ecs$year_begun),]
```

The programme launched officially in 2001 at Dunblane Cathedral in Stirling and by 2005 the project had `r length(ecs_complete_cases[ecs_complete_cases$year_begun < 2006, ])` congregations registered to be a part of the programme and 25 which had completed the curriculum successfully and received an Eco-Congregation award. By 2011, the number of registrations had tripled to `r length(ecs_complete_cases[ecs_complete_cases$year_begun < 2012, ])` and the number of awarded congregations had quadrupled to `sum(ecs$award1 < "01/01/2012", na.rm=TRUE)`. This process of taking registrations and using a tiered award or recognition scheme is common to many voluntary organisations. The ECS curriculum was developed in part by consulting the Eco-Congregation England and Wales materials which had been released just a year earlier in 1999, though it has been subsequently revised, particularly with a major redesign in 2010. In the USA, a number of similar groups take a similar approach including Earth Ministry (earthministry.org) and Green Faith (greenfaith.org). 

In the case of Eco-Congregation Scotland, congregations are invited to begin by "registering" their interest in the programme by completing a basic one-sided form. The next step requires the completion of an award application, which includes a facilitated curriculum called a "church check-up" and after an application is submitted, the site is visited and assessed by third-party volunteer assessors. Sites are invited to complete additional applications for further awards which are incremental (as is the application process). Transition communities, at least in the period reflected on their map, go through a similar process (though this does not involve the use of a supplied curriculum) by which they are marked first as "interested," become "active" and then gain "official" status.[^1554162]

# Representation by Regional Authorities (Council Areas) {.tabset}

```{r import_admin_data, message=FALSE, warning=FALSE, include=FALSE}
# read in polygon for Scottish admin boundaries
# TODO: upload bundle of admin data to new zenodo repository and alter below to use new URLs
# TODO: need to remove readOGR below once st_read is confirmed to be working as sf

if (file.exists("data/scotland_ca_2010.shp") == FALSE) {
download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/Scotland_ca_2010.zip", 
              destfile = "data/Scotland_ca_2010.zip")
unzip("data/Scotland_ca_2010.zip", exdir = "data")
}
admin_lev1 <- readOGR("./data", "scotland_ca_2010")
admin_lev1_sf <- st_read("data/scotland_ca_2010.shp") %>% st_transform(paste0("+init=epsg:",27700))

# read in polygon for intermediate admin boundary layers
if (file.exists("data/scotland_parlcon_2011.shp") == FALSE) {
download.file("http://census.edina.ac.uk/ukborders/easy_download/prebuilt/shape/Scotland_parlcon_2011.zip", 
              destfile = "data/Scotland_parlcon_2011.zip")
unzip("data/Scotland_parlcon_2011.zip", exdir = "data")
}
admin_lev2 <- readOGR("./data", "scotland_parlcon_2011")
admin_lev2_sf <- st_read("data/scotland_parlcon_2011.shp") %>% st_transform(paste0("+init=epsg:",27700))

# Set CRS using epsg code on spdf for symmetry with datasets below
proj4string(admin_lev1) <- bng
proj4string(admin_lev2) <- bng

# Generate new sf shape using bounding box for central belt for map insets below
# Note: coordinates use BNG as CRS (EPSG: 27700)

scotland <- st_bbox(c(xmin = 5513.0000, xmax = 470332.0000,
                      ymin = 530249.0000, ymax = 1220301.5000),
                    crs = st_crs("+init=epsg:27700")) %>% 
  st_as_sfc()
  
centralbelt_region <- st_bbox(c(xmin = 234841, xmax = 346309,
                      ymin = 653542, ymax = 686722),
                    crs = st_crs("+init=epsg:27700")) %>% 
  st_as_sfc()

centralbelt_ratio <- get_asp_ratio(centralbelt_region)
scotland_ratio<- get_asp_ratio(scotland)
```

```{r import_groups_data, message=FALSE, warning=FALSE, include=FALSE}
# read in Transition Towns data and turn it into a SpatialPointsDataFrame
transition_wgs <- read.csv(text=getURL("https://zenodo.org/record/165519/files/SCCAN_1.4.csv"))
coordinates(transition_wgs) <- c("X", "Y")
proj4string(transition_wgs) <- CRS("+init=epsg:4326")
transition <- spTransform(transition_wgs, bng)
transition_sf <- st_as_sf(transition, coords = c("X", "Y"), crs=27700)

# read in all_churches data (data set generated by Jeremy Kidwell to replace PointX data used from Ordnance Survey)
# TODO: need to remove all data points which are outside BNG area to resolve error
# also need to make symmetrical with ECS denominations, add Methodist 
# churches, remove nazarene and salvation army
 
# churches_all <- read.csv("data/all_churches_0.9.csv")
# churches_all_clean<-churches_all[complete.cases(churches_all),]
# churches_all_null<-churches_all[!complete.cases(churches_all),]
# coordinates(churches_all) <- c("X", "Y")
# proj4string(churches_all) <- proj4string(admin_lev1)

# read in pointX data and turn it into a SpatialPointsDataFrame
pow_pointX <- read.csv("./data/poi_2015_12_scot06340459.csv", sep="|")
coordinates(pow_pointX) <- c("feature_easting", "feature_northing")
# TODO: need to alter to draw from wgs84 or bng defined in preamble above
proj4string(pow_pointX) <- CRS("+init=epsg:27700")
pow_pointX_sf <- st_as_sf(pow_pointX, coords = c("X", "Y"), crs=27700)

# read in Scottish Community Dev. trust data and turn it into a SpatialPointsDataFrame
dtas <- read.csv("data/community-dev-trusts-2.6.csv")
coordinates(dtas) <- c("X", "Y")
proj4string(dtas) <- CRS("+init=epsg:27700")
dtas_sf <- st_as_sf(dtas, coords = c("X", "Y"), crs=27700)

# read in permaculture data and turn it into a SpatialPointsDataFrame
permaculture <- read.csv("data/permaculture_scot-0.8.csv")
coordinates(permaculture) <- c("X", "Y")
proj4string(permaculture) <- CRS("+init=epsg:27700")
permaculture_sf <- st_as_sf(permaculture, coords = c("X", "Y"), crs=27700)

# subset point datasets for inset maps below
ecs_sf_centralbelt <- st_intersection(ecs_sf, centralbelt_region)
dtas_sf_centralbelt <- st_intersection(dtas_sf, centralbelt_region)
transition_sf_centralbelt <- st_intersection(transition_sf, centralbelt_region)
pow_sf_centralbelt <- st_intersection(pow_pointX_sf, centralbelt_region)
permaculture_sf_centralbelt <- st_intersection(permaculture_sf, centralbelt_region)

```

```{r process_admin_data}
# Augment existing dataframes to run calculations and add columns with point counts per polygon, 
# percentages, and normalising data.

# This code will generate a table of frequencies for each spatialpointsdataframe in admin
# calculate count of ECS for fields in admin and provide percentages
# JK Note: need to convert from poly.counts, which uses sp() to st_covers() which uses sf() - cf. https://stackoverflow.com/questions/45314094/equivalent-of-poly-counts-to-count-lat-long-pairs-falling-inside-of-polygons-w#45337050

admin_lev1$ecs_count <- poly.counts(ecs,admin_lev1)
admin_lev1$ecs_percent<- prop.table(admin_lev1$ecs_count)
# calculate count of places of worship in PointX db for fields in admin and provide percentages
admin_lev1$pow_count <- poly.counts(pow_pointX,admin_lev1)
admin_lev1$pow_percent<- prop.table(admin_lev1$pow_count)
# calculate count of Transition for fields in admin and provide percentages
admin_lev1$transition_count <- poly.counts(transition,admin_lev1)
admin_lev1$transition_percent<- prop.table(admin_lev1$transition_count)
# calculate count of dtas for fields in admin and provide percentages
admin_lev1$dtas_count <- poly.counts(dtas,admin_lev1)
admin_lev1$dtas_percent<- prop.table(admin_lev1$dtas_count)
# calculate count of permaculture for fields in admin and provide percentages
admin_lev1$permaculture_count <- poly.counts(permaculture,admin_lev1)
admin_lev1$permaculture_percent<- prop.table(admin_lev1$permaculture_count)

# run totals for intermediate boundaries level 2
# This code will generate a table of frequencies for each spatialpointsdataframe in admin_lev2
# calculate count of ECS for fields in admin_lev2 and provide percentages
admin_lev2$ecs_count <- poly.counts(ecs,admin_lev2)
admin_lev2$ecs_percent<- prop.table(admin_lev2$ecs_count)
# calculate count of places of worship in PointX db for fields in admin_lev2 and provide percentages
admin_lev2$pow_count <- poly.counts(pow_pointX,admin_lev2)
admin_lev2$pow_percent<- prop.table(admin_lev2$pow_count)
# calculate count of Transition for fields in admin_lev2 and provide percentages
admin_lev2$transition_count <- poly.counts(transition,admin_lev2)
admin_lev2$transition_percent<- prop.table(admin_lev2$transition_count)
# calculate count of dtas for fields in admin_lev2 and provide percentages
admin_lev2$dtas_count <- poly.counts(dtas,admin_lev2)
admin_lev2$dtas_percent<- prop.table(admin_lev2$dtas_count)
# calculate count of permaculture for fields in admin_lev2 and provide percentages
admin_lev2$permaculture_count <- poly.counts(permaculture,admin_lev2)
admin_lev2$permaculture_percent<- prop.table(admin_lev2$permaculture_count)

# calculate count of ECS for fields in admin_lev2_sf
# TODO: for future migration to sf throughout, remove above content and swap out references.
admin_lev2_sf$ecs_count <- lengths(st_covers(admin_lev2_sf, st_as_sf(ecs_sf, coords = c("long", "lat"), crs = st_crs(27700))))
admin_lev2_sf$pow_count <- lengths(st_covers(admin_lev2_sf, st_as_sf(pow_pointX_sf, coords = c("long", "lat"), crs = st_crs(27700))))
admin_lev2_sf$transition_count <- lengths(st_covers(admin_lev2_sf, st_as_sf(transition_sf, coords = c("long", "lat"), crs = st_crs(27700))))
admin_lev2_sf$dtas_count <- lengths(st_covers(admin_lev2_sf, st_as_sf(dtas_sf, coords = c("long", "lat"), crs = st_crs(27700))))
admin_lev2_sf$permaculture_count <- lengths(st_covers(admin_lev2_sf, st_as_sf(permaculture_sf, coords = c("long", "lat"), crs = st_crs(27700))))

# Import csv with population data for each level of administrative subdivision and join to spatialdataframe

# Placeholder for England parish data for later implementation
# download.file("http://census.edina.ac.uk/ukborders/easy_download/prebuilt/shape/England_cp_1991.zip", destfile = "parishes/parishes-1991.zip")
# unzip("parishes/parishes-1991.zip", exdir = "parishes")
# parishes <- rgdal::readOGR(dsn = "parishes", "england_cp_1991")

# TODO - consider adapting to use ONS mid-year statistics: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/populationestimatesforukenglandandwalesscotlandandnorthernireland

# Load population statistics for normalising data by population and merge with admin_lev1
admin_lev1_pop <- read.csv("./data/scotland_admin_2011pop.csv")
admin_lev1 <- merge(x=admin_lev1, y=admin_lev1_pop, by.x = "code", by.y = "CODE")
# Convert number stored as string with a comma to an integer
admin_lev1$X2011_pop <- as.numeric(gsub(",", "", admin_lev1@data$X2011_pop))
# Calculate counts as percentages of total for normalising below
admin_lev1$pop_percent<- prop.table(admin_lev1$X2011_pop)
admin_lev1$pow_percent<- prop.table(admin_lev1$pow_count)
# Normalise ecs_count using population figures (using ArcGIS method)
admin_lev1$ecs_count_popnorm <- admin_lev1$ecs_count * admin_lev1$pop_percent
# Normalise ecs_count using places of worship counts (using ArcGIS method)
admin_lev1$ecs_count_pownorm <- admin_lev1$ecs_count * admin_lev1$pow_percent

# Preserve scale
admin_lev1$ecs_count_popnorm_scaled <- admin_lev1$ecs_count_popnorm*(sum(admin_lev1$ecs_count)/sum(admin_lev1$ecs_count_popnorm))
admin_lev1$ecs_count_pownorm_scaled <- admin_lev1$ecs_count_pownorm*(sum(admin_lev1$ecs_count)/sum(admin_lev1$ecs_count_pownorm))

# Load population statistics for normalising data by population on admin_lev2
# TODO: source this data file in a replicable way (e.g. zenodo)
admin_lev2_pop <- read.csv("./data/scotland_and_wales_const_scotland_2011pop.csv")
admin_lev2 <- merge(x=admin_lev2, y=admin_lev2_pop, by.x = "code", by.y = "CODE")
# Convert number stored as string with a comma to an integer (using data.zone data)
admin_lev2$Data.zone.Population <- as.numeric(gsub(",", "", admin_lev2$Data.zone.Population)) 
# Calculate counts as percentages of total for normalising below
admin_lev2$pop_percent <- prop.table(admin_lev2$Data.zone.Population)
admin_lev2$pow_percent <- prop.table(admin_lev2$pow_count)
# Normalise ecs_count using population figures (using ArcGIS method)
admin_lev2$ecs_count_popnorm <- admin_lev2$ecs_count * admin_lev2$pop_percent
# Normalise ecs_count using places of worship counts (using ArcGIS method)
admin_lev2$ecs_count_pownorm <- admin_lev2$ecs_count * admin_lev2$pow_percent

# Preserve scale
admin_lev2$ecs_count_popnorm_scaled <- admin_lev2$ecs_count_popnorm*(sum(admin_lev2$ecs_count)/sum(admin_lev2$ecs_count_popnorm))
admin_lev2$ecs_count_pownorm_scaled <- admin_lev2$ecs_count_pownorm*(sum(admin_lev2$ecs_count)/sum(admin_lev2$ecs_count_pownorm))
```

```{r additional_poi_data}
# Load in updated DTAS data set
dtas_new <- read.csv("data/dtas_4.0.csv")
coordinates(dtas_new) <- c("X", "Y")
proj4string(dtas_new) <- CRS("+init=epsg:27700")

# Load in retail data from geolytics dataset
# from here: https://geolytix.co.uk/?retail_points
poi_grocery_wgs <- read.csv("data/retailpoints_version11_dec17.txt", sep = "\t")
# select useful columns
poi_grocery_wgs <- subset(poi_grocery_wgs, select = c("retailer", "store_name", "long_wgs", "lat_wgs"))
# convert to spdf
coordinates(poi_grocery_wgs) <- c("long_wgs", "lat_wgs")
proj4string(poi_grocery_wgs) <- CRS("+init=epsg:4326")
poi_grocery <- spTransform(poi_grocery_wgs, proj4string(admin_lev1))
# filter out non-Scottish data
poi_grocery <- poi_grocery[!is.na(over(poi_grocery, geometry(admin_lev1))),]
poi_grocery_sf <- st_as_sf(poi_grocery, coords = c("long_wgs", "lat_wgs"), crs=paste0("+init=epsg:",27700))

# Load in British pubs from Ordnance survey dataset
poi_pubs <- read.csv("data/poi_pubs.csv", header = FALSE, sep = "|")
# select useful columns 
poi_pubs <- subset(poi_pubs, select = c("V1", "V2", "V3", "V4", "V5"))
# rename columns to tidier names
colnames(poi_pubs) <- c("refnum", "name", "code", "x", "y")
coordinates(poi_pubs) <- c("x", "y")
proj4string(poi_pubs) <- proj4string(admin_lev1)
# filter out non-Scottish pubs
poi_pubs <- poi_pubs[!is.na(over(poi_pubs, geometry(admin_lev1))),]
poi_pubs_sf <- st_as_sf(poi_pubs, coords = c("x", "y"), crs=paste0("+init=epsg:",27700))
```

```{r wilderness_data_prep} 

# 1. Download data for SSSI:

if (file.exists("data/SSSI_SCOTLAND.shp") == FALSE) {
# TODO: upload data to zenodo, uncomment below
# http://gateway.snh.gov.uk/natural-spaces/dataset.jsp?dsid=SSSI
# download.file("", destfile = "data/SSSI_SCOTLAND_ESRI.zip")
unzip("data/SSSI_SCOTLAND_ESRI.zip", exdir = "data")
}
sssi <- st_read("data/SSSI_SCOTLAND.shp") %>% st_transform(paste0("+init=epsg:",27700))
sssi_sp <- readOGR("./data", "SSSI_SCOTLAND") 
# Generate simplified polygon for plots below
sssi_simplified <- st_simplify(sssi)
# sssi_simplified_sp <- rgeos::gSimplify(sssi_sp, tol=3)

# 2. Download wild land areas:

if (file.exists("data/WILDLAND_SCOTLAND.shp") == FALSE) {
# TODO: upload data to zenodo, uncomment below
# https://gateway.snh.gov.uk/natural-spaces/dataset.jsp?dsid=WILDLAND
# download.file("", destfile = "data/WILDLAND_SCOTLAND_ESRI.zip")
unzip("data/WILDLAND_SCOTLAND_ESRI.zip", exdir = "data")
}

wildland <- st_read("data/WILDLAND_SCOTLAND.shp") %>% st_transform(paste0("+init=epsg:",27700))
# Generate simplified polygon for plots below
wildland_simplified <- st_simplify(wildland)

# 3. Download data for National Forest Inventory:

# Note: UK-wide data is here: https://opendata.arcgis.com/datasets/bcd6742a2add4b68962aec073ab44138_0.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D

if (file.exists("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") == FALSE) {
download.file("https://opendata.arcgis.com/datasets/3cb1abc185a247a48b9d53e4c4a8be87_0.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D", 
              destfile = "data/National_Forest_Inventory_Woodland_Scotland_2017.zip")
unzip("data/National_Forest_Inventory_Woodland_Scotland_2017.zip", exdir = "data")
}

forestinv <- st_read("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") %>% st_transform(paste0("+init=epsg:",27700))
# Generate simplified polygon for plots below
forestinv_simplified <- st_simplify(forestinv)

# Download data for scenic areas
# https://opendata.arcgis.com/datasets/d7f6b987c7224a72a185ce012258d500_23.zip
# International Union for Conservation of Nature (IUCN), Category V Protected Landscapes
# for England: https://environment.data.gov.uk/DefraDataDownload/?mapService=NE/AreasOfOutstandingNaturalBeautyEngland&Mode=spatial

# download.file("https://opendata.arcgis.com/datasets/d7f6b987c7224a72a185ce012258d500_23.zip", destfile = "data/ScenicAreas.zip")
# unzip("data/ScenicAreas.zip", exdir = "data")

scenicareas <- st_read("data/SG_NationalScenicAreas_1998.shp") %>% st_transform(paste0("+init=epsg:",27700))
# Generate simplified polygon for plots below
scenicareas_simplified <- st_simplify(scenicareas)

# Set symmetrical CRS for analysis below (inserted here in order to correct errors, may be deprecated later)
# st_crs(sssi) <- 27700
# st_crs(ecs_sf) <- 27700
# st_crs(pow_pointX_sf) <- 27700
# st_crs(dtas_sf) <- 27700
# st_crs(transition_sf) <- 27700
# st_crs(permaculture_sf) <- 27700

# Define buffer and measure number of ECS groups within 0.5 miles of all SSSI
# CRS uses meters for units, so a buffer for 0.5 miles would use 805 meters)

# Define buffers as lines and polygons for reuse below (as some of these operations take > 10 minutes)
# TODO: Working with simplified polygons here to optimise execution, need to confirm calculations are the same as polygons without simplification (likely!)

sssi_buf50 <- st_buffer(sssi_simplified, dist = 50)
sssi_buf500 <- st_buffer(sssi_simplified, dist = 500)
# Lines offer even more optimised representation suitable for plotting
# TODO: resolve error:
# Error in st_cast.sfc(., to = "LINESTRING") : 
#  use smaller steps for st_cast; first cast to MULTILINESTRING or POLYGON?
# sssi_buf50_lines <- st_union(sssi_simplified) %>% st_buffer(50) %>% 
#   st_cast(to = "MULTILINESTRING") %>% 
#   st_cast(to = "LINESTRING")
# sssi_buf500_lines <- st_union(sssi_simplified) %>% st_buffer(500) %>% 
#   st_cast(to = "MULTILINESTRING") %>% 
#   st_cast(to = "LINESTRING")

wildland_buf50 <- st_buffer(wildland_simplified, dist = 50)
wildland_buf500 <- st_buffer(wildland_simplified, dist = 500)
# Lines offer even more optimised representation suitable for plotting
# wildland_buf50_lines = st_union(wildland_simplified) %>% st_buffer(50) %>% 
#   st_cast(to = "LINESTRING")
# wildland_buf500_lines = st_union(wildland_simplified) %>% st_buffer(500) %>% 
#   st_cast(to = "LINESTRING")

forestinv_buf50 <- st_buffer(forestinv_simplified, dist = 50)
forestinv_buf500 <- st_buffer(forestinv_simplified, dist = 500)
# Lines offer even more optimised representation suitable for plotting
# forestinv_buf50_lines = st_union(forestinv_simplified) %>% st_buffer(50) %>% 
#   st_cast(to = "LINESTRING")
# forestinv_buf500_lines = st_union(forestinv_simplified) %>% st_buffer(500) %>% 
#   st_cast(to = "LINESTRING")

scenicareas_buf50 <- st_buffer(scenicareas_simplified, dist = 50)
scenicareas_buf500 <- st_buffer(scenicareas_simplified, dist = 500)

# Calculate number of groups within polygons

# calculate coincidence of ecs points within each polygons and buffers for each
# TODO: possibly use st_difference(sssi_buf50, sssi)
# TODO: consider subsetting instead here, e.g.:

# plot(lnd, col = "lightgrey") # plot the london_sport object
# sel <- lnd$Partic_Per > 25
# plot(lnd[ sel, ], col = "turquoise", add = TRUE) # add selected zones to map
# from https://gotellilab.github.io/Bio381/StudentPresentations/SpatialDataTutorial.html

ecs_sf_sssi <- st_within(ecs_sf, sssi_simplified)
ecs_sf_sssi50m <- st_within(ecs_sf, sssi_buf50)
ecs_sf_sssi500m <- st_within(ecs_sf, sssi_buf500)
ecs_sf_sssibeyond500m <- !(st_within(ecs_sf, sssi_buf500))

ecs_sf_wildland <- st_within(ecs_sf, wildland_simplified)
ecs_sf_wildland50m <- st_within(ecs_sf, wildland_buf50)
ecs_sf_wildland500m <- st_within(ecs_sf, wildland_buf500)
ecs_sf_wildlandbeyond500m <- !(st_within(ecs_sf, wildland_buf500))

ecs_sf_forestinv <- st_within(ecs_sf, forestinv_simplified)
ecs_sf_forestinv50m <- st_within(ecs_sf, forestinv_buf50)
ecs_sf_forestinv500m <- st_within(ecs_sf, forestinv_buf500)
ecs_sf_forestinvbeyond500m <- !(st_within(ecs_sf, forestinv_buf500))

ecs_sf_scenicareas <- st_within(ecs_sf, scenicareas_simplified)
ecs_sf_scenicareas50m <- st_within(ecs_sf, scenicareas_buf50)
ecs_sf_scenicareas500m <- st_within(ecs_sf, scenicareas_buf500)
ecs_sf_scenicareasbeyond500m <- !(st_within(ecs_sf, scenicareas_buf500))

# TODO: implement more efficient code using do.call() function or sapply() as here https://stackoverflow.com/questions/3642535/creating-an-r-dataframe-row-by-row

# Generate dataframe based on SSSI buffers

# Calculate incidence of ecs within SSSI and within buffers at 50/500m
ecs_sssi_row <- c(sum(apply(st_within(ecs_sf, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, sssi_buf500, sparse=FALSE), 1, any)))

pow_sssi_row <- c(sum(apply(st_within(pow_pointX_sf, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, sssi_buf500, sparse=FALSE), 1, any)))

dtas_sssi_row <- c(sum(apply(st_within(dtas_sf, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, sssi_buf500, sparse=FALSE), 1, any)))

transition_sssi_row <- c(sum(apply(st_within(transition_sf, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, sssi_buf500, sparse=FALSE), 1, any)))

permaculture_sssi_row <- c(sum(apply(st_within(permaculture_sf, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, sssi_buf500, sparse=FALSE), 1, any)))

grocery_sssi_row <- c(sum(apply(st_within(poi_grocery_sf, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, sssi_buf500, sparse=FALSE), 1, any)))

pubs_sssi_row <- c(sum(apply(st_within(poi_pubs_sf, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, sssi_buf500, sparse=FALSE), 1, any)))


# Generate dataframe from rows based on counts
sssi_counts <- rbind(ecs_sssi_row, pow_sssi_row)
sssi_counts <- rbind(sssi_counts, dtas_sssi_row)
sssi_counts <- rbind(sssi_counts, transition_sssi_row)
sssi_counts <- rbind(sssi_counts, permaculture_sssi_row)
sssi_counts <- rbind(sssi_counts, grocery_sssi_row)
sssi_counts <- rbind(sssi_counts, pubs_sssi_row)
sssi_counts <- as.data.frame(sssi_counts)
colnames(sssi_counts) <- c("Within SSSIs", "...50m", "...500m")

# Generate dataframe from rows based on percentages of totals
# TODO: Work with ecs_sf etc., but first sort out why length is different between shapefiles
ecs_sssi_row_pct <- ecs_sssi_row/length(ecs)
pow_sssi_row_pct <- pow_sssi_row/length(pow_pointX)
dtas_sssi_row_pct <- dtas_sssi_row/length(dtas)
transition_sssi_row_pct <- transition_sssi_row/length(transition)
permaculture_sssi_row_pct <- permaculture_sssi_row/length(permaculture)
grocery_sssi_row_pct <- grocery_sssi_row/length(poi_grocery_sf)
pubs_sssi_row_pct <- pubs_sssi_row/length(poi_pubs_sf)

sssi_counts_pct <- rbind(ecs_sssi_row_pct, pow_sssi_row_pct)
sssi_counts_pct <- rbind(sssi_counts_pct, dtas_sssi_row_pct)
sssi_counts_pct <- rbind(sssi_counts_pct, transition_sssi_row_pct)
sssi_counts_pct <- rbind(sssi_counts_pct, permaculture_sssi_row_pct)
sssi_counts_pct <- rbind(sssi_counts_pct, grocery_sssi_row_pct)
sssi_counts_pct <- rbind(sssi_counts_pct, pubs_sssi_row_pct)
colnames(sssi_counts_pct) <- c("% Within SSSIs", "% within 50m", "% within 500m")

# Merge into larger dataframe
sssi_counts_merged <- cbind(sssi_counts, sssi_counts_pct)

# Generate dataframe based on wildland buffers

ecs_wildland_row <- c(sum(apply(st_within(ecs_sf, wildland_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, wildland_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, wildland_buf500, sparse=FALSE), 1, any)))

pow_wildland_row <- c(sum(apply(st_within(pow_pointX_sf, wildland_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, wildland_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, wildland_buf500, sparse=FALSE), 1, any)))
wildland_counts <- rbind(ecs_wildland_row, pow_wildland_row)

dtas_wildland_row <- c(sum(apply(st_within(dtas_sf, wildland_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, wildland_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, wildland_buf500, sparse=FALSE), 1, any)))
wildland_counts <- rbind(wildland_counts, dtas_wildland_row)

transition_wildland_row <- c(sum(apply(st_within(transition_sf, wildland_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, wildland_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, wildland_buf500, sparse=FALSE), 1, any)))
wildland_counts <- rbind(wildland_counts, transition_wildland_row)

permaculture_wildland_row <- c(sum(apply(st_within(permaculture_sf, wildland_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, wildland_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, wildland_buf500, sparse=FALSE), 1, any)))
wildland_counts <- rbind(wildland_counts, permaculture_wildland_row)

grocery_wildland_row <- c(sum(apply(st_within(poi_grocery_sf, wildland_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, wildland_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, wildland_buf500, sparse=FALSE), 1, any)))

pubs_wildland_row <- c(sum(apply(st_within(poi_pubs_sf, wildland_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, wildland_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, wildland_buf500, sparse=FALSE), 1, any)))


colnames(wildland_counts) <- c("Within Wildland Areas", "...50m", "...500m")

# Generate dataframe from rows based on percentages of totals

ecs_wildland_row_pct <- ecs_wildland_row/length(ecs_sf)
pow_wildland_row_pct <- pow_wildland_row/length(pow_pointX)
dtas_wildland_row_pct <- dtas_wildland_row/length(dtas)
transition_wildland_row_pct <- transition_wildland_row/length(transition)
permaculture_wildland_row_pct <- permaculture_wildland_row/length(permaculture)

wildland_counts_pct <- rbind(ecs_wildland_row_pct, pow_wildland_row_pct)
wildland_counts_pct <- rbind(wildland_counts_pct, dtas_wildland_row_pct)
wildland_counts_pct <- rbind(wildland_counts_pct, transition_wildland_row_pct)
wildland_counts_pct <- rbind(wildland_counts_pct, permaculture_wildland_row_pct)
colnames(wildland_counts_pct) <- c("% Within wildlands", "% within 50m", "% within 500m")

# Merge into larger dataframe
wildland_counts_merged <- cbind(wildland_counts, wildland_counts_pct)

# Generate dataframe based on forestinv buffers

ecs_forestinv_row <- c(sum(apply(st_within(ecs_sf, forestinv_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, forestinv_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, forestinv_buf500, sparse=FALSE), 1, any)))

pow_forestinv_row <- c(sum(apply(st_within(pow_pointX_sf, forestinv_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, forestinv_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, forestinv_buf500, sparse=FALSE), 1, any)))
forestinv_counts <- rbind(ecs_forestinv_row, pow_forestinv_row)

dtas_forestinv_row <- c(sum(apply(st_within(dtas_sf, forestinv_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, forestinv_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, forestinv_buf500, sparse=FALSE), 1, any)))
forestinv_counts <- rbind(forestinv_counts, dtas_forestinv_row)

transition_forestinv_row <- c(sum(apply(st_within(transition_sf, forestinv_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, forestinv_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, forestinv_buf500, sparse=FALSE), 1, any)))
forestinv_counts <- rbind(forestinv_counts, transition_forestinv_row)

permaculture_forestinv_row <- c(sum(apply(st_within(permaculture_sf, forestinv_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, forestinv_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, forestinv_buf500, sparse=FALSE), 1, any)))
forestinv_counts <- rbind(forestinv_counts, permaculture_forestinv_row)

grocery_forestinv_row <- c(sum(apply(st_within(poi_grocery_sf, forestinv_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, forestinv_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, forestinv_buf500, sparse=FALSE), 1, any)))

pubs_forestinv_row <- c(sum(apply(st_within(poi_pubs_sf, forestinv_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, forestinv_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, forestinv_buf500, sparse=FALSE), 1, any)))

colnames(forestinv_counts) <- c("Within Woodlands", "...50m", "...500m")

# Generate dataframe from rows based on percentages of totals
ecs_forestinv_row_pct <- ecs_forestinv_row/length(ecs_sf)
pow_forestinv_row_pct <- pow_forestinv_row/length(pow_pointX)
dtas_forestinv_row_pct <- dtas_forestinv_row/length(dtas)
transition_forestinv_row_pct <- transition_forestinv_row/length(transition)
permaculture_forestinv_row_pct <- permaculture_forestinv_row/length(permaculture)

forestinv_counts_pct <- rbind(ecs_forestinv_row_pct, pow_forestinv_row_pct)
forestinv_counts_pct <- rbind(forestinv_counts_pct, dtas_forestinv_row_pct)
forestinv_counts_pct <- rbind(forestinv_counts_pct, transition_forestinv_row_pct)
forestinv_counts_pct <- rbind(forestinv_counts_pct, permaculture_forestinv_row_pct)
colnames(forestinv_counts_pct) <- c("% Within Woodlands", "% within 50m", "% within 500m")

# Merge into larger dataframe
forestinv_counts_merged <- cbind(forestinv_counts, forestinv_counts_pct)

# Generate dataframe based on scenicareas buffers

ecs_scenicareas_row <- c(sum(apply(st_within(ecs_sf, scenicareas_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, scenicareas_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, scenicareas_buf500, sparse=FALSE), 1, any)))

pow_scenicareas_row <- c(sum(apply(st_within(pow_pointX_sf, scenicareas_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, scenicareas_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, scenicareas_buf500, sparse=FALSE), 1, any)))
scenicareas_counts <- rbind(ecs_scenicareas_row, pow_scenicareas_row)

dtas_scenicareas_row <- c(sum(apply(st_within(dtas_sf, scenicareas_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, scenicareas_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, scenicareas_buf500, sparse=FALSE), 1, any)))
scenicareas_counts <- rbind(scenicareas_counts, dtas_scenicareas_row)

transition_scenicareas_row <- c(sum(apply(st_within(transition_sf, scenicareas_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, scenicareas_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, scenicareas_buf500, sparse=FALSE), 1, any)))
scenicareas_counts <- rbind(scenicareas_counts, transition_scenicareas_row)

permaculture_scenicareas_row <- c(sum(apply(st_within(permaculture_sf, scenicareas_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, scenicareas_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, scenicareas_buf500, sparse=FALSE), 1, any)))
scenicareas_counts <- rbind(scenicareas_counts, permaculture_scenicareas_row)

grocery_scenicareas_row <- c(sum(apply(st_within(poi_grocery_sf, scenicareas_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, scenicareas_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_grocery_sf, scenicareas_buf500, sparse=FALSE), 1, any)))

pubs_scenicareas_row <- c(sum(apply(st_within(poi_pubs_sf, scenicareas_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, scenicareas_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(poi_pubs_sf, scenicareas_buf500, sparse=FALSE), 1, any)))

colnames(scenicareas_counts) <- c("Within Scenic Areas", "...50m", "...500m")

# Generate dataframe from rows based on percentages of totals

ecs_scenicareas_row_pct <- ecs_scenicareas_row/length(ecs_sf)
pow_scenicareas_row_pct <- pow_scenicareas_row/length(pow_pointX)
dtas_scenicareas_row_pct <- dtas_scenicareas_row/length(dtas)
transition_scenicareas_row_pct <- transition_scenicareas_row/length(transition)
permaculture_scenicareas_row_pct <- permaculture_scenicareas_row/length(permaculture)

scenicareas_counts_pct <- rbind(ecs_scenicareas_row_pct, pow_scenicareas_row_pct)
scenicareas_counts_pct <- rbind(scenicareas_counts_pct, dtas_scenicareas_row_pct)
scenicareas_counts_pct <- rbind(scenicareas_counts_pct, transition_scenicareas_row_pct)
scenicareas_counts_pct <- rbind(scenicareas_counts_pct, permaculture_scenicareas_row_pct)
colnames(scenicareas_counts_pct) <- c("% Within scenicareass", "% within 50m", "% within 500m")

# Merge into larger dataframe
scenicareas_counts_merged <- cbind(scenicareas_counts, scenicareas_counts_pct)


```

# Proximity to "Wilderness"

Chasing down a curiosity, I decided to try and calculate whether proximity to "wilderness" or "scenic nature" or just trees might have some impact on generating more mobilised communities. I realised that there would be several problems with this kind of calculation up front, first being that "nature" is a deeply problematic construct, reviled by geographers and philosophers alike. With this in mind, I identified several different ways of reckoning wilderness, starting with the highly anachronistic "Scenic Land" designation from the 1970s. Then I pursued the more carefully calculated "core wild areas" generated by SNH just a few years ago. However, even the core wile areas concept has been criticised heavily, so I also expanded out my search to include all sites of special scientific interest and then went even wider to include the Scottish Forestry Service's "Native Woodland" and finally, the most generic possible measurement, any land identified as forested at the last Forest Inventory.

Proximity to these areas was the next concern, because many of these designations deliberately exclude human habitat, so it was necessary to measure the number of sites within proximity. There is a question which lies here regarding aesthetics, namely, what sort of proximity might generate an affective connection? From my own experience, I decided upon the distance represented by a short walk, i.e. a half-kilometre. However, with the more generic measurements, such as SSSI and forestation, this wouldn't do, as there are so many of these sites that a buffer of 500 meters encapsulates almost all of inhabited Scotland. So for these sites I also calculated a count within 50 metres.

So what did I discover? The results were inconclusive. First, it is important to note that on the whole, Eco-Congregations tend to be more urban than place of worship taken generally at a rate of nearly 3:1 (5.4% of Eco-Congregations lie in areas currently designated as "Very Remote Rural Areas" whereas nearly 15% of places of worship lie in these areas), so what I was testing for was whether this gap was smaller when specifying these various forms of "wild" remoteness. For our narrowest measurements, there were so few sites captured as to render measurement unreliable. There are, for obvious reasons, `st_within(ecs_sf, wildland)` sites located within any of SNG's core wild areas. Similarly, there are very few of our activist communities located within SSSI's (only `st_within(pow_pointX_sf, sssi)` places of worship out of `r length(pow_pointX)`, `st_within(transition_sf, sssi)` transition towns, (or 2%) and `st_within(dtas_sf, sssi)` community development trusts (3%)). However, expanding this out makes things a bit more interesting, within 50 metres of SSSI's in Scotland lie `st_within(ecs_sf, sssi_buf50)` Eco-Congregations (or just under 1%), which compares favourably with the `st_within(pow_pointX_sf, sssi_buf50)` places of worship (or just 1.5%) far exceeding our ratio (1:1.5 vs. 1:3). This is the same with our more anachronistic measure of "scenic areas," there are 7 eco-congregations within these areas, and 175 places of worship, making for a ratio of nearly 1:2 (2.1% vs. 4.3%). Taking our final measure, of forested areas, this is hard to calculate, as only `st_within(ecs_sf, forestinv)` Eco-Congregation lies within either native or generally forested land.


```{r 13_wilderness_tables} 

# Output mmd tables using kable

sssi_counts_merged %>%
  kable(format = "html", col.names = colnames(sssi_counts_merged), caption = "Group counts within SSSIs") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", full_width = F, "responsive"))

wildland_counts_merged %>%
  kable(format = "html", col.names = colnames(wildland_counts_merged), caption = "Group counts within Wildland Areas") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", full_width = F, "responsive"))

forestinv_counts_merged %>%
  kable(format = "html", col.names = colnames(forestinv_counts_merged), caption = "Group counts within Woodlands") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", full_width = F, "responsive"))

scenicareas_counts_merged %>%
  kable(format = "html", col.names = colnames(scenicareas_counts_merged), caption = "Group counts within Scenic Areas") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", full_width = F, "responsive"))


# Output CSV files for tables above
write.csv(sssi_counts_merged, "derivedData/sssi_counts_merged.csv", row.names=TRUE)
write.csv(wildland_counts_merged, "derivedData/wildland_counts_merged.csv", row.names=TRUE)
write.csv(forestinv_counts_merged, "derivedData/forestinv_counts_merged.csv", row.names=TRUE)
write.csv(scenicareas_counts_merged, "derivedData/scenicareas_counts_merged.csv", row.names=TRUE)

```

# Appendix A

```{r admin_table}
# Output CSV files for various levels of admin
write.csv(admin_lev1, "derivedData/admin_lev1.csv", row.names=FALSE)
write.csv(admin_lev2, "derivedData/admin_lev2.csv", row.names=FALSE)
write.csv(ecs, "derivedData/ecs.csv", row.names=FALSE)
write.csv(transition, "derivedData/transition.csv", row.names=FALSE)
write.csv(permaculture, "derivedData/permaculture.csv", row.names=FALSE)
write.csv(dtas, "derivedData/dtas.csv", row.names=FALSE)
write.csv(simd, "derivedData/simd.csv", row.names=FALSE)

# Output mmd tables using kable
admin_lev1_df <- as_data_frame(admin_lev1[,c(3,5,7,11,13)])

admin_lev1_df %>%
  kable(format = "html", col.names = colnames(admin_lev1_df)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
    kableExtra::scroll_box(width = "100%", height = "800px")
```

# Appendix B

```{r admin_table_withproportions}

admin_lev1_prop_df <- as_data_frame(admin_lev1[,c(3,5:14,22:25)])

admin_lev1_prop_df %>%
  kable(format = "html", col.names = colnames(admin_lev1_prop_df)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
    kableExtra::scroll_box(width = "100%", height = "800px")
```

# Appendix C - Data by Urban / Rural Classification

```{r urbanrural_table}
# Output CSV files for urbanrural tables
write.csv(urbanrural, "derivedData/urbanrural.csv", row.names=FALSE)

# Output mmd tables using kable

urbanrural_table <- as_data_frame(urbanrural[,c(1,7:16)])

urbanrural_table %>%
  kable(format = "html", col.names = colnames(urbanrural_table)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
    kableExtra::scroll_box(width = "100%")

```

# Citations

[^15541312]: This research was jointly funded by the AHRC/ESRC under project numnbers AH/K005456/1 and AH/P005063/1.
[^158261118]: This is not to say that there have been no collaborations before 2000, noteworthy in this respect is the WWF who helped to found the Alliance of Religion and Conservation (ARC) in 1985.
[^159141043]: This suggestion should be qualified - RSPB would greatly exceed ECS both in terms of the number of individual subscribers and budget. The RSPB trustee's report for 2013-2014 suggests that their member base was 1,114,938 people across Britain with a net income of £127m - the latter of which exceeds the Church of Scotland. If we adjust this based on the Scottish share of the population of the United Kingdom as of the 2011 census (8.3%) this leaves us with an income of £9.93m. The British charity commission requires charities to self-report the number of volunteers and staff, and from their most recent statistics we learn that RSPB engaged with 17,600 volunteers and employed 2,110 members of staff. Again, adjusted for population, this leaves 1,460 volunteers in Scotland and 176 staff. However, if we measure environmental groups based on the number of sites they maintain, RSPB has only 40 reserves with varying levels of local community engagement. For comparison, as of Sep 14 2015, Friends of the Earth Scotland had only 10 local groups (concentrated mostly in large urban areas). Depending on how one measures "volunteerism," it may be possible that ECS has more engaged volunteers in Scotland as well - if each ECS group had only 4 "volunteers" then this would exceed RSPB.
[^15541313]: Kidwell, Jeremy. (2016). Eco-Congregation Scotland, 2014-2016. University of Edinburgh. http://dx.doi.org/10.7488/ds/1357.
[^15541342]:My dataset on transition towns will be made available later in 2016. Initial data was aquired from the Transition Scotland website http://www.transitionscotland.org/transition-in-scotland on December 10, 2014. We are currently in the process of collaboratively generating a more up-to-date dataset which will reflect their collaboration with SCCAN.
[^177171536]: For further detail on Dataset generation, see Kidwell, Forthcoming, 2018.
[^158261232]:Data was acquired from the Development Trusts Association website, http://www.dtascot.org.uk, accessed on 20 July 2015. As above, we are currently in the process of active collaboration with volunteers from the DTAS to co-generate a new dataset.
[^15541614]:PointX data for "Landscape Data" items is sourced from Ordnance Survey Land-Line and MasterMap(R) and the data points are augmented with additional information provided through research by PointX staff, and data aquired from unidentified "local data companie(s)" and the "118 Information" database (see: http://www.118information.co.uk). This data is under license and cannot be made available for use. It is important to note that I became aware of inaccuracies in this dataset over the course of use and subsequently generated my own dataset in collaboration with churches in Scotland. This will be made available later in 2016. I am in active conversation with OS about improving the quality of the data in PointX regarding places of worship.
[^15826124]:Interview with Margaret Warnock, 29 Aug 2014.
[^158261210]:From http://www.forthenvironmentlink.org, accessed 12 July 2015.
[^1554162]:From the Transition map key, "Green pins are 'official' groups
Blue pins are active communities who are connected to the Scottish Transition network Yellow pins show interest in this area"
[^15571030]:This was calculated by calculating a 10m wide footprint for every postcode in Scotland, areas which are not within 10m of a postcode (as of May 2014) are counted as uninhabited.
[^159142242]: Fiona Tweedie, *Ecumenical Audit: Questionnaire Findings* (2014). 
[^15914204]:See note above regarding the data used from the PointX POI database. Note, for our research,we filtered out religious groups not represented within the Eco-Congregation footprint. We discuss representation by tradition and religion further below.adition and religion further below.
