Day 1: A few examples
ENS-215
Winter 2023
A few interesting examples
library(tidyverse)
library(stationaRy)
library(sf)
library(tmap)
library(spData)
library(AOI)
library(climateR)
library(raster)
library(rasterVis)
library(patchwork)Precip and temp for specified location (point)
# Get the data
#loc2use <- 'Death Valley National Park'
loc2use <- 'Union College'
AOI = AOI::geocode(loc2use,
pt = TRUE)
ts = getGridMET(AOI, param = c("tmax","tmin", "prcp"),
startDate = "2022-01-01",
endDate = "2022-12-31")# convert units
ts <- ts %>%
mutate(tmax = ((tmax - 273.15)*(9/5) +32),
tmin = ((tmin - 273.15)*(9/5) +32),
prcp = prcp/25.4
)# Summary table
ts %>%
summarize(pcrp_max = max(prcp),
temp_max = max(tmax),
temp_min = min(tmin),
temp_mean = mean((tmin+tmax)/2),
temp_max_change = max(tmax - tmin),
n_days = n()
)## pcrp_max temp_max temp_min temp_mean temp_max_change n_days
## 1 2.30315 96.53 -7.33 50.14252 46.26 365# Make the figures
fig_01 <- ts %>%
ggplot() +
geom_line(aes(x = date, y = tmax), color = "red") +
geom_line(aes(x = date, y = tmin), color = "blue") +
labs(x = "",
y = "Temperature (F)"
) +
geom_hline(yintercept = 32, linetype = "dashed") +
theme_bw()
fig_02 <- ts %>%
ggplot() +
geom_col(aes(x = date, y = prcp), color = "black") +
labs(y = "Precipitation (inches)") +
theme_bw()# Make the figures
(fig_01 / fig_02) +
plot_annotation(title = paste(loc2use,": Temperature and precipitation", sep = ""),
caption = "Data source: GridMET",
tag_levels = "a"
) 
Obtain and analyze precipitation data for Hurricane Henry
henri = getGridMET(aoi_get(state = c("NY")),
param = "prcp",
startDate = "2021-08-16", endDate = "2021-08-24")r = raster::stack(henri)r_inches <- r/25.4Static map
tmap_mode("plot")
fig_map <- r_inches %>%
tm_shape() +
tm_raster(style = "cont") +
tm_shape(spData::us_states) +
tm_borders()
fig_map
Interactive map
tmap_mode("view")
fig_map tmap_mode("plot")Get the maximum precipitation amount for each day of the storm
round(raster::maxValue(r_inches),2) %>% as.data.frame()## .
## 1 1.37
## 2 3.79
## 3 4.24
## 4 3.63
## 5 2.03
## 6 4.92
## 7 6.57
## 8 2.97
## 9 0.45Get climate data for specified region of the US
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched offtemperature_US = getGridMET(aoi_get(state = "conus"),
param = "tmin",
startDate = "2022-12-24", endDate = "2022-12-24")
#> Spherical geometry (s2) switched off
#> Spherical geometry (s2) switched ontemperature_US <- raster::stack(temperature_US)temperature_US <- temperature_US - 273.15temperature_US <- (temperature_US*9/5) + 32tmap_mode("plot")## tmap mode set to plottingtemperature_US %>%
tm_shape() +
tm_raster(style = "cont", palette = "-RdBu", midpoint = 32) +
tm_shape(spData::us_states) +
tm_borders()
Map of springs in New Mexico
library(osmdata)
library(tigris)
library(sf)
library(osmplotr)
library(tmaptools)
library(OpenStreetMap)
Sys.setenv(MAPBOX_API_KEY = "pk.eyJ1Ijoic3RhaGxtIiwiYSI6ImNrZnJiMDMxbDA0aGsyenFlajhvMzZ4bXUifQ.I7l7fJBAHCQWRwyYozq4ZQ")loc2use <- "New Mexico"
bb_values <- getbb(loc2use)
bb_values## min max
## x -109.05022 -103.00223
## y 31.33221 37.00015springs_data <- opq(bb_values) %>%
add_osm_feature(key = 'natural', value = 'spring') %>%
osmdata_sf()loc_border <- spData::us_states %>%
filter(NAME == "New Mexico")tmap_mode("view")## tmap mode set to interactive viewingmap_springs <-
tm_shape(loc_border) +
tm_borders(col = "black") +
tm_shape(springs_data$osm_points) +
tm_dots(col = "blue") map_springsObtain streamflow data from the USGS (Colorado River at Lee’s Ferry example)
library(dataRetrieval)
library(lubridate)df_stream_data <- readNWISdv(siteNumbers = "09380000",
parameterCd = c("00060"),
statCd = "00003") %>%
renameNWISColumns()Daily flows
df_stream_data %>%
ggplot(aes(x = Date, y = Flow)) +
geom_line() +
theme_classic()
Annual flow statistics
table_flows <- df_stream_data %>%
mutate(Year = year(Date)) %>%
group_by(Year) %>%
summarize(mean_flow = mean(Flow, na.rm= T),
min_flow = min(Flow, na.rm = T),
max_flow = max(Flow, na.rm = T),
n_meas = n()) %>%
filter(n_meas > 350)
table_flows## # A tibble: 101 x 5
## Year mean_flow min_flow max_flow n_meas
## <dbl> <dbl> <dbl> <dbl> <int>
## 1 1922 22196. 3700 116000 365
## 2 1923 23407. 4680 96200 365
## 3 1924 16099. 1000 72800 366
## 4 1925 17052. 1500 52300 365
## 5 1926 18041. 2500 84000 365
## 6 1927 24190. 2700 119000 365
## 7 1928 20265. 2700 113000 366
## 8 1929 27060. 4100 111000 365
## 9 1930 17119. 2900 71400 365
## 10 1931 8589. 1800 33500 365
## # ... with 91 more rowsfig_max <- table_flows %>%
ggplot(aes(x = Year)) +
geom_line(aes(y = max_flow), size = 1, color = "blue") +
#geom_line(aes(y = min_flow), size = 1, color = "red") +
#geom_line(aes(y = mean_flow), size = 1, color = "black") +
theme_classic()
fig_min <- table_flows %>%
ggplot(aes(x = Year)) +
geom_line(aes(y = min_flow), size = 1, color = "red") +
#geom_line(aes(y = mean_flow), size = 1, color = "black") +
theme_classic()
fig_mean <- table_flows %>%
ggplot(aes(x = Year)) +
geom_line(aes(y = mean_flow), size = 1, color = "black") +
theme_classic()fig_mean/fig_max/fig_min
Obtain US Census data
https://api.census.gov/data/2015/acs/acs5/variables
library(tidycensus)census_example_df <- get_acs(geography = "county",
variables = c(medincome = "B19013_001"),
state = "NY",
year = 2021)## Getting data from the 2017-2021 5-year ACScensus_example_df %>%
arrange(-estimate)## # A tibble: 62 x 5
## GEOID NAME variable estimate moe
## <chr> <chr> <chr> <dbl> <dbl>
## 1 36059 Nassau County, New York medincome 126576 1363
## 2 36103 Suffolk County, New York medincome 111660 1368
## 3 36079 Putnam County, New York medincome 111617 4424
## 4 36119 Westchester County, New York medincome 105387 1643
## 5 36087 Rockland County, New York medincome 99707 2498
## 6 36061 New York County, New York medincome 93956 1493
## 7 36091 Saratoga County, New York medincome 90800 2122
## 8 36085 Richmond County, New York medincome 89427 2247
## 9 36027 Dutchess County, New York medincome 87112 2254
## 10 36071 Orange County, New York medincome 85640 2000
## # ... with 52 more rowscensus_example_df %>%
mutate(NAME = gsub(" County, New York", "", NAME)) %>%
ggplot(aes(x = estimate, y = reorder(NAME, estimate))) +
geom_errorbarh(aes(xmin = estimate - moe, xmax = estimate + moe)) +
geom_point(color = "red", size = 3) +
labs(title = "Household income by county in New York",
subtitle = "2021 American Community Survey",
y = "",
x = "ACS estimate (bars represent margin of error)") +
theme_bw()
Schdy <- get_acs(
state = "NY",
county = "Schenectady",
geography = "tract",
variables = "B19013_001",
geometry = TRUE,
year = 2020
)##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|============================ | 40%
|
|============================= | 41%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|===================================== | 52%
|
|===================================== | 53%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================= | 64%
|
|============================================== | 65%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================= | 70%
|
|================================================== | 71%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 85%
|
|============================================================ | 86%
|
|============================================================= | 88%
|
|============================================================== | 89%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|======================================================================| 100%tmap_mode("plot")## tmap mode set to plottingmap_income <- tm_shape(Schdy) +
tm_polygons(col = "estimate", alpha = 1)
map_income
map_income_interactive <- tm_shape(Schdy) +
tm_polygons(col = "estimate", alpha = 0.6)tmap_mode("view")## tmap mode set to interactive viewingmap_income_interactive