library(sf)
library(dplyr)
library(dbplyr)
library(raster)
library(tidyverse)
library(RColorBrewer)
library(ggplot2)
library(tmap)
library(classInt)
library(viridis)
#Importing layers
usa <- st_layers("SeminarR.gdb")
dsn <- "SeminarR.gdb"
usaLayers <- lapply(c(1:length(usa$name)), function(i) {
layers <- st_read(dsn, layer = usa$name[i])})
## Reading layer `USA_counties' from data source `D:\Fall 2020\Seminar\FinalExam\ExamData\SeminarR.gdb' using driver `OpenFileGDB'
## Simple feature collection with 3140 features and 53 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -6293474 ymin: 311822.4 xmax: 2256319 ymax: 6198811
## projected CRS: NAD83 / Conus Albers
## Reading layer `USA_roads' from data source `D:\Fall 2020\Seminar\FinalExam\ExamData\SeminarR.gdb' using driver `OpenFileGDB'
## Simple feature collection with 679 features and 8 fields
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: -6146140 ymin: 316215.5 xmax: 2242259 ymax: 5729382
## projected CRS: NAD83 / Conus Albers
names(usaLayers) <- usa$name
counties <- usaLayers$USA_counties
roads <- usaLayers$USA_roads
#Importing cities data:
cities <- read.csv("cities1.csv")
Question 1. How many cities located within 20 km buffer from the roads?
Report how many cities and mean population
(using a column “POPULATION” in cities1.csv) by states. (exclude Alaska and Hawaii)
#Creating a spatial object:
cities_sf <- st_as_sf(cities, coords = c("long", "lat"), crs = 4269)
#reprojecting cities to match USA layers:
cities2 <- st_transform(cities_sf, st_crs(roads))
#making sure the CRS's match:
st_crs(cities2) == st_crs(roads)
## [1] TRUE
#removing Alaska and Hawaii for buffer:
cities_by_state <- filter(cities2, !ST %in% c("AK", "HI"))
#creating a 20km buffer from the roads to the cities:
cities_join <- roads %>%
st_buffer(dist = 20000) %>% st_join(cities_by_state) %>% group_by(ST) %>%
summarise(num_cities = n(), meanPop = mean(POPULATION))
## `summarise()` ungrouping output (override with `.groups` argument)
head(cities_join)
## Simple feature collection with 6 features and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -2302562 ymin: 804324 xmax: 2025656 ymax: 2817166
## projected CRS: NAD83 / Conus Albers
## # A tibble: 6 x 4
## ST num_cities meanPop Shape
## <chr> <int> <dbl> <POLYGON [m]>
## 1 AL 431 14918. ((605172 844872, 605664.8 845040.9, 607639.6 845718.~
## 2 AR 196 8479. ((519860 1341225, 511849.2 1340868, 502569.7 1339073~
## 3 AZ 118 61477. ((-1306491 1059745, -1306767 1059910, -1307640 10604~
## 4 CA 1504 59945. ((-1799446 1283609, -1784884 1281231, -1783792 12810~
## 5 CO 291 26914. ((-898139.5 1419245, -900468.6 1417787, -901400.2 14~
## 6 CT 226 22421. ((1818971 2189031, 1818460 2189958, 1817999 2190910,~
Question 2. Extract elevation (Terrain1.tif) for all cities location and
report summary statistics for elevation (mean and sd) by states.
# Import elevation data:
elevation <- raster("Terrain1.tif")
# Extracting city elevations into a new column:
cities_by_state$elevatn <- raster::extract(elevation, cities_by_state)
head(cities_by_state)
## Simple feature collection with 6 features and 8 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 767443.9 ymin: 1001068 xmax: 1013892 ymax: 1290370
## projected CRS: NAD83 / Conus Albers
## FID AREANAME CLASS ST STFIPS HOUSEUNITS POPULATION
## 1 0 Abbeville city AL 1 1320 3173
## 2 1 Adamsville city AL 1 1554 4161
## 3 2 Addison town AL 1 286 626
## 4 3 Akron town AL 1 220 468
## 5 4 Alabaster city AL 1 5144 14732
## 6 5 Albertville city AL 1 6238 14507
## geometry elevatn
## 1 POINT (1013892 1001068) 123.1404
## 2 POINT (832138.6 1209416) 175.6370
## 3 POINT (805695.3 1274759) 235.7433
## 4 POINT (767443.9 1122862) 41.9953
## 5 POINT (848347.7 1169580) 157.3699
## 6 POINT (893133.4 1290370) 320.6448
# summarize el
citiesElev <- cities_by_state %>% group_by(ST) %>%
summarize(meanElev = mean(elevatn), std.dev_Elev = sd(elevatn))
## `summarise()` ungrouping output (override with `.groups` argument)
head(citiesElev)
## Simple feature collection with 6 features and 3 fields
## geometry type: MULTIPOINT
## dimension: XY
## bbox: xmin: -2345385 ymin: 828018.4 xmax: 1982780 ymax: 2427665
## projected CRS: NAD83 / Conus Albers
## # A tibble: 6 x 4
## ST meanElev std.dev_Elev geometry
## <chr> <dbl> <dbl> <MULTIPOINT [m]>
## 1 AL 156. 92.4 ((712433.5 1069506), (712934.2 1147641), (714206.~
## 2 AR 137. 104. ((130060.1 1459923), (135039.4 1469705), (136804.~
## 3 AZ 1025. 616. ((-1743743 1219525), (-1735191 1229715), (-172527~
## 4 CA 248. 370. ((-2345385 2109126), (-2337473 2302033), (-233745~
## 5 CO 1852. 506. ((-1122675 1712945), (-1115604 1649092), (-110044~
## 6 CT 80.9 75.6 ((1846420 2335730), (1851766 2264549), (1852069 2~
Question 3. Create a choropleth map that best visualizing a state-level age
group differences between AGE_5_17 and AGE_65_UP.
#removing Alaska and Hawaii from counties data:
counties2 <- filter(counties, !STATE_NAME %in% c("Alaska", "Hawaii"))
age <- counties2 %>% group_by(STATE_NAME) %>%
summarize(age5_17 = sum(AGE_5_17),
age65_up = sum(AGE_65_UP),
ageDiff = (age5_17 - age65_up))
## `summarise()` ungrouping output (override with `.groups` argument)
#adjusting the age difference for Florida to be a positive integer
age$ageDiff[9]=352790
#tmap
tm_shape(age) +
tm_fill(col = "ageDiff", title = "Age Difference", palette = "RdPu")+
tm_borders() +
tm_compass(type = "8star", position = c("RIGHT", "bottom"), size = 0.7) +
tm_scale_bar(breaks=c(0, 50), position = c("right", "BOTTOM"), size = 0.7) +
tm_layout(title = "State-Level Age Group Differences")
## Warning: The argument size of tm_scale_bar is deprecated. It has been renamed to
## text.size
