Libraries
# SPDS
library(tidyverse)
library(sf)
library(units)
# Data
library(USAboundaries)
library(rnaturalearthdata)
# Visualization
library(gghighlight)
library(ggrepel)
library(knitr)
library(readr)
Question 1
# 1.1- Define A Projection
eqdc = '+proj=eqdc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'
# 1.2- Get State Boundaries
states = USAboundaries::us_states(resolution = "low") %>%
filter(!stusps %in% c("PR", "HI", "AK")) %>%
st_transform(eqdc)
# 1.3 - Get country boundaries for Mexico, the United States of America, and Canada
worldboundaries = rnaturalearthdata::countries110 %>%
st_as_sf(worldboundaries) %>%
filter(admin %in% c("United States of America", "Canada", "Mexico")) %>%
st_transform(eqdc)
# 1.4 - Get city locations from the CSV file
uscities = readr::read_csv("../data/uscities.csv") %>%
st_as_sf(coords = c("lng", "lat"), crs = 4269) %>%
filter(!state_name %in% c("Hawaii", "Puerto Rico", "Alaska")) %>%
st_transform(eqdc)
##
## -- Column specification --------------------------------------------------------
## cols(
## city = col_character(),
## city_ascii = col_character(),
## state_id = col_character(),
## state_name = col_character(),
## county_fips = col_double(),
## county_name = col_character(),
## lat = col_double(),
## lng = col_double(),
## population = col_double(),
## density = col_double(),
## source = col_character(),
## military = col_logical(),
## incorporated = col_logical(),
## timezone = col_character(),
## ranking = col_double(),
## zips = col_character(),
## id = col_double()
## )
Question 2
# 2.1 - Distance to USA Border (coastline or national) (km)
USAboundaries::us_states() %>%
filter(!state_name %in% c("Puerto Rico", "Hawaii", "Alaska"))
## Simple feature collection with 49 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.7258 ymin: 24.49813 xmax: -66.9499 ymax: 49.38436
## Geodetic CRS: WGS 84
## First 10 features:
## statefp statens affgeoid geoid stusps name lsad
## 1 23 01779787 0400000US23 23 ME Maine 00
## 2 04 01779777 0400000US04 04 AZ Arizona 00
## 3 05 00068085 0400000US05 05 AR Arkansas 00
## 4 10 01779781 0400000US10 10 DE Delaware 00
## 5 13 01705317 0400000US13 13 GA Georgia 00
## 6 27 00662849 0400000US27 27 MN Minnesota 00
## 7 06 01779778 0400000US06 06 CA California 00
## 8 11 01702382 0400000US11 11 DC District of Columbia 00
## 9 12 00294478 0400000US12 12 FL Florida 00
## 10 16 01779783 0400000US16 16 ID Idaho 00
## aland awater state_name state_abbr jurisdiction_type
## 1 79885221885 11748755195 Maine ME state
## 2 294198560125 1027346486 Arizona AZ state
## 3 134771517596 2960191698 Arkansas AR state
## 4 5047194742 1398720828 Delaware DE state
## 5 149169848456 4741100880 Georgia GA state
## 6 206232257655 18929176411 Minnesota MN state
## 7 403501101370 20466718403 California CA state
## 8 158364992 18633403 District of Columbia DC district
## 9 138924199212 31386038155 Florida FL state
## 10 214042908012 2398669593 Idaho ID state
## geometry
## 1 MULTIPOLYGON (((-68.92401 4...
## 2 MULTIPOLYGON (((-114.7997 3...
## 3 MULTIPOLYGON (((-94.61792 3...
## 4 MULTIPOLYGON (((-75.77379 3...
## 5 MULTIPOLYGON (((-85.60516 3...
## 6 MULTIPOLYGON (((-97.22904 4...
## 7 MULTIPOLYGON (((-118.594 33...
## 8 MULTIPOLYGON (((-77.11976 3...
## 9 MULTIPOLYGON (((-81.81169 2...
## 10 MULTIPOLYGON (((-117.243 44...
border = USAboundaries::us_states() %>%
filter(!state_name %in% c("Puerto Rico", "Hawaii", "Alaska")) %>%
st_union %>%
st_cast("MULTILINESTRING")%>%
st_transform(eqdc)
plot(border)
distance1.1 = uscities %>%
mutate(distance_to_border = st_distance(uscities, border)) %>%
select(city, distance_to_border, state_name) %>%
arrange(-distance_to_border) %>%
slice(1:5)
knitr::kable(distance1.1, caption = "5 Cities Furthest From USA Border")
5 Cities Furthest From USA Border
| Dresden |
1012398 [m] |
Kansas |
POINT (-356590.2 -33214.32) |
| Herndon |
1007763 [m] |
Kansas |
POINT (-384537.3 41.28381) |
| Hill City |
1005143 [m] |
Kansas |
POINT (-311583.2 -63679.48) |
| Atwood |
1004754 [m] |
Kansas |
POINT (-405635.6 -9814.497) |
| Jennings |
1003646 [m] |
Kansas |
POINT (-346143.8 -27290.37) |
# 2.2 - Distance to States (km)
stateborders = USAboundaries::us_states() %>%
filter(!state_name %in% c("Puerto Rico", "Hawaii", "Alaska")) %>%
st_combine %>%
st_cast("MULTILINESTRING")%>%
st_transform(eqdc)
plot(stateborders)
distance2.1 = uscities %>%
mutate(distance_to_border = st_distance(uscities, border),
distance_to_border = set_units(distance_to_border, "km")) %>%
select(city, distance_to_border, state_name) %>%
arrange(-distance_to_border) %>%
slice(1:5)
knitr::kable(distance2.1, caption = "5 Cities Furthest From USA Border")
5 Cities Furthest From USA Border
| Dresden |
1012.398 [km] |
Kansas |
POINT (-356590.2 -33214.32) |
| Herndon |
1007.763 [km] |
Kansas |
POINT (-384537.3 41.28381) |
| Hill City |
1005.143 [km] |
Kansas |
POINT (-311583.2 -63679.48) |
| Atwood |
1004.754 [km] |
Kansas |
POINT (-405635.6 -9814.497) |
| Jennings |
1003.646 [km] |
Kansas |
POINT (-346143.8 -27290.37) |
# 2.3 - Distance to Mexico (km)
mexico = rnaturalearthdata::countries110 %>%
st_as_sf(worldboundaries) %>%
filter(admin %in% c("Mexico"))%>%
st_transform(eqdc)
distance3.1 = uscities %>%
mutate(distance_to_border = st_distance(uscities, mexico),
distance_to_border = set_units(distance_to_border, "km")) %>%
select(city, distance_to_border, state_name) %>%
arrange(-distance_to_border) %>%
slice(1:5)
knitr::kable(distance3.1, caption = "5 Cities Furthest From Mexican Border")
5 Cities Furthest From Mexican Border
| Caribou |
3250.334 [km] |
Maine |
POINT (1981398 1070430) |
| Presque Isle |
3234.570 [km] |
Maine |
POINT (1987753 1051523) |
| Calais |
3134.348 [km] |
Maine |
POINT (2093840 904408.7) |
| Passamaquoddy Pleasant Point |
3127.869 [km] |
Maine |
POINT (2111748 890089.8) |
| Eastport |
3125.624 [km] |
Maine |
POINT (2115970 885912.2) |
# 2.4 - Distance to Canada (km)
canada = rnaturalearthdata::countries110 %>%
st_as_sf(worldboundaries) %>%
filter(admin %in% c("Canada"))%>%
st_transform(eqdc)
distance4.1 = uscities %>%
mutate(distance_to_border = st_distance(uscities, canada),
distance_to_border = set_units(distance_to_border, "km")) %>%
select(city, distance_to_border, state_name) %>%
arrange(-distance_to_border) %>%
slice(1:5)
knitr::kable(distance4.1, caption = "5 Cities Furthest From Canadian Border")
5 Cities Furthest From Canadian Border
| Guadalupe Guerra |
2206.455 [km] |
Texas |
POINT (-298391.5 -1502124) |
| Fronton |
2204.790 [km] |
Texas |
POINT (-297784.3 -1500504) |
| Fronton Ranchettes |
2202.118 [km] |
Texas |
POINT (-293005.2 -1500698) |
| Evergreen |
2202.020 [km] |
Texas |
POINT (-292736.4 -1500773) |
| Ramos |
2201.882 [km] |
Texas |
POINT (-293014.6 -1500398) |
Question 3
# 3.1 Data
top10 = uscities %>%
select(city, population) %>%
slice_max(population, n = 10)
continents = rnaturalearthdata::countries110 %>%
st_as_sf(continents) %>%
filter(admin %in% c("United States of America", "Canada", "Mexico")) %>%
select(name, geometry) %>%
st_transform(eqdc)
plot(continents['name'], key.pos = 1)
combined_c = st_combine(continents)
ggplot() +
geom_sf(data = combined_c, lty = 3) +
geom_sf(data = stateborders, lty = 3) +
geom_sf(data = top10, lty = .5, col = "red") +
ggrepel::geom_label_repel(
data = top10,
aes(label = city, geometry = geometry),
stat = "sf_coordinates",
size = 3)
# 3.2 City Distance from the Border
allcities = uscities %>%
mutate(distance_to_border = st_distance(uscities, border),
distance_to_border = set_units(distance_to_border, "km")) %>%
select(city, distance_to_border, state_name) %>%
arrange(-distance_to_border)
ggplot()+
geom_sf(data = stateborders) +
geom_sf(data = allcities, aes(col = as.numeric(distance_to_border)))+
scale_colour_gradient(low = "#132B43",
high = "#56B1F7",
space = "Lab",
na.value = "grey50",
guide = "colourbar",
aesthetics = "colour") +
geom_sf(data = distance2.1, col = "red", lty = 3) +
ggrepel::geom_label_repel(
data = distance2.1,
aes(label = city, geometry = geometry),
stat = "sf_coordinates",
size = 3) +
labs(col = "Distance To Border",
title = "Cities Distance From Border",
x = "",
y = "") +
ggthemes::theme_map()
#3.3 City Distance from Nearest State
statecities = uscities %>%
mutate(distance_state_border = st_distance(uscities, stateborders),
distance_state_border = set_units(distance_state_border, "km")) %>%
select(city, distance_state_border, state_name) %>%
arrange(-distance_state_border)
ggplot()+
geom_sf(data = stateborders) +
geom_sf(data = statecities, aes(col = as.numeric(distance_state_border)))+
scale_colour_gradient(low = "99009", high = "663366") +
geom_sf(data = distance2.1, col = "red", lty = 3) +
ggrepel::geom_label_repel(
data = distance2.1,
aes(label = city, geometry = geometry),
stat = "sf_coordinates",
size = 3) +
labs(col = "Distance To State Border",
title = "Cities Distance To State Border",
x = "",
y = "") +
ggthemes::theme_map()
#3.4 Equidistance boundary from Mexico and Canada
equi = uscities %>%
mutate(equidistance = st_distance(canada, mexico),
equidistance = set_units(equidistance, "km")) %>%
select(city, equidistance, state_name) %>%
arrange(-equidistance)
ggplot()+
geom_sf(data = continents) +
geom_sf(data = equi, aes(col = as.numeric(equidistance)))+
scale_colour_gradient(low = "99009", high = "663366") +
geom_sf(data = distance2.1, col = "red", lty = 3) +
ggrepel::geom_label_repel(
data = distance2.1,
aes(label = city, geometry = geometry),
stat = "sf_coordinates",
size = 3) +
labs(col = "Distance To State Border",
title = "Equidistance",
x = "",
y = "") +
ggthemes::theme_map()
Question 4
#4.1 Quantifing Border Zone
# How many cities are in this 100 mile zone? (100 miles ~ 160 kilometers)
top100 = uscities %>%
mutate(distance_to_border = st_distance(uscities, border),
distance_to_border = set_units(distance_to_border, "km")) %>%
select(city, distance_to_border, state_name) %>%
arrange(distance_to_border) %>%
slice(1:12008)
count(top100)
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -2224960 ymin: -1601803 xmax: 2115970 ymax: 1275643
## CRS: +proj=eqdc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
## # A tibble: 1 x 2
## n geometry
## * <int> <MULTIPOINT [m]>
## 1 12008 ((-2224960 220719.4), (-2222534 227183.6), (-2219713 260735.3), (-22192~
# 12008
# How many people live in a city within 100 miles of the border?
pop100 = uscities %>%
mutate(distance_to_border = st_distance(uscities, border),
distance_to_border = set_units(distance_to_border, "km")) %>%
select(population, state_name, city)
count(pop100)
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -2224960 ymin: -1601803 xmax: 2115970 ymax: 1275643
## CRS: +proj=eqdc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
## # A tibble: 1 x 2
## n geometry
## * <int> <MULTIPOINT [m]>
## 1 27620 ((-2224960 220719.4), (-2222534 227183.6), (-2219713 260735.3), (-22192~
# 27620
centroid = st_centroid(states)
## Warning in st_centroid.sf(states): st_centroid assumes attributes are constant
## over geometries of x
ggplot()+
geom_sf(data = stateborders)+
geom_sf(data = top100, aes(col = as.numeric(distance_to_border)))+
# Extra Credit
geom_sf(data = centroid, col = "red", lty = 3)+
scale_color_gradient(low = "red",
high = "orange") +
labs(col = "Top 100 Cities in Danger Zone",
title = "Top 100 Cities",
x = "",
y = "") +
ggthemes::theme_map()
