# Load in packages
library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.7.71
library(tidyterra)
## Warning: package 'tidyterra' was built under R version 4.3.3
##
## Attaching package: 'tidyterra'
## The following object is masked from 'package:stats':
##
## filter
library(sp)
## Warning: package 'sp' was built under R version 4.3.2
library(ggplot2)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.3.2
library(urbnmapr)
library(sf)
## Warning: package 'sf' was built under R version 4.3.2
## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE
# Map 1: Raster Map
#Load in raster data
setwd("C:/Users/kajsa/OneDrive/Documents/Binghamton/Senior - Spring/DIDA 370/NE1_50M_SR_W")
natural_earth <- rast("NE1_50M_SR_W.tif")
natural_earth <- terra::project(natural_earth, "EPSG:32116", method = "bilinear")
##
|---------|---------|---------|---------|
=========================================
#Inspect raster object
natural_earth
## class : SpatRaster
## dimensions : 32498, 20166, 3 (nrow, ncol, nlyr)
## resolution : 1656.569, 1656.569 (x, y)
## extent : -16453285, 16953079, -31346642, 22488528 (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / New York Central (EPSG:32116)
## source : spat_6088af1f55_24712.tif
## names : NE1_50M_SR_W_1, NE1_50M_SR_W_2, NE1_50M_SR_W_3
## min values : 62.21814, 90.4782, 92.21101
## max values : 255.00000, 255.0000, 255.00000
#Plot raster object
plot(natural_earth)
counties <- get_urbn_map("counties", sf = TRUE)
#filter the data to get just NYS
counties <- counties %>%
filter(state_abbv == "NY") %>%
st_transform("EPSG:32116")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## Warning in CPL_crs_from_input(x): GDAL Message 1: CRS EPSG:2163 is deprecated.
## Its non-deprecated replacement EPSG:9311 will be used instead. To use the
## original CRS, set the OSR_USE_NON_DEPRECATED configuration option to NO.
# plot the data
ggplot()+
geom_spatraster(data = natural_earth)+
geom_sf(counties, mapping = aes(), fill = 'transparent', color = "gray")+
theme_map()+
scale_fill_whitebox_c(
palette = "soft",
na.value = "white")+
labs(title = "Landcover in NYS",
fill = "Fill")+
theme(plot.title = element_text(hjust = 0.5, size = 16),
legend.position="right")
## <SpatRaster> resampled to 501084 cells for plotting
## ! `tidyterra::geom_spatraster()`: Plotting 3 overlapping layers: NE1_50M_SR_W_1, NE1_50M_SR_W_2, and NE1_50M_SR_W_3. Either:
## Use `facet_wrap(~lyr)` for faceting or
## Use `aes(fill = <name_of_layer>)` for displaying single layers
#filter the data to get just NYS
counties <- counties %>%
filter(state_abbv == "NY") %>%
st_transform("EPSG:32116")
rast <- crop(natural_earth, counties)
ggplot()+
geom_spatraster(data = rast)+
geom_sf(counties, mapping = aes(), fill = 'transparent', color = "white")+
theme_map()+
scale_fill_whitebox_c(
palette = "atlas",
na.value = "white")+
labs(title = "Landcover in New York State",
fill = "Fill")+
theme(plot.title = element_text(hjust = 0.5, size = 16),
legend.position="right")
## ! `tidyterra::geom_spatraster()`: Plotting 3 overlapping layers: NE1_50M_SR_W_1, NE1_50M_SR_W_2, and NE1_50M_SR_W_3. Either:
## Use `facet_wrap(~lyr)` for faceting or
## Use `aes(fill = <name_of_layer>)` for displaying single layers
# Map 2: Leaflet Map
# Load in packages
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.3.3
library(sf)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
##
## intersect, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(spData)
## Warning: package 'spData' was built under R version 4.3.2
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(urbnmapr)
library(ggplot2)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:terra':
##
## extract
counties <- get_urbn_map("counties", sf = TRUE)
counties <- counties %>%
filter(state_abbv == "NY") %>%
st_transform("EPSG:32116")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
counties <- counties %>%
filter(state_abbv == "NY") %>%
st_transform("EPSG:32116")
counties
## Simple feature collection with 62 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13302.14 ymin: 57720 xmax: 647370.2 ymax: 561701.5
## Projected CRS: NAD83 / New York Central
## First 10 features:
## county_fips state_abbv state_fips county_name fips_class state_name
## 1 36021 NY 36 Columbia County H1 New York
## 2 36087 NY 36 Rockland County H1 New York
## 3 36101 NY 36 Steuben County H1 New York
## 4 36123 NY 36 Yates County H1 New York
## 5 36033 NY 36 Franklin County H1 New York
## 6 36093 NY 36 Schenectady County H1 New York
## 7 36035 NY 36 Fulton County H1 New York
## 8 36049 NY 36 Lewis County H1 New York
## 9 36003 NY 36 Allegany County H1 New York
## 10 36057 NY 36 Montgomery County H1 New York
## geometry
## 1 MULTIPOLYGON (((470186.4 23...
## 2 MULTIPOLYGON (((447233.3 12...
## 3 MULTIPOLYGON (((155545 2676...
## 4 MULTIPOLYGON (((185775.2 29...
## 5 MULTIPOLYGON (((396740 5568...
## 6 MULTIPOLYGON (((438010.2 30...
## 7 MULTIPOLYGON (((397382.1 34...
## 8 MULTIPOLYGON (((309716 4316...
## 9 MULTIPOLYGON (((108241.8 28...
## 10 MULTIPOLYGON (((400329.3 33...
# download NYS income data
setwd("C:/Users/kajsa/OneDrive/Documents/Binghamton/Senior - Spring/DIDA 370/Assignment - 3")
income <- read.csv("personal_income_per_captia.csv")
merged_data <- left_join(counties, income, by = c("county_name" = "NYS.County"))
merged_data
## Simple feature collection with 62 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -13302.14 ymin: 57720 xmax: 647370.2 ymax: 561701.5
## Projected CRS: NAD83 / New York Central
## First 10 features:
## county_fips state_abbv state_fips county_name fips_class state_name
## 1 36021 NY 36 Columbia County H1 New York
## 2 36087 NY 36 Rockland County H1 New York
## 3 36101 NY 36 Steuben County H1 New York
## 4 36123 NY 36 Yates County H1 New York
## 5 36033 NY 36 Franklin County H1 New York
## 6 36093 NY 36 Schenectady County H1 New York
## 7 36035 NY 36 Fulton County H1 New York
## 8 36049 NY 36 Lewis County H1 New York
## 9 36003 NY 36 Allegany County H1 New York
## 10 36057 NY 36 Montgomery County H1 New York
## Income geometry
## 1 53436 MULTIPOLYGON (((470186.4 23...
## 2 58133 MULTIPOLYGON (((447233.3 12...
## 3 43528 MULTIPOLYGON (((155545 2676...
## 4 38601 MULTIPOLYGON (((185775.2 29...
## 5 36617 MULTIPOLYGON (((396740 5568...
## 6 49389 MULTIPOLYGON (((438010.2 30...
## 7 42020 MULTIPOLYGON (((397382.1 34...
## 8 42814 MULTIPOLYGON (((309716 4316...
## 9 33942 MULTIPOLYGON (((108241.8 28...
## 10 41070 MULTIPOLYGON (((400329.3 33...
counties.sf = st_transform(merged_data, "EPSG:4326")
leaflet() %>%
leaflet::addPolygons(data = counties.sf, color = "green")
p_popup <- paste0("<strong>income: </strong>", merged_data$Income)
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)
pal_fun
## function (x)
## {
## binsToUse <- quantile(x, probs, na.rm = TRUE, names = FALSE)
## ints <- cut(x, binsToUse, labels = FALSE, include.lowest = TRUE,
## right = right)
## if (any(is.na(x) != is.na(ints)))
## warning("Some values were outside the color scale and will be treated as NA")
## colorFunc(ints)
## }
## <bytecode: 0x000001d69b3fafc0>
## <environment: 0x000001d69b3fa0a8>
## attr(,"colorType")
## [1] "quantile"
## attr(,"colorArgs")
## attr(,"colorArgs")$probs
## [1] 0.0 0.2 0.4 0.6 0.8 1.0
##
## attr(,"colorArgs")$na.color
## [1] "#808080"
leaflet() %>%
leaflet::addPolygons(data = counties.sf,
stroke = FALSE,
fillColor = ~pal_fun(Income),
fillOpacity = 1,
smoothFactor = 0.5,
popup = p_popup) %>%
addLegend("bottomright",
pal=pal_fun,
values=counties.sf$inequality,
title = 'Income Levels') %>%
addProviderTiles(providers$CartoDB.Positron)