Chapter 5 Notes Part 1
Harold Nelson
2023-02-13
Setup
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
## Loading required package: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-7
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(ggplot2)
library(dplyr)
##
## 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)
library(scales)
library(RColorBrewer)
library(units)
## udunits database from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/units/share/udunits/udunits2.xml
Import the Data
okcounty <- st_read("ok_counties.shp", quiet = TRUE)
tpoint <- st_read("ok_tornado_point.shp", quiet = TRUE)
tpath <- st_read("ok_tornado_path.shp", quiet = TRUE)
class(okcounty)
## [1] "sf" "data.frame"
## [1] "sf" "data.frame"
glimpse(okcounty)
## Rows: 77
## Columns: 8
## $ STATEFP <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "…
## $ COUNTYFP <chr> "077", "025", "011", "107", "105", "153", "001", "053", "059"…
## $ COUNTYNS <chr> "01101826", "01101800", "01101793", "01101841", "01101840", "…
## $ AFFGEOID <chr> "0500000US40077", "0500000US40025", "0500000US40011", "050000…
## $ GEOID <chr> "40077", "40025", "40011", "40107", "40105", "40153", "40001"…
## $ NAME <chr> "Latimer", "Cimarron", "Blaine", "Okfuskee", "Nowata", "Woodw…
## $ LSAD <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "…
## $ geometry <POLYGON [°]> POLYGON ((-95.50766 35.0292..., POLYGON ((-103.0025 3…
Some WA Data
Go to https://geo.wa.gov/datasets/wadnr::wa-county-boundaries/explore?location=47.245765%2C-120.817600%2C8.12
Download the shapefile and import it. Move the files to your Chapter
Five Project folder first. Then use geom_sf() to look at it.
Solution
wa_counties_shp = st_read("WA_County_Boundaries.shp")
## Reading layer `WA_County_Boundaries' from data source
## `/Users/haroldnelson/Library/CloudStorage/Dropbox/GDSWR/gdswr_data/Chapter5/WA_County_Boundaries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 39 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -13899440 ymin: 5707530 xmax: -13014940 ymax: 6275274
## Projected CRS: WGS 84 / Pseudo-Mercator
## Rows: 39
## Columns: 12
## $ OBJECTID <int> 2457, 2461, 6627, 6989, 7610, 7660, 7679, 7683, 7715, 7724,…
## $ JURISDICT_ <int> 25, 37, 33, 8, 4698173, 35, 34, 3, 20, 32, 39, 29, 4698174,…
## $ JURISDIC_1 <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ JURISDIC_2 <chr> "Grant", "Benton", "Garfield", "Island", "Stevens", "Walla …
## $ JURISDIC_3 <chr> "Grant County", "Benton County", "Garfield County", "Island…
## $ JURISDIC_4 <int> 13, 3, 12, 15, 33, 36, 7, 28, 14, 2, 30, 11, 26, 37, 5, 25,…
## $ JURISDIC_5 <int> 53025, 53005, 53023, 53029, 53065, 53071, 53013, 53055, 530…
## $ JURISDIC_6 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EDIT_DATE <date> 2018-03-15, 2022-07-13, 2022-06-23, 2018-03-15, 2022-05-18…
## $ EDIT_STATU <int> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1,…
## $ EDIT_WHO <chr> "TSTE490", "TSTE490", "TSTE490", "TSTE490", "TSTE490", "JDU…
## $ geometry <POLYGON [m]> POLYGON ((-13245041 6100462..., POLYGON ((-13292341…
wa_counties_shp %>%
ggplot() +
geom_sf()
The CSV File
Try the same thing with the csv file.
Solution
wa_counties_csv = st_read("WA_County_Boundaries.csv")
## Reading layer `WA_County_Boundaries' from data source
## `/Users/haroldnelson/Library/CloudStorage/Dropbox/GDSWR/gdswr_data/Chapter5/WA_County_Boundaries.csv'
## using driver `CSV'
## Warning: no simple feature geometries present: returning a data.frame or tbl_df
## Rows: 39
## Columns: 11
## $ OBJECTID <chr> "2457", "2461", "6627", "6989", "7610", "7660",…
## $ JURISDICT_SYST_ID <chr> "25", "37", "33", "8", "4698173", "35", "34", "…
## $ JURISDICT_TYPE_CD <chr> "4", "4", "4", "4", "4", "4", "4", "4", "4", "4…
## $ JURISDICT_LABEL_NM <chr> "Grant", "Benton", "Garfield", "Island", "Steve…
## $ JURISDICT_NM <chr> "Grant County", "Benton County", "Garfield Coun…
## $ JURISDICT_DESG_CD <chr> "13", "3", "12", "15", "33", "36", "7", "28", "…
## $ JURISDICT_FIPS_DESG_CD <chr> "53025", "53005", "53023", "53029", "53065", "5…
## $ JURISDICT_VACATED_FLG <chr> "", "", "", "", "", "", "", "", "", "", "", "",…
## $ EDIT_DATE <chr> "2018/03/15 00:26:49+00", "2022/07/13 23:45:43+…
## $ EDIT_STATUS <chr> "1", "1", "1", "1", "0", "1", "1", "1", "1", "1…
## $ EDIT_WHO <chr> "TSTE490", "TSTE490", "TSTE490", "TSTE490", "TS…
wa_counties_csv %>%
ggplot() +
geom_sf()
## Error in `check_required_aesthetics()`:
## ! stat_sf requires the following missing aesthetics: geometry
Solution
wa_counties_kml = st_read("WA_County_Boundaries.kml")
## Reading layer `WA_County_Boundaries' from data source
## `/Users/haroldnelson/Library/CloudStorage/Dropbox/GDSWR/gdswr_data/Chapter5/WA_County_Boundaries.kml'
## using driver `KML'
## Simple feature collection with 39 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -124.8608 ymin: 45.54373 xmax: -116.9152 ymax: 49.00243
## Geodetic CRS: WGS 84
## Rows: 39
## Columns: 3
## $ Name <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""…
## $ Description <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""…
## $ geometry <POLYGON [°]> POLYGON ((-118.9822 47.9615..., POLYGON ((-119.407…
wa_counties_kml %>%
ggplot() +
geom_sf()
Solution
wa_counties_geojson = st_read("WA_County_Boundaries.geojson")
## Reading layer `WA_County_Boundaries' from data source
## `/Users/haroldnelson/Library/CloudStorage/Dropbox/GDSWR/gdswr_data/Chapter5/WA_County_Boundaries.geojson'
## using driver `GeoJSON'
## Simple feature collection with 39 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -124.8608 ymin: 45.54373 xmax: -116.9152 ymax: 49.00243
## Geodetic CRS: WGS 84
glimpse(wa_counties_geojson)
## Rows: 39
## Columns: 12
## $ OBJECTID <int> 2457, 2461, 6627, 6989, 7610, 7660, 7679, 7683,…
## $ JURISDICT_SYST_ID <int> 25, 37, 33, 8, 4698173, 35, 34, 3, 20, 32, 39, …
## $ JURISDICT_TYPE_CD <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ JURISDICT_LABEL_NM <chr> "Grant", "Benton", "Garfield", "Island", "Steve…
## $ JURISDICT_NM <chr> "Grant County", "Benton County", "Garfield Coun…
## $ JURISDICT_DESG_CD <int> 13, 3, 12, 15, 33, 36, 7, 28, 14, 2, 30, 11, 26…
## $ JURISDICT_FIPS_DESG_CD <int> 53025, 53005, 53023, 53029, 53065, 53071, 53013…
## $ JURISDICT_VACATED_FLG <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ EDIT_DATE <dttm> 2018-03-14 17:26:49, 2022-07-13 16:45:43, 2022…
## $ EDIT_STATUS <int> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1,…
## $ EDIT_WHO <chr> "TSTE490", "TSTE490", "TSTE490", "TSTE490", "TS…
## $ geometry <POLYGON [°]> POLYGON ((-118.9822 47.9615..., POLYGON…
wa_counties_geojson %>%
ggplot() +
geom_sf()
Dissecting a Graphic
The following chunk runs a graphic from the text. First, just run it
as written.
Note that the years are four-digit, not 2.
tpoint_16_21 = tpoint %>%
filter(yr >= 2016)
ggplot() +
geom_sf(data = tpoint_16_21,
aes(color = as.factor(yr))) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()
Why factor(yr)
Remove this and see what happens.
Note that scale_color_discrete becomes invalid. You need
scale_color_continuous instead.
tpoint_16_21 = tpoint %>%
filter(yr >= 2016)
ggplot() +
geom_sf(data = tpoint_16_21,
aes(color = yr)) +
geom_sf(data = okcounty, fill = NA) +
scale_color_continuous(name = "Year") +
coord_sf(datum = NA) +
theme_void()
coord_sf()?
What happens if you remove datum = NA?
Solution
tpoint_16_21 = tpoint %>%
filter(yr >= 2016)
ggplot() +
geom_sf(data = tpoint_16_21,
aes(color = as.factor(yr))) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf() +
theme_void()
## Remove Completely
tpoint_16_21 = tpoint %>%
filter(yr >= 2016)
ggplot() +
geom_sf(data = tpoint_16_21,
aes(color = as.factor(yr))) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
theme_void()
What does it do?
CRS
Check the CRS of both sources? Are they different?
Solution
## Coordinate Reference System:
## User input: NAD83
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]]
## Coordinate Reference System:
## User input: NAD83
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]]
They are the same. So coord_sf() is not needed.