diff --git a/mapping_draft.Rmd b/mapping_draft.Rmd index ba9d7c3..9e4b12c 100644 --- a/mapping_draft.Rmd +++ b/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 %>% - kable(format = "html", col.names = colnames(sssi_counts), "Group counts within SSSIs") %>% +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 %>% - kable(format = "html", col.names = colnames(wildland_counts), caption = "Group counts within Wildland Areas") %>% +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")) -forest_inventory_counts %>% - kable(format = "html", col.names = colnames(forest_inventory_counts), caption = "Group counts within woodlands") %>% +forest_inventory_counts_merged %>% + 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) + ``` diff --git a/testing_sf_plot_times.R b/testing_sf_plot_times.R index b0541c3..1086e85 100644 --- a/testing_sf_plot_times.R +++ b/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 -if (dir.exists("data") == FALSE) { - dir.create("data") -} +setwd("~/Downloads/test") +# load 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) { -# 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") -} +unzip("SSSI_SCOTLAND_ESRI.zip") -if (file.exists("data/National_Forest_Inventory_Woodland_Scotland_2017.shp") == FALSE) { -download.file("https://opendata.arcgis.com/datasets/3cb1abc185a247a48b9d53e4c4a8be87_0.zip", destfile = "data/National_Forest_Inventory_Woodland_Scotland_2017.zip") +# # Download data as ESRI Shapefile from page at: http://data-forestry.opendata.arcgis.com/datasets/3cb1abc185a247a48b9d53e4c4a8be87_0 -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") -forest_inventory_sp <- readOGR("./data", "National_Forest_Inventory_Woodland_Scotland_2017") +sssi_sf <- st_read("SSSI_SCOTLAND.shp") +sssi_sp <- readOGR("./", "SSSI_SCOTLAND") -st_crs(sssi) <- 27700 -st_crs(ecs_sf) <- 27700 +forest_inventory_sf <- st_read("National_Forest_Inventory_Woodland_Scotland_2017.shp") +forest_inventory_sp <- readOGR("./", "National_Forest_Inventory_Woodland_Scotland_2017") + +# First test out plots using spatialfeatures and spdf with core R system.time( -ggplot() + - geom_sf(data = forest_inventory_sf) +plot(sssi_sf) ) system.time( -ggplot() + - geom_polygon(data = forest_inventory_sp) +plot(sssi_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() \ No newline at end of file