I am going to use maps to display data for each Unit. I will need their boundaries in order to draw those maps. Will also draw the major highways for geographical references.
library(rgdal,quietly = T) #for reading/writing geo files
## rgdal: version: 1.2-20, (SVN revision 725)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
## Linking to sp version: 1.2-7
library(rgeos,quietly = T) #for simplification and true centroid
## rgeos version: 0.3-26, (SVN revision 560)
## GEOS runtime version: 3.6.1-CAPI-1.10.1 r0
## Linking to sp version: 1.2-7
## Polygon checking: TRUE
library(magrittr,quietly = T)
library(ggplot2,quietly = T) #for the fortify function
library(dplyr,quietly = T)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# https://hub.arcgis.com/items/93c485b0f27e4ff68c066a4b23cbaad3
I’m using a 2015 version, and I’ve saved it in the datasets directory
shapefile <- rgdal::readOGR("~/_code/colorado-dow/datasets/CPW_GMUBoundaries/BigGameGMUBoundaries03172015.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/psarnow/_code/colorado-dow/datasets/CPW_GMUBoundaries/BigGameGMUBoundaries03172015.shp", layer: "BigGameGMUBoundaries03172015"
## with 185 features
## It has 12 fields
## Integer64 fields read as strings: GMUID
shapefile$GMUID <- as.character(shapefile$GMUID)
# convert utm to longlat
shapefile <- spTransform(shapefile, CRS("+proj=longlat +datum=WGS84"))
Will use when drawing maps
Unitboundaries <- shapefile %>% fortify(region = "GMUID")
Unitboundaries2 <- merge(Unitboundaries, shapefile@data, by.x = 'id', by.y = 'GMUID')
Unitboundaries2$Unit <- as.character(Unitboundaries2$id)
Will use for labeling units, and looking up historical Unit weather
# use gCentroid in the rgeos package to calculate the true centroid
shapefile_centroids <- shapefile
shapefile_centroids@data$centroids <- as.data.frame(coordinates(shapefile_centroids))
shapefile_centroids2 <- shapefile_centroids %>% fortify(region = "GMUID")
shapefile_centroids2 <- merge(shapefile_centroids2, shapefile_centroids@data, by.x = 'id', by.y = 'GMUID')
shapefile_centroids2$Unit <- as.character(shapefile_centroids2$id)
shapefile_centroids2$longitude <- shapefile_centroids2$centroids$V1
shapefile_centroids2$latitude <- shapefile_centroids2$centroids$V2
data_centroids <- select(shapefile_centroids2,Unit, group,longitude,latitude)
data_centroids <- data_centroids[!duplicated(data_centroids$Unit),]
data_centroids <- select(data_centroids, longitude, latitude, group, Unit)
data_centroids$Unit <- as.character(data_centroids$Unit)
data_centroids
Will display these on maps for visual references
# Download a statemap with some roads on it (I have saved off into datasets directory)
roaddata <- rgdal::readOGR("~/_code/colorado-dow/datasets/ne_10m_roads/ne_10m_roads.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/psarnow/_code/colorado-dow/datasets/ne_10m_roads/ne_10m_roads.shp", layer: "ne_10m_roads"
## with 56601 features
## It has 29 fields
## Integer64 fields read as strings: scalerank question
USAroads <- roaddata %>% subset(.,sov_a3 == "USA" & type == "Major Highway")
# Convert to data frames so that I can use the data with ggplot2
USAroads <- fortify(USAroads)
# Get min/max of long/lat and zoom to Colorado
longset <- c(min(Unitboundaries2$long),max(Unitboundaries2$long))
latset <- c(min(Unitboundaries2$lat),max(Unitboundaries2$lat))
COroads <- filter(USAroads, long > longset[1] & long < longset[2])
COroads <- filter(COroads, lat > latset[1] & lat < latset[2])
COroads