Setup
library(sf) # Objects and functions for geospatial data
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(rgdal) # Functions for spatial data input/output
## Loading required package: sp
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/proj
library(ggplot2) # Graphing functions
library(dplyr) # Functions for processing tabular data
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr) # Functions for processing tabular data
library(scales) # Additional graphics functions
library(RColorBrewer) # Color ramps for graphs and maps
library(gridExtra) # Functions for arranging multiple plots on a page
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
##
## Attaching package: 'readr'
## The following object is masked from 'package:scales':
##
## col_factor
Get SD Counties
Read the shapefile.
sdcounty <- st_read(dsn = "SD_counties.shp")
## Reading layer `SD_counties' from data source `/Users/haroldnelson/Dropbox/Wimberly/gdswr_chapter4_data/SD_counties.shp' using driver `ESRI Shapefile'
## Simple feature collection with 66 features and 18 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -167412.2 ymin: -264522.5 xmax: 456081.5 ymax: 118174.9
## Projected CRS: Custom
Get the Climate Data
Data is in a csv file.
sdclimate <- read_csv(file = "sd_climate.csv")
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double()
## )
## ℹ Use `spec()` for the full column specifications.
## fips t3001 t3002 t3003
## Min. :46003 Min. :-676.9 Min. :-358.5140 Min. :316.8
## 1st Qu.:46036 1st Qu.:-478.3 1st Qu.:-186.6987 1st Qu.:473.5
## Median :46068 Median :-289.0 Median : 0.0325 Median :609.9
## Mean :46068 Mean :-289.5 Mean : -2.6165 Mean :601.1
## 3rd Qu.:46100 3rd Qu.:-123.4 3rd Qu.: 156.9818 3rd Qu.:718.9
## Max. :46137 Max. : 169.4 Max. : 438.3400 Max. :900.5
## t3004 t3005 t3006 t3007 t3008
## Min. :1144 Min. :1692 Min. :2242 Min. :2719 Min. :2661
## 1st Qu.:1380 1st Qu.:2067 1st Qu.:2564 1st Qu.:2970 1st Qu.:2892
## Median :1501 Median :2135 Median :2667 Median :3068 Median :2983
## Mean :1477 Mean :2130 Mean :2650 Mean :3050 Mean :2971
## 3rd Qu.:1581 3rd Qu.:2217 3rd Qu.:2739 3rd Qu.:3141 3rd Qu.:3064
## Max. :1633 Max. :2272 Max. :2800 Max. :3241 Max. :3175
## t3009 t3010 t3011 t3012 t30sum
## Min. :2089 Min. :1417 Min. :388.9 Min. :-360.8 Min. :1227
## 1st Qu.:2317 1st Qu.:1599 1st Qu.:543.0 1st Qu.:-209.8 1st Qu.:1334
## Median :2431 Median :1698 Median :659.7 Median : -40.8 Median :1462
## Mean :2417 Mean :1696 Mean :651.3 Mean : -41.3 Mean :1443
## 3rd Qu.:2519 3rd Qu.:1797 3rd Qu.:769.4 3rd Qu.: 106.9 3rd Qu.:1540
## Max. :2603 Max. :1871 Max. :869.4 Max. : 312.0 Max. :1646
## p3001 p3002 p3003 p3004
## Min. : 906.8 Min. : 952.8 Min. :1877 Min. :4192
## 1st Qu.:1005.9 1st Qu.:1266.3 1st Qu.:3102 1st Qu.:4942
## Median :1144.7 Median :1380.7 Median :3454 Median :5365
## Mean :1198.6 Mean :1379.6 Mean :3542 Mean :5599
## 3rd Qu.:1279.5 3rd Qu.:1461.7 3rd Qu.:4099 3rd Qu.:6391
## Max. :2227.6 Max. :2291.8 Max. :4920 Max. :7146
## p3005 p3006 p3007 p3008 p3009
## Min. : 6570 Min. : 6921 Min. :5225 Min. :3355 Min. :2788
## 1st Qu.: 7246 1st Qu.: 7820 1st Qu.:6755 1st Qu.:4953 1st Qu.:3589
## Median : 7967 Median : 8465 Median :7381 Median :5673 Median :4757
## Mean : 8028 Mean : 8532 Mean :7339 Mean :5840 Mean :4729
## 3rd Qu.: 8612 3rd Qu.: 9153 3rd Qu.:8123 3rd Qu.:6918 3rd Qu.:5858
## Max. :10019 Max. :10184 Max. :8923 Max. :8077 Max. :6588
## p3010 p3011 p3012 p30sum
## Min. :3252 Min. :1437 Min. : 884.8 Min. :397.8
## 1st Qu.:3878 1st Qu.:1738 1st Qu.:1004.1 1st Qu.:478.9
## Median :4256 Median :2199 Median :1062.9 Median :540.5
## Mean :4256 Mean :2348 Mean :1141.4 Mean :538.8
## 3rd Qu.:4663 3rd Qu.:2898 3rd Qu.:1214.9 3rd Qu.:605.0
## Max. :5118 Max. :3770 Max. :2531.5 Max. :673.8
Join Climate and Geography
sdclimate <- mutate(sdclimate, fips = as.character(fips))
sdcounty <- left_join(sdcounty, sdclimate, by = c("FIPS" = "fips"))
names(sdcounty)
## [1] "STATE" "COUNTY" "FIPS" "C2002" "C2003" "C2004"
## [7] "C2005" "C2006" "C2007" "P2002" "P2003" "P2004"
## [13] "P2005" "P2006" "Perimeter" "Area" "Acres" "Hectares"
## [19] "t3001" "t3002" "t3003" "t3004" "t3005" "t3006"
## [25] "t3007" "t3008" "t3009" "t3010" "t3011" "t3012"
## [31] "t30sum" "p3001" "p3002" "p3003" "p3004" "p3005"
## [37] "p3006" "p3007" "p3008" "p3009" "p3010" "p3011"
## [43] "p3012" "p30sum" "geometry"
Graph
Create a graph with county colored according to C2003. Use ggplot2 and geom_sf()
ggplot(data = sdcounty) +
geom_sf(aes(fill = C2003), color = "white", size = 0.25)
## Use tmap()
The package tmap is an alternative to ggplot2 for maps.
library(tmap)
tm_shape(sdcounty) +
tm_polygons("C2003", palette = "Blues",
style = "cont") +
tm_layout(legend.outside = TRUE,frame = F) +
tm_legend(text.size = 0.7)
## More ggplot2()
ggplot(data = sdcounty) +
geom_sf(aes(fill = C2003), size = 0.25) +
scale_fill_distiller(name="Cases",
palette = "YlGn",
breaks = pretty_breaks()) +
theme_bw() +
labs(title="2003 West Nile Virus Cases in South Dakota")
Graph Centroids
Size is proportional to Population.
ggplot() +
geom_sf(data = sdcounty, aes(fill = C2003), color = "black", size = 0.25) +
scale_fill_distiller(name="Cases", palette = "YlGn", breaks = pretty_breaks())+
geom_sf(data = sdcntrd, aes(size = P2003), color = "black", show.legend = "point") +
scale_size(name = "Population") +
theme_bw() +
labs(title="2003 West Nile Virus Cases in South Dakota")
Alternative with tmap.
tm_shape(sdcounty) +
tm_polygons("C2003",
palette = "YlGn",
style = "cont") +
tm_shape(sdcntrd) +
tm_bubbles(size = "P2003") +
tm_layout(main.title = "2003 West Nile Virus Cases in South Dakota",
legend.outside = TRUE,
frame = F,
legend.outside.size = .5) +
tm_legend(text.size = 0.7)
## Pivot Longer
monthly_prec <- sdcounty %>%
pivot_longer(sdcounty,
cols = p3001:p3012,
values_to = "precip",
names_to = "monthvar") %>%
select(FIPS, monthvar, precip, geometry) %>%
mutate(monthvar = factor(monthvar, labels = month.name))
## Warning in gsub(paste0("^", names_prefix), "", names(cols)): argument 'pattern'
## has length > 1 and only the first element will be used
## Warning in val_cols[col_id] <- unname(as.list(data[cols])): number of items to
## replace is not a multiple of replacement length
Make monthly_prec an sf file.
monthly_prec <- st_as_sf(monthly_prec)
class(monthly_prec)
## [1] "sf" "tbl_df" "tbl" "data.frame"
Graph Monthly Precipitation.
ggplot() +
geom_sf(data = monthly_prec, aes(fill = precip), color = "black", size = 0.25) +
scale_fill_distiller(name="Cases", palette = "Blues", breaks = pretty_breaks())+
facet_wrap(~ monthvar, ncol = 4) +
labs(title="Monthly Climatology of Precipitaton in South Dakota")
