OrtizLab4
Joe Ortiz
2025-02-22
Lesson 4 Part A
The first part of this lab examines point and line features representing tornado touchdown locations and travel paths respectively. The simple features and ggplot libraries were used to map and visualize the tornado dataset. Figure 1 plots a composite of touchdown points and paths symbolized with varying color by for each year from 2016-2021. Figure 2 maps tornado touchdown density by county for each year. Figure 3 shows an attempt at setting custom symbology breaks. This attempt was unsuccessful.
#Load package libraries
library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(RColorBrewer)
library(units)
library(cowplot)library(here)## here() starts at C:/PSU/GEOG588/WS_Lesson4#Read in data
okcounty <- st_read(here("data", "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(here("data", "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(here("data", "ok_tornado_path.shp"), quiet = TRUE)#Map tornado paths & points then save as objects
pathmap <- ggplot() +
geom_sf(data = tpath,
aes(color = as.factor(yr)), show.legend = FALSE) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()
pointmap <- ggplot() +
geom_sf(data = tpoint,
aes(color = as.factor(yr)), show.legend = FALSE) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()#Create composite plot
plot_grid(pathmap, pointmap,
labels = c("A) Tornado paths 1953-2021",
"B) Tornado points 1953-2021",
label_size = 12),
ncol = 1,
hjust = 0,
label_x = 0,
align = "hv")
Figure 1.) A composite plot of tornado touchdown
points and paths with year values from 2016-2021 symbolized with unique
color grouping. Here map legends are hidden for aesthetics.
#Iterate over annual data, creating a density layer for each year from 2016 to 2021
# filter year 2016
tpoint_16 <- tpoint %>%
filter(yr == 2016) %>%
select(om, yr, date)
# summarize county w/ points for 2016 (will have no data values)
countypnt16 <- st_join(tpoint_16, okcounty)
countypnt16 <- st_drop_geometry(countypnt16)
countysum16 <- countypnt16 %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# build county layer summarized for 2016 (fill no data values from previous step)
countymap16 <- okcounty %>%
left_join(countysum16, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
# append year column lost during summary with appropriate value
countymap16['yr'] <- '2016'
# filter year 2017
tpoint_17 <- tpoint %>%
filter(yr == 2017) %>%
select(om, yr, date)
# summarize county w/ points for 2017 (will have no data values)
countypnt17 <- st_join(tpoint_17, okcounty)
countypnt17 <- st_drop_geometry(countypnt17)
countysum17 <- countypnt17 %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# build county layer summarized for 2017 (fill no data values from previous step)
countymap17 <- okcounty %>%
left_join(countysum17, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
# append year column lost during summary with appropriate value
countymap17['yr'] <- '2017'
# filter year 2018
tpoint_18 <- tpoint %>%
filter(yr == 2018) %>%
select(om, yr, date)
# summarize county w/ points for 2018 (will have no data values)
countypnt18 <- st_join(tpoint_18, okcounty)
countypnt18 <- st_drop_geometry(countypnt18)
countysum18 <- countypnt18 %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# build county layer summarized for 2018 (fill no data values from previous step)
countymap18 <- okcounty %>%
left_join(countysum18, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
# append year column lost during summary with appropriate value
countymap18['yr'] <- '2018'
# filter year 2019
tpoint_19 <- tpoint %>%
filter(yr == 2019) %>%
select(om, yr, date)
# summarize county w/ points for 2019 (will have no data values)
countypnt19 <- st_join(tpoint_19, okcounty)
countypnt19 <- st_drop_geometry(countypnt19)
countysum19 <- countypnt19 %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# build county layer summarized for 2019 (fill no data values from previous step)
countymap19 <- okcounty %>%
left_join(countysum19, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
# append year column lost during summary with appropriate value
countymap19['yr'] <- '2019'
# filter year 2020
tpoint_20 <- tpoint %>%
filter(yr == 2020) %>%
select(om, yr, date)
# summarize county w/ points for 2020 (will have no data values)
countypnt20 <- st_join(tpoint_20, okcounty)
countypnt20 <- st_drop_geometry(countypnt20)
countysum20 <- countypnt20 %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# build county layer summarized for 2020 (fill no data values from previous step)
countymap20 <- okcounty %>%
left_join(countysum20, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
# append year column lost during summary with appropriate value
countymap20['yr'] <- '2020'
# filter year 2021
tpoint_21 <- tpoint %>%
filter(yr == 2021) %>%
select(om, yr, date)
# summarize county w/ points for 2021 (will have no data values)
countypnt21 <- st_join(tpoint_21, okcounty)
countypnt21 <- st_drop_geometry(countypnt21)
countysum21 <- countypnt21 %>%
group_by(GEOID) %>%
summarize(tcnt = n())
# build county layer summarized for 2021 (fill no data values from previous step)
countymap21 <- okcounty %>%
left_join(countysum21, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
# append year column lost during summary with appropriate value
countymap21['yr'] <- '2021'
# compile county layers using rbind()
countymap <- rbind(countymap16, countymap17, countymap18, countymap19, countymap20, countymap21)
# sum dens by county & year
ggplot(data = countymap) +
geom_sf(aes(fill = tdens)) +
facet_wrap(vars(yr), ncol = 2) +
coord_sf(datum = NA) +
theme_void()
Figure 2.) Here tornado touchdown points are joined with county boundaries to create a summarized density by county and year.
#Create 4 choropleth maps #First break manipulation
# filter year range ----
tpoint_16_21 <- tpoint %>%
filter(yr >= 2016 & yr <= 2021) %>%
select(om, yr, date)
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
group_by(GEOID) %>%
summarize(tcnt = n())
countymap <- okcounty %>%
left_join(countysum, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
numclas <- 4
qbrks <- seq(0, 1, length.out = numclas + 1)
qbrks## [1] 0.00 0.25 0.50 0.75 1.00countymap <- countymap %>%
mutate(tdens_c1 = cut(tdens,
breaks = quantile(tdens, breaks = qbrks),
include.lowest = T))
densTry1 <- ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c1)) +
scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd") +
theme_void() +
theme(legend.position = "bottom")#Second manipulation
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
group_by(GEOID) %>%
summarize(tcnt = n())
countymap <- okcounty %>%
left_join(countysum, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
numclas <- 6
qbrks <- seq(0, .5, length.out = numclas + 1)
qbrks## [1] 0.00000000 0.08333333 0.16666667 0.25000000 0.33333333 0.41666667 0.50000000countymap <- countymap %>%
mutate(tdens_c1 = cut(tdens,
breaks = quantile(tdens, breaks = qbrks),
include.lowest = T))
densTry2 <- ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c1)) +
scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd") +
theme_void() +
theme(legend.position = "bottom")#Plot density attempts
plot_grid(densTry1, densTry2,
labels = c("Density Attempt 1",
"Density Attempt 2",
label_size = 12),
ncol = 1,
hjust = 0,
label_x = 0,
align = "hv")
Figure 3.) This plot attempts to show variation
among breaks. In each attempt, the number class and quantile break
values were manipulated but I did not observe variation within the
plots.
Lesson 4 Part B
In the second part of this lab, point features from the California Department of Transportation (CALTRANS) are mapped using leaflet library. The point features represent highway rest areas across California. The underlying basemap is Open Street Map and points are represented by transparent circle markers with hover labels and popup icons.
#Load libraries
library(leaflet)
library(sf)
library(tidyverse)#Create map
calTransRestStops = read_csv(here("data", "Rest_Areas.csv")) %>%
st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs = 4326)
leaflet(data = calTransRestStops) %>%
setView(lng= -120.1157639, lat= 36.7416723, zoom = 7) %>%
addTiles() %>%
addCircleMarkers(popup = ~REST_AREA_NAME, label = ~CITY, radius = 4, color = "#FFA500") %>%
addLegend(labels = "CALTRANS Rest Area", colors = "#FFA500")
Figure 4.) A leaflet map showing CALTRANS
highway rest areas with hover labels showing local city names and popups
showing facility names.
Lesson 4 Part C
The final portion of this lab consists of exposure to Coordinate Reference Systems (CRS) in Michael C. Wimberly’s 2023 book, Geographic Data Science with R.
#Load libraries
library(terra)
library(colorspace)# read in data ----
county <- st_read(here("data", "cb_2018_us_county_20m.shp"), quiet = TRUE)
county <- county %>%
mutate(state = as.numeric(as.character(STATEFP))) %>%
filter(state != 2, state != 15, state < 60)
st_crs(county)## 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]]# plot map
ggplot(data = county) +
geom_sf(fill = NA) +
theme_bw()
# select key elements
st_crs(county)$Name## [1] "NAD83"st_crs(county)$epsg## [1] 4269st_crs(county)$proj4string## [1] "+proj=longlat +datum=NAD83 +no_defs"writeLines(st_crs(county)$WktPretty)## GEOGCS["NAD83",
## DATUM["North_American_Datum_1983",
## SPHEROID["GRS 1980",6378137,298.257222101]],
## PRIMEM["Greenwich",0],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AXIS["Latitude",NORTH],
## AXIS["Longitude",EAST],
## AUTHORITY["EPSG","4269"]]# reproject
county_aea <- st_transform(county, 5070)
writeLines(st_crs(county_aea)$WktPretty)## PROJCS["NAD83 / Conus Albers",
## GEOGCS["NAD83",
## DATUM["North_American_Datum_1983",
## SPHEROID["GRS 1980",6378137,298.257222101]],
## PRIMEM["Greenwich",0],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4269"]],
## PROJECTION["Albers_Conic_Equal_Area"],
## PARAMETER["latitude_of_center",23],
## PARAMETER["longitude_of_center",-96],
## PARAMETER["standard_parallel_1",29.5],
## PARAMETER["standard_parallel_2",45.5],
## PARAMETER["false_easting",0],
## PARAMETER["false_northing",0],
## UNIT["metre",1],
## AXIS["Easting",EAST],
## AXIS["Northing",NORTH],
## AUTHORITY["EPSG","5070"]]ggplot(data = county_aea) +
geom_sf(fill = NA) +
theme_bw()
# manually enter CRS WKT
badcrs_wkt <-
'PROJCS["BadAlbers",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["latitude_of_center",37.5],
PARAMETER["longitude_of_center",-96],
PARAMETER["standard_parallel_1",75],
PARAMETER["standard_parallel_2",80],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]'
mybadcrs <- st_crs(badcrs_wkt)
county_bad <- st_transform(county, mybadcrs)
ggplot(data = county_bad) +
geom_sf(fill = NA) +
theme_bw()
# reproject using existing file CRS
county_fixed <- st_transform(county_bad, st_crs(county_aea))
ggplot(data = county_fixed) +
geom_sf(fill = NA) +
theme_bw()
# NY example
newyork <- filter(county_aea, state == 36)
ggplot(data = newyork) +
geom_sf(fill = NA) +
theme_bw()
ny_utm <- st_transform(newyork, 26918)
writeLines(st_crs(ny_utm)$WktPretty)## PROJCS["NAD83 / UTM zone 18N",
## GEOGCS["NAD83",
## DATUM["North_American_Datum_1983",
## SPHEROID["GRS 1980",6378137,298.257222101]],
## PRIMEM["Greenwich",0],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4269"]],
## PROJECTION["Transverse_Mercator"],
## PARAMETER["latitude_of_origin",0],
## PARAMETER["central_meridian",-75],
## PARAMETER["scale_factor",0.9996],
## PARAMETER["false_easting",500000],
## PARAMETER["false_northing",0],
## UNIT["metre",1],
## AXIS["Easting",EAST],
## AXIS["Northing",NORTH],
## AUTHORITY["EPSG","26918"]]ggplot(data = ny_utm) +
geom_sf(fill = NA) +
theme_bw()
# load raster
landcov <- rast(here("data", "NLCD_2016_Land_Cover_Walton.tiff"))
nrow(landcov)## [1] 1579ncol(landcov)## [1] 1809res(landcov)## [1] 30 30ext(landcov)[1:4]## xmin xmax ymin ymax
## 1096815 1151085 1237395 1284765writeLines(crs(landcov))## PROJCRS["Albers_Conical_Equal_Area",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101004,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["Albers Equal Area",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",23,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-96,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",29.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",45.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["northing",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]crs(landcov, describe = TRUE)$name## [1] "Albers_Conical_Equal_Area"writeLines(st_crs(landcov)$WktPretty)## PROJCS["Albers_Conical_Equal_Area",
## GEOGCS["NAD83",
## DATUM["North_American_Datum_1983",
## SPHEROID["GRS 1980",6378137,298.257222101004]],
## PRIMEM["Greenwich",0],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4269"]],
## PROJECTION["Albers_Conic_Equal_Area"],
## PARAMETER["latitude_of_center",23],
## PARAMETER["longitude_of_center",-96],
## PARAMETER["standard_parallel_1",29.5],
## PARAMETER["standard_parallel_2",45.5],
## PARAMETER["false_easting",0],
## PARAMETER["false_northing",0],
## UNIT["metre",1,
## AUTHORITY["EPSG","9001"]],
## AXIS["Easting",EAST],
## AXIS["Northing",NORTH]]# reclassify raster ----
oldclas <- unique(landcov)
newclas <- c(1, 2, 2, 2, 2, 3, 4, 4, 4, 5, 5, 5, 6, 7, 7)
lookup <- data.frame(oldclas, newclas)
landcov_rc <- classify(landcov, lookup)
newnames <- c("Water",
"Developed",
"Barren",
"Forest",
"GrassShrub",
"Cropland",
"Wetland")
newcols <- c("mediumblue",
"red2",
"gray60",
"darkgreen",
"yellow2",
"orange4",
"paleturquoise2")
newcols2 <- desaturate(newcols, amount = 0.3)
source(here("data", "rasterdf.R"))
landcov_df <- rasterdf(landcov_rc)
ggplot(data = landcov_df) +
geom_raster(aes(x = x,
y = y,
fill = as.character(value))) +
scale_fill_manual(name = "Land cover",
values = newcols2,
labels = newnames,
na.translate = FALSE) +
coord_sf(expand = FALSE) +
theme_void()
# reproject raster
landcov_utm <- project(landcov_rc, "epsg:26917",
method = "near",
res = 30)
crs(landcov_utm, describe = TRUE)$name## [1] "NAD83 / UTM zone 17N"# explore
nrow(landcov_utm)## [1] 1828ncol(landcov_utm)## [1] 2047res(landcov_utm)## [1] 30 30# plot reprojected raster
landcovutm_df <- rasterdf(landcov_utm)
ggplot(data = landcovutm_df) +
geom_raster(aes(x = x,
y = y,
fill = as.character(value))) +
scale_fill_manual(name = "Land cover",
values = newcols2,
labels = newnames,
na.translate = FALSE) +
coord_sf(expand = FALSE) +
theme_void()## Warning: Removed 883216 rows containing missing values or values outside the scale range
## (`geom_raster()`).
# reproject raster from existing raster projection
landcov_goback <- project(landcov_utm,
landcov,
method = "near")
crs(landcov_goback, describe = TRUE)$name## [1] "Albers_Conical_Equal_Area"crs(landcov_rc, describe = TRUE)$name## [1] "Albers_Conical_Equal_Area"compareGeom(landcov_rc, landcov_goback)## [1] TRUE# examine raster for changes post reproject
freq(landcov_rc)## layer value count
## 1 1 1 37165
## 2 1 2 687865
## 3 1 3 7928
## 4 1 4 1296342
## 5 1 5 691107
## 6 1 6 6629
## 7 1 7 129375freq(landcov_goback)## layer value count
## 1 1 1 37182
## 2 1 2 687909
## 3 1 3 7916
## 4 1 4 1296069
## 5 1 5 691183
## 6 1 6 6631
## 7 1 7 129411# change matrices ----
projcomp <- crosstab(c(landcov_rc, landcov_goback))
shortnames <- c("Wat",
"Dev",
"Bare",
"For",
"Grass",
"Crop",
"Wet")
rownames(projcomp) <- shortnames
colnames(projcomp) <- shortnames