tweaking calculations

mapping_draft.Rmd

@@ -960,21 +960,16 @@ st_crs(permaculture_sf) <- 27700
 
 # Calculate incidence of ecs within SSSI and within buffers at 50/500m
 ecs_sssi_row <- c(sum(apply(st_within(ecs_sf, sssi, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, st_buffer(sssi, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, st_buffer(sssi, dist = 500), sparse=FALSE), 1, any)))
-# Calculate row based on percentage of total
-ecs_sssi_row_pct <- ecs_sssi_row/length(ecs)
 
 pow_sssi_row <- c(sum(apply(st_within(pow_pointX_sf, sssi, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, st_buffer(sssi, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, st_buffer(sssi, dist = 500), sparse=FALSE), 1, any)))
-pow_sssi_row_pct <- pow_sssi_row/length(pow_pointX)
 
 dtas_sssi_row <- c(sum(apply(st_within(dtas_sf, sssi, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, st_buffer(sssi, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, st_buffer(sssi, dist = 500), sparse=FALSE), 1, any)))
-dtas_sssi_row_pct <- dtas_sssi_row/length(dtas)
 
 transition_sssi_row <- c(sum(apply(st_within(transition_sf, sssi, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, st_buffer(sssi, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, st_buffer(sssi, dist = 500), sparse=FALSE), 1, any)))
-transition_sssi_row_pct <- transition_sssi_row/length(transition)
 
 permaculture_sssi_row <- c(sum(apply(st_within(permaculture_sf, sssi, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, st_buffer(sssi, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, st_buffer(sssi, dist = 500), sparse=FALSE), 1, any)))
-permaculture_sssi_row_pct <- permaculture_sssi_row/length(permaculture)
 
+# 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)
@@ -982,12 +977,20 @@ sssi_counts <- rbind(sssi_counts, permaculture_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
+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)
+
 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)
 colnames(sssi_counts_pct) <- c("% Within SSSIs", "% within 50m", "% within 500m")
 
+# Merge into larger dataframe
 sssi_counts_merged <- cbind(sssi_counts, sssi_counts_pct)
 
 
@@ -996,50 +999,70 @@ sssi_counts_merged <- cbind(sssi_counts, sssi_counts_pct)
 titles <- c("Within Wildland Areas", "...50m", "...500m")
 
 ecs_wildland_row <- c(sum(apply(st_within(ecs_sf, wildland, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, st_buffer(wildland, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, st_buffer(wildland, dist = 500), sparse=FALSE), 1, any)))
-
-# TODO: add additional column with counts converted to percentages - use prop.table()?
 ecs_wildland_row <- rbind(titles, ecs_wildland_row)
 
 pow_wildland_row <- c(sum(apply(st_within(pow_pointX_sf, wildland, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, st_buffer(wildland, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, st_buffer(wildland, dist = 500), sparse=FALSE), 1, any)))
-
 wildland_counts <- rbind(ecs_wildland_row, pow_wildland_row)
 
 dtas_wildland_row <- c(sum(apply(st_within(dtas_sf, wildland, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, st_buffer(wildland, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, st_buffer(wildland, dist = 500), sparse=FALSE), 1, any)))
-
 wildland_counts <- rbind(wildland_counts, dtas_wildland_row)
 
 transition_wildland_row <- c(sum(apply(st_within(transition_sf, wildland, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, st_buffer(wildland, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, st_buffer(wildland, dist = 500), sparse=FALSE), 1, any)))
-
 wildland_counts <- rbind(wildland_counts, transition_wildland_row)
 
 permaculture_wildland_row <- c(sum(apply(st_within(permaculture_sf, wildland, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, st_buffer(wildland, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, st_buffer(wildland, dist = 500), sparse=FALSE), 1, any)))
-
 wildland_counts <- rbind(wildland_counts, permaculture_wildland_row)
 
+# Generate dataframe from rows based on percentages of totals
+ecs_wildland_row_pct <- ecs_wildland_row/length(ecs)
+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 forest_inventory buffers
 
 titles <- c("Within woodlands", "...50m", "...500m")
 
 ecs_forest_inventory_row <- c(sum(apply(st_within(ecs_sf, forest_inventory, sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, st_buffer(forest_inventory, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(ecs_sf, st_buffer(forest_inventory, dist = 500), sparse=FALSE), 1, any)))
-
-# TODO: add additional column with counts converted to percentages - use prop.table()?
 ecs_forest_inventory_row <- rbind(titles, ecs_forest_inventory_row)
 
 pow_forest_inventory_row <- c(sum(apply(st_within(pow_pointX_sf, forest_inventory, sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, st_buffer(forest_inventory, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(pow_pointX_sf, st_buffer(forest_inventory, dist = 500), sparse=FALSE), 1, any)))
-
 forest_inventory_counts <- rbind(ecs_forest_inventory_row, pow_forest_inventory_row)
 
 dtas_forest_inventory_row <- c(sum(apply(st_within(dtas_sf, forest_inventory, sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, st_buffer(forest_inventory, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(dtas_sf, st_buffer(forest_inventory, dist = 500), sparse=FALSE), 1, any)))
-
 forest_inventory_counts <- rbind(forest_inventory_counts, dtas_forest_inventory_row)
 
 transition_forest_inventory_row <- c(sum(apply(st_within(transition_sf, forest_inventory, sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, st_buffer(forest_inventory, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(transition_sf, st_buffer(forest_inventory, dist = 500), sparse=FALSE), 1, any)))
-
 forest_inventory_counts <- rbind(forest_inventory_counts, transition_forest_inventory_row)
 
 permaculture_forest_inventory_row <- c(sum(apply(st_within(permaculture_sf, forest_inventory, sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, st_buffer(forest_inventory, dist = 50), sparse=FALSE), 1, any)), sum(apply(st_within(permaculture_sf, st_buffer(forest_inventory, dist = 500), sparse=FALSE), 1, any)))
+forest_inventory_counts <- rbind(forest_inventory_counts, permaculture_forest_inventory_row)
 
-forest_inventory_counts <- rbind(forest_inventory_counts, permaculture_wildland_row)
+# Generate dataframe from rows based on percentages of totals
+ecs_forest_inventory_row_pct <- ecs_forest_inventory_row/length(as.data.frame(ecs))
+pow_forest_inventory_row_pct <- pow_forest_inventory_row/length(pow_pointX)
+dtas_forest_inventory_row_pct <- dtas_forest_inventory_row/length(dtas)
+transition_forest_inventory_row_pct <- transition_forest_inventory_row/length(transition)
+permaculture_forest_inventory_row_pct <- permaculture_forest_inventory_row/length(permaculture)
+
+forest_inventory_counts_pct <- rbind(ecs_forest_inventory_row_pct, pow_forest_inventory_row_pct)
+forest_inventory_counts_pct <- rbind(forest_inventory_counts_pct, dtas_forest_inventory_row_pct)
+forest_inventory_counts_pct <- rbind(forest_inventory_counts_pct, transition_forest_inventory_row_pct)
+forest_inventory_counts_pct <- rbind(forest_inventory_counts_pct, permaculture_forest_inventory_row_pct)
+colnames(forest_inventory_counts_pct) <- c("% Within Woodlands", "% within 50m", "% within 500m")
+
+# Merge into larger dataframe
+forest_inventory_counts_merged <- cbind(forest_inventory_counts, forest_inventory_counts_pct)
 
 ```
 
@@ -1055,18 +1078,23 @@ So what did I discover? The results were inconclusive. First, it is important to
 
 # Output mmd tables using kable
 
-sssi_counts %>%
+sssi_counts_merged %>%
-  kable(format = "html", col.names = colnames(sssi_counts), "Group counts within SSSIs") %>%
+  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 %>%
+wildland_counts_merged %>%
-  kable(format = "html", col.names = colnames(wildland_counts), caption = "Group counts within Wildland Areas") %>%
+  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"))
 
-forest_inventory_counts %>%
+forest_inventory_counts_merged %>%
-  kable(format = "html", col.names = colnames(forest_inventory_counts), caption = "Group counts within woodlands") %>%
+  kable(format = "html", col.names = colnames(forest_inventory_counts_merged), caption = "Group counts within woodlands") %>%
   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(forest_inventory_counts_merged, "derivedData/forest_inventory_counts_merged.csv", row.names=TRUE)
+
 ```
 
 

testing_sf_plot_times.R

@@ -3,12 +3,11 @@ require(sf) # new simplefeature data class, supercedes sp in many ways
 # using GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
 require(sp) # needed for proj4string, deprecated by sf()
 require(rgdal) # version version: 1.3-6
+require(rgeos) # used for buffering below
 require(devtools)
 
-# Set up local workspace and load data
+setwd("~/Downloads/test")
-if (dir.exists("data") == FALSE) {
+# load data
-  dir.create("data") 
-}
 
 transition_wgs <- read.csv(text=getURL("https://zenodo.org/record/165519/files/SCCAN_1.4.csv"))
 coordinates(transition_wgs) <- c("X", "Y")
@@ -16,37 +15,50 @@ proj4string(transition_wgs) <- CRS(wgs84)
 transition_sp <- spTransform(transition_wgs, bng)
 transition_sf <- st_as_sf(transition, coords = c("X", "Y"), crs=27700)
 
+# Download data as ESRI Shapefile from page at: https://gateway.snh.gov.uk/natural-spaces/dataset.jsp?dsid=SSSI
 
-if (file.exists("data/SSSI_SCOTLAND.shp") == FALSE) {
+unzip("SSSI_SCOTLAND_ESRI.zip")
-# TODO: get reliable URL for data download
-# 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")
-}
 
-if (file.exists("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") == FALSE) {
+# # Download data as ESRI Shapefile from page at: http://data-forestry.opendata.arcgis.com/datasets/3cb1abc185a247a48b9d53e4c4a8be87_0
-download.file("https://opendata.arcgis.com/datasets/3cb1abc185a247a48b9d53e4c4a8be87_0.zip", destfile = "data/National_Forest_Inventory_Woodland_Scotland_2017.zip")
 
-unzip("data/National_Forest_Inventory_Woodland_Scotland_2017.zip", exdir = "data")
+unzip("National_Forest_Inventory_Woodland_Scotland_2017.zip")
-}
 
-forest_inventory_sf <- st_read("data/National_Forest_Inventory_Woodland_Scotland_2017.shp")
+sssi_sf <- st_read("SSSI_SCOTLAND.shp")
-forest_inventory_sp <- readOGR("./data", "National_Forest_Inventory_Woodland_Scotland_2017") 
+sssi_sp <- readOGR("./", "SSSI_SCOTLAND") 
 
-st_crs(sssi) <- 27700
+forest_inventory_sf <- st_read("National_Forest_Inventory_Woodland_Scotland_2017.shp")
-st_crs(ecs_sf) <- 27700
+forest_inventory_sp <- readOGR("./", "National_Forest_Inventory_Woodland_Scotland_2017") 
+
+# First test out plots using spatialfeatures and spdf with core R
 
 system.time(
-ggplot() +
+plot(sssi_sf)
-  geom_sf(data = forest_inventory_sf)
   )
 
 system.time(
-ggplot() +
+plot(sssi_sp)
-  geom_polygon(data = forest_inventory_sp)
   )
 
-count_data <- sum(apply(st_within(points_sf, st_buffer(sssi, dist = 50), sparse=FALSE), 1, any))
+# First test out plots using spatialfeatures and spdf with ggplot2
 
+system.time(
+  ggplot() +
+    geom_sf(data = sssi_sf)
+)
+
+system.time(
+  ggplot() +
+    geom_polygon(data = sssi_sp)
+)
+# Now try to run a count within a buffer:
+
+st_crs(sssi_sf) <- 27700
+st_crs(transition_sf) <- 27700
+
+# CRS uses meters for units, so buffer here should be a modest 50m:
+
+count_data_sf <- sum(apply(st_within(points_sf, st_buffer(sssi, dist = 50), sparse=FALSE), 1, any))
+
+# count_data_sf <- sum(apply(gWithin(points_sf, gBuffer(sssi,width=50)
   
 sessioninfo::session_info()