final fixes to conversion over to sf() and tidyverse

mapping_draft-hpc_optimised.Rmd

@@ -386,7 +386,7 @@ tm_shape(admin_lev2) +
 
 
 
-```{r 02_admin_ecs_normed_choropleth, fig.width=4, fig.show="hold", fig.cap="Figure 3"}
+```{r 02_admin_ecs_normed_choropleth, fig.width=4, fig.show="hold", fig.cap="Figure 2"}
 
 # Plot out first figure with normalised data: 
 
@@ -403,7 +403,7 @@ tm_shape(admin_lev2) +
               just = c("left", "bottom"),
               align = "left") +
   tm_layout(asp = NA,
-            title = "Figure 3a",
+            title = "Figure 2a",
             frame = FALSE, 
             title.size = .7, 
             legend.title.size = .7,
@@ -414,7 +414,7 @@ tm_shape(admin_lev2) +
 
 tm_shape(admin_lev2) + 
   tm_fill(col = "ecs_count_popnorm_scaled", palette = "Oranges", n = 5, 
-          title = "Concentration of ECS groups, data normalised by places of worship") +
+          title = "Concentration of ECS groups, data normalised by population density") +
   tm_borders(alpha=.5, lwd=0.1) +
   tm_shape(admin_lev1_simplified) +
   tm_borders(lwd=0.6) +
@@ -427,14 +427,14 @@ tm_shape(admin_lev2) +
               align = "left") +
   tm_layout(asp = NA,
             frame = FALSE, 
-            title = "Figure 3b",
+            title = "Figure 2b",
             title.size = .7, 
             legend.title.size = .7,
             inner.margins = c(0.1, 0.1, 0.05, 0.05)
             )
 ```
 
-Given the way population and places of worship are unevenly distributed across Scotland it is important to represent data in terms of relative distribution. For this study, we attempted to "normalise" our data in two different ways, (1) as shown by Figure 2 above, by taking population figures from the 2011 census (see data sheet in Appendix A) and (2) by adjusting relative to the number of places of worship in each council region.[^15914204] The latter of these two can yield particularly unexpected results. Thus, of the `r length(pow_pointX)` "places of worship" in Scotland, the highest concentration is actually the `r as.character(admin_lev1$NAME_2[which.max(admin_lev1$pow_count)])` region, with `r max(admin_lev1$pow_count)`, second is `r max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)] )` (`r as.character(admin_lev1$NAME_2[which.max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)])] )`). Rank of Council Areas by population and number of places of worship is also included in Appendix A.
+Given the way population and places of worship are unevenly distributed across Scotland it is important to represent data in terms of relative distribution. For this study, we attempted to "normalise" our data in two different ways, (1) as shown in Figure 2a above, by adjusting relative to the number of places of worship in each council region  and (2) as shown by Figure 2b above, by taking population figures from the 2011 census (see data sheet in Appendix A).[^15914204] The latter of these two can yield particularly unexpected results. Thus, of the `r length(pow_pointX)` "places of worship" in Scotland, the highest concentration is actually the `r as.character(admin_lev1$NAME_2[which.max(admin_lev1$pow_count)])` region, with `r max(admin_lev1$pow_count)`, second is `r max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)] )` (`r as.character(admin_lev1$NAME_2[which.max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)])] )`). Rank of Council Areas by population and number of places of worship is also included in Appendix A.
 
 ```{r create_admin_proportions}
 # Calculate factors by which ECS representation exceeds rep by population and total pow counts
@@ -450,7 +450,7 @@ Turning to the total of `r length(pow_pointX)` "places of worship" in Scotland,
 
 Whereas our initial measurements indicated a prominent lead for Edinburgh, by normalising our data in this way we can highlight the stronger-than-expected presence of several others that might otherwise escape notice because they lie in a region with significantly lower population or numerically less places of worship. Taking the PointX data on "places of worship" in Scotland, we find a less dramatic picture, but also a slightly different one. The positive outliers include East Renfrewshire (3.4x) Edinburgh (2.9x), Stirling (2.2), West Lothian (1.9x) and Aberdeen (1.5x). Again, negative outliers are far less dramatic, with only Midlothian possessing a ratio of more than 100% negative difference from the number of "places of worship" at 1.5x *fewer*.
 
-```{r 03_admin_barplot, fig.width=4, fig.cap="Figure 4"}
+```{r 03_admin_barplot, fig.width=4, fig.cap="Figure 3"}
 # comvert admin back to dataframe for analysis
 admin.df <- data.frame(admin_lev1)
 
