In GRID3/polyCheck: Verifies delineated polygons for the GRID3 project.

# Set knitr options
knitr::opts_chunk$set(echo = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      out.width = "100%",
                      fig.align = "center")

# This chunk will only run if reloadData is set as TRUE. 
# This allows the document to quickly be run as a child document of a parent R Markdown document

# You may need to install bookdown if you haven't used it before
# install.packages("bookdown")
library(raster)
library(sp)
library(sf)
library(leaflet)
library(ggplot2)
library(rgeos)

# --- Load Data

# Settlement Seed
settlements <- params$settlement_path %>%
  shapefile()

# Survey Points
survey_points <- params$survey_points %>%
  shapefile()

survey_points$FID <- 0:(nrow(survey_points)-1) # Add FID to object

# load the ORNL data
ornl <- patams$ornl %>%
  shapefile()

Shapefile Details

shp_path <- params$filepath     # Select the shapefile
shp_name <- basename(shp_path)  # Extract the name
shp <- raster::shapefile(x = shp_path)  # Load the shapefile
shp

Location of Polygon

The map below shows the location of the polygon, relative to the survey point. Note, you may need to zoom out to view the basemaps as ESRI basemaps are not always available at a fully zoomed level.

id <- extract_id(shp_path)   # Extract the ID from the file name
survey_points_id <- survey_points[survey_points$FID == id,]
coords <- coordinates(survey_points_id)

leaflet(shp, width = "100%") %>%
  addProviderTiles(provider = "Esri.WorldImagery") %>%
  addPolygons() %>%
  addMarkers(lng = coords[1], lat = coords[2])

# Convert projection from degrees to metres
crs_m <- "+init=epsg:2062"
shp_m <- spTransform(shp, CRS(crs_m))
survey_points_id_m <- spTransform(survey_points_id, CRS(crs_m))

# Distance in kilometers
dist <- rgeos::gDistance(spgeom1 = shp_m, survey_points_id_m)/1000

cat("**Note**: it appears that the polygon is not in the same location as the polygon. Make sure that the location is not an error")

# Calculate hte intesection
points <- raster::intersect(settlements, shp)
num_points <- nrow(points)

We can verify the number of settlement points contained within the polygon. The region covers r nrow(points). This value should be between 450 and 550.

GRID3/polyCheck documentation built on May 26, 2019, 8:35 a.m.