@@ -465,8 +465,8 @@ admin_barplot <-
     geom_bar(position="dodge", stat="identity") + 
     coord_flip() + 
     labs(
-     title = "Figure 4", 
+     title = "Figure 3", 
-     subtitle ="Comparison of Groups by Admin Region", 
+     subtitle ="Comparison of Groups by Administrative Region", 
      fill = "Groups") +
     theme(
      legend.position="bottom", 
@@ -476,13 +476,13 @@ admin_barplot <-
 # TODO, need better export method now that tmap() and graphic device are fixed
 pdf("figures/03_admin_barplot.pdf")
 print(admin_barplot)
-# dev.off()
+dev.off()
 
 ```
 
 ## Other environmental groups shown by concentration in administrative regions (NUTS)
 
-```{r 04_choropleth_others, fig.width=4, fig.show="hold", fig.cap="Figure 5"}
+```{r 04_choropleth_others, fig.width=4, fig.show="hold", fig.cap="Figure 4"}
 
 # TODO: consider switching to two-dimensional kernel densities instead of dots as shown here: https://github.com/mtennekes/tmap/tree/master/demo/LondonCrimes
 
@@ -503,7 +503,7 @@ admin_lev2_scotland_ecs_plot <-
               align = "left") +
   tm_layout(asp = NA,
             frame = FALSE, 
-            title = "Figure 5a", 
+            title = "Figure 4a", 
             title.size = .7,
             legend.title.size = .7,
             inner.margins = c(0.1, 0.1, 0.05, 0.05),
@@ -530,7 +530,7 @@ pdf("figures/03_admin_lev2_scotland_ecs_plot.pdf")
 print(admin_lev2_centralbelt_ecs_plot, 
       vp = vp_admin_lev2_centralbelt_ecs_plot)
 # TODO, need better export method now that tmap() and graphic device are fixed
-# dev.off()
+dev.off()
 
 # Second plot, revealing transition towns
 
@@ -550,7 +550,7 @@ tm_shape(admin_lev2) +
               align = "left") +
   tm_layout(asp = NA,
             frame = FALSE, 
-            title = "Figure 5b", 
+            title = "Figure 4b", 
             title.size = .7,
             legend.title.size = .7,
             inner.margins = c(0.1, 0.1, 0.05, 0.05),
@@ -573,7 +573,7 @@ tm_shape(admin_lev2) +
               align = "left") +
   tm_layout(asp = NA,
             frame = FALSE, 
-            title = "Figure 5c", 
+            title = "Figure 4c", 
             title.size = .7,
             legend.title.size = .7,
             inner.margins = c(0.1, 0.1, 0.05, 0.05),
@@ -597,7 +597,7 @@ tm_shape(admin_lev2) +
               align = "left") +
   tm_layout(asp = NA,
             frame = FALSE, 
-            title = "Figure 5d", 
+            title = "Figure 4d", 
             title.size = .7,
             legend.title.size = .7,
             inner.margins = c(0.1, 0.1, 0.05, 0.05),
@@ -607,7 +607,7 @@ tm_shape(admin_lev2) +
 
 ## Cartogram Comparisons
 
-```{r 05_cartograms, fig.width=4, fig.show="hold", fig.cap="Figure 6"}
+```{r 05_cartograms, fig.width=4, fig.show="hold", fig.cap="Figure 5"}
 
 # # TODO: plot as animated chorogram: 
 # # https://www.r-graph-gallery.com/331-basic-cartogram/
@@ -636,7 +636,7 @@ So why provide this kind of data (i.e. at the level of individual churches) when
 ```{r 06_ecs_denomination_table}
 # TODO: fix to work with new sf() object
 ecs_dt <- (st_set_geometry(ecs, NULL))
-knitr::kable(table(ecs_dt$denomination), caption = 'ECS by denomination')
+knitr::kable(table(ecs_dt$denomination), caption = 'Figure 5. Table of ECS by denomination')
 # TODO: Add dataframe with overall church numbers for the UK from "ScottishChurches" dataset by JK
 ```
 
@@ -694,7 +694,7 @@ The key question which this analysis seeks to answer is whether ECS, or the othe
 
 Of all the groups surveyed in this study, Eco-Congregation Scotland is the most heavily concentrated in large urban areas (33.53%), exceeding by almost 50% the rate for all places of worship (22.96% in large urban areas). Transition is a much more modest 20% and development trusts a bit lower at 15%. It is interesting to note that the rate of ECS concentration in these large urban areas matches the level of overall population distribution (34.5%). On the other end of the scale, Eco-Congregation Scotland is the least concentrated in remote rural areas (with 3.93% on level 7 and 5.44% on level 8 on the urban-rural scale), though again, they correlate roughly to the general population distribution (3.2% and 2.9% respectively). Places of worship outpace both the population of Scotland and the footprint of Eco-Congregation Scotland, with 14.98% in very remote rural areas, but this is exceeded by transition at 16.47% and both by Scottish community development trusts at 32.14%. So while Eco-Congregation Scotland correlates roughly with Scottish population distribution across the urban-rural scale, it has a considerably more urban profile than either of the other two groups surveyed.
 
-```{r 07_ur_barplot, fig.width=4, fig.cap="Figure 7"}
+```{r 07_ur_barplot, fig.width=4, fig.cap="Figure 6"}
 # Need to flatten urbanrural based on all the count columns and generate using ggplot
 
 urbanrural_gathered <- gather(data.frame(urbanrural), key="group_type", value="number", ecs_count, transition_count, dtas_count, permaculture_count)
@@ -706,13 +706,13 @@ ggplot(urbanrural_gathered,
   geom_bar(position="dodge", stat="identity") + 
   coord_flip() + 
   labs(
-    title = "Figure 7", 
+    title = "Figure 6", 
     subtitle="Comparison of Groups by UrbanRural category", 
     fill = "Groups")
 ```
 ## Eco-Congregation Scotland\nconcentrations in Urban Rural 8-fold classifications
 
-```{r 08_urbanrural_ecs_chart_choropleth, message=FALSE, warning=FALSE, fig.width=4, fig.cap="Figure 8"}
+```{r 08_urbanrural_ecs_chart_choropleth, message=FALSE, warning=FALSE, fig.width=4, fig.cap="Figure 7"}
 
 # TODO: Double check data licenses for tm_credits
 # Generate code for inset map of central belt
@@ -732,7 +732,7 @@ urbanrural_uk_ecs_choropleth_plot <- tm_shape(urbanrural_simplified) +
               just = c("left", "bottom"),
               align = "left") +
   tm_layout(asp = NA,
-            title = "Figure 8", 
+            title = "Figure 7", 
             frame = FALSE, 
             title.size = .7, 
             legend.title.size = .7,
@@ -794,8 +794,8 @@ simd_shapes <- st_read("data/sc_dz_11.shp")
 simd_indicators <- read.csv("./data/simd2016_withinds.csv")
 simd_indicators_min <- simd_indicators[c(1,6:17)]
 # Combine shapes with data
+simd_wgs = left_join(simd_shapes, simd_indicators_min, by = c("DataZone" = "Data_Zone"))
 # Original data is in wgs, so need to convert from wgs to bng to work with analysis across data sets
-simd_wgs <- merge(x=simd_shapes, y=simd_indicators_min, by.x = "DataZone", by.y = "Data_Zone")
 simd <- simd_wgs %>% st_transform(27700)
  # Create simplified version of geometry for visualisations, with a tolerance of 100m
 # Workaround to address invalid spherical geometry issue 
@@ -806,19 +806,21 @@ round(c(object.size(simd_shapes), object.size(simd_simplified)) / 1024)
 
 # Augment each dataset with relevant (geolocated) columns from SIMD
 # assign combined table with SIMD columns to attribute table slot of ecs table
+ecs <- st_join(ecs, simd )
 
-# TODO: need to resolve issue here with simd columns not being added, also issues with geometry
-ecs=cbind(ecs,st_intersects(simd,ecs))
 # assign combined table with SIMD columns to attribute table slot of ecs table
-pow_pointX=cbind(pow_pointX,st_intersects(pow_pointX,simd))
+pow_pointX <- st_join(pow_pointX, simd )
+
 # assign combined table with SIMD columns to attribute table slot of transition table
 # TODO: fix bad geometry here causing fails as it lies outside simd polygons
-# transition=cbind(transition,st_intersects(transition,simd))
+transition <- st_join(transition, simd )
+
 # assign combined table with SIMD columns to attribute table slot of permaculture table
-permaculture=cbind(permaculture,st_intersects(permaculture,simd))
+permaculture <- st_join(permaculture, simd )
+
 # assign combined table with SIMD columns to attribute table slot of dtas table
 # TODO: fix bad geometry here causing fails as it lies outside simd polygons
-# dtas=cbind(dtas,st_intersects(dtas,simd))
+dtas <- st_join(dtas, simd )
 
 # Augment SIMD with group counts
 
@@ -829,19 +831,19 @@ simd$permaculture_count <- lengths(st_covers(simd, permaculture))
 simd$pow_count <- lengths(st_covers(simd, pow_pointX))
 
 # Generate simplified dataframes with counts for each group
-# ecs_simd <- data.frame(ecs[c(1,16,25,30:36)])
+ecs_simd <- data.frame(select(ecs, name, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank))
-# TODO: adapt code below to include equivalent rows from above now that we are using tidyverse
-ecs_simd <- data.frame(select(ecs, row.id, col.id))
-
 colnames(ecs_simd)[1] <- "group_name"
 ecs_simd$group_type <- "ecs"
-transition_simd <- data.frame(transition[c(2,5,14,19:25)])
+
+transition_simd <- data.frame(select(transition, Name.x, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank))
 colnames(transition_simd)[1] <- "group_name"
 transition_simd$group_type <- "transition"
-dtas_simd <- data.frame(dtas[c(1,5,14,19:25)])
+
+dtas_simd <- data.frame(select(dtas, Name.x, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank))
 colnames(dtas_simd)[1] <- "group_name"
 dtas_simd$group_type <- "dtas"
-permaculture_simd <- data.frame(permaculture[c(1,11,20,25:31)])
+
+permaculture_simd <- data.frame(select(permaculture, Link, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank))
 colnames(permaculture_simd)[1] <- "group_name"
 permaculture_simd$group_type <- "permaculture"
 
@@ -865,7 +867,7 @@ allgroups_gathered <- gather(allgroups_simd, key = "simd_category", value = "ran
 
 ### Barplot
 
-```{r 11_simd_barplot, fig.width=4, fig.cap="Figure 9"}
+```{r 11_simd_barplot, fig.width=4, fig.cap="Figure 8"}
 # Run plots
 
 # Faceted stacked bar plot
@@ -889,14 +891,13 @@ ggplot(data=allgroups_gathered,
   labs(x = NULL, 
        y = "SIMD Rank (in bins by quantile)", 
        fill = "Groups",
-       title = "Figure 9", 
+       title = "Figure 8", 
        subtitle="Distribution of Groups across IMD Domains by Rank",
        caption = paste("Jeremy H. Kidwell :: jeremykidwell.info",
                        "Data: UK Data Service (OGL) & Jeremy H. Kidwell",
                        "You may redistribute this graphic under the terms of the CC-by-SA 4.0 license.",
                        sep = "\n"))
   
-
 # write.csv(simd_percents_only, "derivedData/simd_percents_only.csv", row.names=FALSE)
 ```
 
@@ -934,7 +935,6 @@ if (file.exists("data/SSSI_SCOTLAND.shp") == FALSE) {
 unzip("data/SSSI_SCOTLAND_ESRI.zip", exdir = "data")
 }
 sssi <- st_read("data/SSSI_SCOTLAND.shp")
-sssi_sp <- readOGR("./data", "SSSI_SCOTLAND") 
  # Create simplified version of geometry for visualisations, with a tolerance of 100m
 sssi_simplified <- st_simplify(sssi, preserveTopology = TRUE, dTolerance = 100)
 # compare memory usage of objects
@@ -1069,7 +1069,54 @@ ecs_scenicareasbeyond500m <- !(st_within(ecs, scenicareas_buf500))
 # See: https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html
 
 # Generate dataframe based on SSSI buffers
-# TODO: add pubs and grocery store data (using PointX)
+
+# Load in POI tables so that we can calculate coincidence of other features
+
+# Sadly, this dataset is paywalled by digimap, though most academic researchers should be able to aquire it through 
+# an institutional subscription from Digimap, Ordnance Survey Data Download. The POI set is located under 
+# "Boundary and Location Data" under "Points of Interest," "select visible area" for whole UK then download.
+#
+# I'm going to take some shortcuts here, preventing this code from being truly reproducible
+# because the poi dataset is really huge - I've already downloaded and filtered to specific codes for features
+# so code immediately below is untested. Will need to adjust below to the actual code I'm working with eventually.
+
+## Skip ==>
+# poi <- read.csv("data/poi_csv", header = FALSE, sep = "|")
+# coordinates(poi) <- c("Feature Easting", "Feature Northing")
+# proj4string(poi) <- proj4string(CRS(bng))
+## transform CRS on churches (from BNG) to wgs84 for comparisons below
+# poi_wgs <- spTransform(poi, CRS(wgs84))
+
+## subsets
+#poi_pubs <- poi[!is.na(poi[3] = "01020034")]
+#poi_chequecashing <- poi[!is.na(poi[3] = "02090142")]
+#poi_pawnbrokers <- poi[!is.na(poi[3] = "02090151")]
+#poi_worship <- poi[!is.na(poi[3] = "06340459")]
+# ==> to here for now
+
+# 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 <- select(poi_grocery_wgs, retailer, store_name, long_wgs, lat_wgs)
+# convert to sf
+poi_grocery_wgs <- st_as_sf(poi_grocery_wgs, coords = c("long_wgs", "lat_wgs"))
+st_crs(poi_grocery_wgs) <- 4326
+poi_grocery <- poi_grocery_wgs %>% st_transform(27700)
+# filter out non-Scottish data
+ poi_grocery_scotland <- st_intersection(poi_grocery, admin_lev1)
+
+# Load in British pubs from Ordnance survey dataset
+poi_pubs <- read.csv("data/poi_pubs.csv", header = FALSE, sep = "|")
+# select useful columns
+poi_pubs <- select(poi_pubs, V1, V2, V3, V4, V5)
+# rename columns to tidier names
+colnames(poi_pubs) <- c("refnum", "name", "code", "x", "y")
+# convert to sf
+poi_pubs <- st_as_sf(poi_pubs, coords = c("x", "y"))
+st_crs(poi_pubs) <- 27700
+# filter out non-Scottish data
+poi_pubs <- st_intersection(poi_pubs, admin_lev1)
 
 # Calculate incidence of ecs within SSSI and within buffers at 50/500m
 ecs_sssi_row <- c(sum(apply(st_within(ecs, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(ecs, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(ecs, sssi_buf500, sparse=FALSE), 1, any)))
@@ -1188,6 +1235,7 @@ colnames(forestinv_counts) <- c("Within Woodlands", "...50m", "...500m")
 
 # Generate dataframe from rows based on percentages of totals
 # TODO: fix error generated by ecs_forestinv_row_pct using ecs as sf().
+# Shift to using table summary function
 ecs_forestinv_row_pct <- ecs_forestinv_row/length(ecs)
 pow_forestinv_row_pct <- pow_forestinv_row/length(pow_pointX)
 dtas_forestinv_row_pct <- dtas_forestinv_row/length(dtas)
@@ -1369,7 +1417,8 @@ 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 <- st_set_geometry(admin_lev1)
+admin_lev1_df <- admin_lev1_df[,c(3,5,7,11,13)]
 
 admin_lev1_df %>%
   kable(format = "html", col.names = colnames(admin_lev1_df)) %>%
@@ -1381,8 +1430,7 @@ admin_lev1_df %>%
 
 ```{r admin_table_withproportions}
 
-admin_lev1_prop_df <- as_data_frame(admin_lev1[,c(3,5:14,22:25)])
+admin_lev1_prop_df <- admin_lev1_df[,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")) %>%
@@ -1396,8 +1444,8 @@ admin_lev1_prop_df %>%
 write.csv(urbanrural, "derivedData/urbanrural.csv", row.names=FALSE)
 
 # Output mmd tables using kable
-
+urbanrural_table <- st_set_geometry(urbanrural)
-urbanrural_table <- as_data_frame(urbanrural[,c(1,7:16)])
+urbanrural_table <- urbanrural_table[,c(1,7:16)]
 
 urbanrural_table %>%
   kable(format = "html", col.names = colnames(urbanrural_table)) %>%

sbatch_wild.sh

@@ -1,18 +0,0 @@
-#!/bin/bash
-#SBATCH --mail-type ALL
-#SBATCH --cpus-per-task 1
-#SBATCH --time 180:0
-#SBATCH --ntasks 20
-#SBATCH --qos bbdefault
-
-module purge; module load bluebear
-module load R/3.5.0-iomkl-2018a-X11-20180131
-module load Pandoc/2.3.1
-module load plotly/4.7.1-iomkl-2018a-R-3.5.0
-module load rgeos/0.3-28-iomkl-2018a-R-3.5.0
-module load tmap/2.2-iomkl-2018a-R-3.5.0
-
-Rscript -e 'library(rmarkdown); rmarkdown::render("mapping_draft-hpc_optimised_wilderness.Rmd", "html_document")'
-# cp mapping_draft.html /rds/projects/2016/kidwellj-01/mapping_environmental_action/
-# cp figures/* /rds/projects/2016/kidwellj-01/mapping_environmental_action/figures
-# cp derivedData/* /rds/projects/2016/kidwellj-01/mapping_environmental_action/derivedData