#install.packages("leaflet")
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.5.3
# Install packages once if needed:
# install.packages(c("sf", "tidyverse", "leaflet", "htmlwidgets"))
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.5.3
## Warning: package 'dplyr' was built under R version 4.5.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(leaflet)
library(htmlwidgets)
leaflet() |>
addTiles() |>
setView(
lng = -115.55,
lat = 33.00,
zoom = 9
)
#install.packages("terra")
library(terra)
## Warning: package 'terra' was built under R version 4.5.3
## terra 1.9.34
##
## Attaching package: 'terra'
## The following object is masked from 'package:tidyr':
##
## extract
#install.packages("geodata")
library(geodata)
## Warning: package 'geodata' was built under R version 4.5.3
library(geodata)
# Download the Imperial County boundary
imperial <- gadm(
country = "USA",
level = 2,
path = "data"
)
imperial <- imperial[imperial$NAME_2 == "Imperial", ]
imperial
## class : SpatVector
## geometry : polygons
## dimensions : 1, 13 (geometries, attributes)
## extent : -116.1063, -114.4603, 32.61827, 33.43729 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : GID_2 GID_0 COUNTRY GID_1 NAME_1 NL_NAME_1 NAME_2 VARNAME_2 NL_NAME_2 TYPE_2 (and 3 more)
## type : <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## values : USA.5.13_1 USA United States USA.5_1 California NA Imperial NA NA County
library(sf)
imperial_sf <- st_as_sf(imperial)
plot(imperial_sf)
## Warning: plotting the first 9 out of 13 attributes; use max.plot = 13 to plot
## all
# Get monthly TerraClimate maximum temperature data
library(geodata)
climate <- worldclim_global(
var = "tmax",
res = 10,
path = "data"
)
library(terra)
# Check what we currently have available
list.files("data", recursive = TRUE)
## [1] "climate/wc2.1_10m/wc2.1_10m_tmax_01.tif"
## [2] "climate/wc2.1_10m/wc2.1_10m_tmax_02.tif"
## [3] "climate/wc2.1_10m/wc2.1_10m_tmax_03.tif"
## [4] "climate/wc2.1_10m/wc2.1_10m_tmax_04.tif"
## [5] "climate/wc2.1_10m/wc2.1_10m_tmax_05.tif"
## [6] "climate/wc2.1_10m/wc2.1_10m_tmax_06.tif"
## [7] "climate/wc2.1_10m/wc2.1_10m_tmax_07.tif"
## [8] "climate/wc2.1_10m/wc2.1_10m_tmax_08.tif"
## [9] "climate/wc2.1_10m/wc2.1_10m_tmax_09.tif"
## [10] "climate/wc2.1_10m/wc2.1_10m_tmax_10.tif"
## [11] "climate/wc2.1_10m/wc2.1_10m_tmax_11.tif"
## [12] "climate/wc2.1_10m/wc2.1_10m_tmax_12.tif"
## [13] "gadm/gadm41_USA_2_pk.rds"
#install.packages("daymetr")
library(daymetr)
## Warning: package 'daymetr' was built under R version 4.5.3
imperial_lat <- 33.04
imperial_lon <- -115.48
daymet_data <- download_daymet(
lat = imperial_lat,
lon = imperial_lon,
start = 2005,
end = 2021,
internal = TRUE
)
## Downloading DAYMET data for: Daymet at 33.04/-115.48 latitude/longitude !
## Done !
head(daymet_data)
## $site
## [1] "Daymet"
##
## $tile
## [1] 11013
##
## $latitude
## [1] 33.04
##
## $longitude
## [1] -115.48
##
## $altitude
## [1] -46
##
## $tile
## [1] 11013
head(daymet_data)
## $site
## [1] "Daymet"
##
## $tile
## [1] 11013
##
## $latitude
## [1] 33.04
##
## $longitude
## [1] -115.48
##
## $altitude
## [1] -46
##
## $tile
## [1] 11013
names(daymet_data)
## [1] "site" "tile" "latitude" "longitude" "altitude" "tile"
## [7] "data"
head(daymet_data$data)
## year yday dayl..s. prcp..mm.day. srad..W.m.2. swe..kg.m.2. tmax..deg.c.
## 1 2005 1 35507.18 0.00 267.91 0 18.85
## 2 2005 2 35538.23 0.00 283.30 0 18.89
## 3 2005 3 35571.82 6.14 103.33 0 13.58
## 4 2005 4 35607.95 2.37 125.44 0 13.82
## 5 2005 5 35646.58 0.00 283.98 0 15.95
## 6 2005 6 35687.68 0.00 260.39 0 14.94
## tmin..deg.c. vp..Pa.
## 1 7.01 361.47
## 2 6.02 336.01
## 3 8.11 940.23
## 4 7.13 840.63
## 5 3.18 357.33
## 6 3.85 418.26
names(daymet_data$data)
## [1] "year" "yday" "dayl..s." "prcp..mm.day."
## [5] "srad..W.m.2." "swe..kg.m.2." "tmax..deg.c." "tmin..deg.c."
## [9] "vp..Pa."
library(tidyverse)
annual_temp <- daymet_data$data |>
group_by(year) |>
summarise(
mean_tmax = mean(tmax..deg.c., na.rm = TRUE)
)
annual_temp
## # A tibble: 17 × 2
## year mean_tmax
## <int> <dbl>
## 1 2005 32.0
## 2 2006 32.5
## 3 2007 32.3
## 4 2008 32.3
## 5 2009 32.3
## 6 2010 31.3
## 7 2011 32.0
## 8 2012 32.9
## 9 2013 32.4
## 10 2014 33.1
## 11 2015 32.9
## 12 2016 33.1
## 13 2017 33.1
## 14 2018 32.7
## 15 2019 31.2
## 16 2020 33.4
## 17 2021 33.2
annual_temp <- annual_temp |>
mutate(
mean_tmax_F = mean_tmax * 9/5 + 32
)
annual_temp
## # A tibble: 17 × 3
## year mean_tmax mean_tmax_F
## <int> <dbl> <dbl>
## 1 2005 32.0 89.7
## 2 2006 32.5 90.5
## 3 2007 32.3 90.2
## 4 2008 32.3 90.2
## 5 2009 32.3 90.1
## 6 2010 31.3 88.4
## 7 2011 32.0 89.5
## 8 2012 32.9 91.2
## 9 2013 32.4 90.2
## 10 2014 33.1 91.7
## 11 2015 32.9 91.1
## 12 2016 33.1 91.6
## 13 2017 33.1 91.6
## 14 2018 32.7 90.9
## 15 2019 31.2 88.2
## 16 2020 33.4 92.2
## 17 2021 33.2 91.8
ggplot(annual_temp, aes(x = year, y = mean_tmax_F)) +
geom_line() +
geom_point() +
labs(
title = "Annual Mean Daily Maximum Temperature",
subtitle = "Imperial County, California — Daymet location",
x = "Year",
y = "Temperature (°F)"
) +
theme_minimal()
library(leaflet)
leaflet() |>
addTiles() |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
)
imperial_sf
## Simple feature collection with 1 feature and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -116.1063 ymin: 32.61827 xmax: -114.4603 ymax: 33.43729
## Geodetic CRS: WGS 84
## GID_2 GID_0 COUNTRY GID_1 NAME_1 NL_NAME_1 NAME_2
## 1 USA.5.13_1 USA United States USA.5_1 California <NA> Imperial
## VARNAME_2 NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2 HASC_2
## 1 <NA> <NA> County County <NA> US.CA.IM
## geometry
## 1 POLYGON ((-114.7206 32.7179...
#install.packages("FedData")
library(FedData)
## Warning: package 'FedData' was built under R version 4.5.3
## You have loaded FedData v4.
## As of FedData v4 we have retired
## dependencies on the `sp` and `raster` packages.
## All functions in FedData v4 return `terra` (raster)
## or `sf` (vector) objects by default, and there may be
## other breaking changes.
library(leaflet)
library(htmlwidgets)
annual_temp$label <- paste0(
annual_temp$year,
": ",
round(annual_temp$mean_tmax_F, 1),
" °F"
)
m <- leaflet() |>
addTiles() |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
) |>
addCircleMarkers(
lng = -115.48,
lat = 33.04,
radius = 15,
popup = annual_temp$label[1]
)
m
library(leaflet)
# Create one map layer for each year
m <- leaflet() |>
addTiles() |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
)
for (i in seq_len(nrow(annual_temp))) {
m <- m |>
addCircleMarkers(
lng = -115.48,
lat = 33.04,
radius = 15,
layerId = paste0("year_", annual_temp$year[i]),
group = as.character(annual_temp$year[i]),
label = paste0(
annual_temp$year[i], ": ",
round(annual_temp$mean_tmax_F[i], 1), " °F"
),
popup = paste0(
"<b>Year:</b> ", annual_temp$year[i],
"<br><b>Mean daily maximum:</b> ",
round(annual_temp$mean_tmax_F[i], 1), " °F"
)
)
}
m |>
addLayersControl(
overlayGroups = as.character(annual_temp$year),
options = layersControlOptions(collapsed = FALSE)
)
library(leaflet)
# Create one map layer for each year
m <- leaflet() |>
addTiles() |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
)
for (i in seq_len(nrow(annual_temp))) {
m <- m |>
addCircleMarkers(
lng = -115.48,
lat = 33.04,
radius = 15,
layerId = paste0("year_", annual_temp$year[i]),
group = as.character(annual_temp$year[i]),
label = paste0(
annual_temp$year[i], ": ",
round(annual_temp$mean_tmax_F[i], 1), " °F"
),
popup = paste0(
"<b>Year:</b> ", annual_temp$year[i],
"<br><b>Mean daily maximum:</b> ",
round(annual_temp$mean_tmax_F[i], 1), " °F"
)
)
}
m |>
addLayersControl(
overlayGroups = as.character(annual_temp$year),
options = layersControlOptions(collapsed = FALSE)
)
library(leaflet)
library(htmlwidgets)
# Temperature values
years <- annual_temp$year
temps <- round(annual_temp$mean_tmax_F, 1)
# Start a fresh map
slider_map <- leaflet() |>
addTiles() |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
) |>
addCircleMarkers(
lng = -115.48,
lat = 33.04,
radius = 12,
color = "red",
fillColor = "red",
fillOpacity = 0.8
)
# Add slider control
slider_map <- htmlwidgets::onRender(
slider_map,
sprintf(
"
function(el, x) {
var years = [%s];
var temps = [%s];
var control = L.control({position: 'bottomleft'});
control.onAdd = function(map) {
var div = L.DomUtil.create('div', 'info');
div.style.background = 'white';
div.style.padding = '10px';
div.style.borderRadius = '6px';
div.innerHTML =
'<b>Year: <span id=\"selectedYear\">' + years[0] + '</span></b><br>' +
'<span id=\"selectedTemp\">' + temps[0] + ' °F</span><br>' +
'<input id=\"tempSlider\" type=\"range\" ' +
'min=\"0\" max=\"' + (years.length - 1) + '\" ' +
'value=\"0\" step=\"1\" style=\"width:250px\">';
return div;
};
control.addTo(map);
setTimeout(function() {
var slider = document.getElementById('tempSlider');
var yearText = document.getElementById('selectedYear');
var tempText = document.getElementById('selectedTemp');
slider.addEventListener('input', function() {
var i = Number(this.value);
yearText.innerHTML = years[i];
tempText.innerHTML = temps[i] + ' °F';
});
}, 200);
}
",
paste(years, collapse = ","),
paste(temps, collapse = ",")
)
)
slider_map
library(leaflet)
library(htmlwidgets)
years <- annual_temp$year
temps <- annual_temp$mean_tmax_F
# Temperature color function
temp_color <- function(t) {
if (t < 89) {
"blue"
} else if (t < 91) {
"orange"
} else {
"red"
}
}
slider_map <- leaflet() |>
addTiles() |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
) |>
addCircleMarkers(
lng = -115.48,
lat = 33.04,
radius = 20,
color = temp_color(temps[1]),
fillColor = temp_color(temps[1]),
fillOpacity = 0.7,
popup = paste0(
"<b>Year:</b> ", years[1],
"<br><b>Mean maximum temperature:</b> ",
round(temps[1], 1), " °F"
)
)
slider_map
#install.packages("leaflet.extras2")
library(leaflet.extras2)
## Warning: package 'leaflet.extras2' was built under R version 4.5.3
library(leaflet)
library(leaflet.extras2)
# Make a date for each annual observation
annual_temp$date <- as.Date(
paste0(annual_temp$year, "-01-01")
)
# Build the map
time_map <- leaflet() |>
addTiles() |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
) |>
addCircleMarkers(
lng = -115.48,
lat = 33.04,
radius = 18,
fillOpacity = 0.8,
color = "red",
fillColor = "red",
popup = paste0(
"<b>Year:</b> ", annual_temp$year,
"<br><b>Mean maximum temperature:</b> ",
round(annual_temp$mean_tmax_F, 1),
" °F"
)
)
time_map
library(terra)
# Create a temperature raster for Imperial County
# using the WorldClim maximum-temperature data we already downloaded
tmax_files <- list.files(
"data",
pattern = "tmax.*\\.tif$",
recursive = TRUE,
full.names = TRUE
)
tmax_files
## [1] "data/climate/wc2.1_10m/wc2.1_10m_tmax_01.tif"
## [2] "data/climate/wc2.1_10m/wc2.1_10m_tmax_02.tif"
## [3] "data/climate/wc2.1_10m/wc2.1_10m_tmax_03.tif"
## [4] "data/climate/wc2.1_10m/wc2.1_10m_tmax_04.tif"
## [5] "data/climate/wc2.1_10m/wc2.1_10m_tmax_05.tif"
## [6] "data/climate/wc2.1_10m/wc2.1_10m_tmax_06.tif"
## [7] "data/climate/wc2.1_10m/wc2.1_10m_tmax_07.tif"
## [8] "data/climate/wc2.1_10m/wc2.1_10m_tmax_08.tif"
## [9] "data/climate/wc2.1_10m/wc2.1_10m_tmax_09.tif"
## [10] "data/climate/wc2.1_10m/wc2.1_10m_tmax_10.tif"
## [11] "data/climate/wc2.1_10m/wc2.1_10m_tmax_11.tif"
## [12] "data/climate/wc2.1_10m/wc2.1_10m_tmax_12.tif"
library(terra)
tmax <- rast(tmax_files)
tmax
## class : SpatRaster
## size : 1080, 2160, 12 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : wc2.1_10m_tmax_01.tif
## wc2.1_10m_tmax_02.tif
## wc2.1_10m_tmax_03.tif
## ... and 9 more sources
## names : wc2.1~ax_01, wc2.1~ax_02, wc2.1~ax_03, wc2.1~ax_04, wc2.1~ax_05, wc2.1~ax_06, ...
## min values : -42.419, -39.58325, -53.400002, -59.548749, -59.8395, -60.360001, ...
## max values : 42.157, 40.2645, 41.48825, 43.17525, 44.8155, 46.6155, ...
# Make sure the county and raster use the same coordinate system
imperial_sf <- st_transform(imperial_sf, crs(tmax))
# Convert county boundary to a terra vector
imperial_vect <- vect(imperial_sf)
# Crop and mask the temperature raster to Imperial County
tmax_imperial <- crop(tmax, imperial_vect)
tmax_imperial <- mask(tmax_imperial, imperial_vect)
# Display the first month's temperature surface
plot(
tmax_imperial[[1]],
main = "Imperial County — Maximum Temperature"
)
lines(imperial_vect)
library(leaflet)
# Temperature range in our data
min_temp <- floor(min(annual_temp$mean_tmax_F, na.rm = TRUE))
max_temp <- ceiling(max(annual_temp$mean_tmax_F, na.rm = TRUE))
# Red/orange/yellow heat palette
heat_palette <- colorNumeric(
palette = c("#FFFFCC", "#FED976", "#FD8D3C", "#E31A1C", "#800026"),
domain = c(min_temp, max_temp)
)
heat_map <- leaflet() |>
addProviderTiles("CartoDB.Positron") |>
setView(
lng = -115.48,
lat = 33.04,
zoom = 8
) |>
# Temperature point
addCircleMarkers(
lng = -115.48,
lat = 33.04,
radius = 18,
color = "black",
weight = 2,
fillColor = heat_palette(annual_temp$mean_tmax_F[1]),
fillOpacity = 0.9,
popup = paste0(
"<b>Imperial County Daymet location</b><br>",
"<b>Year:</b> ", annual_temp$year[1], "<br>",
"<b>Mean daily maximum:</b> ",
round(annual_temp$mean_tmax_F[1], 1), " °F"
)
) |>
# Temperature legend
addLegend(
position = "bottomright",
pal = heat_palette,
values = annual_temp$mean_tmax_F,
title = "Mean Daily Maximum<br>Temperature (°F)",
opacity = 0.9
) |>
# Explanation box
addControl(
html = paste0(
"<div style='background:white; padding:12px; ",
"border-radius:6px; max-width:320px; ",
"box-shadow:0 1px 5px rgba(0,0,0,0.4);'>",
"<b>What am I looking at?</b><br>",
"This map shows the mean daily maximum temperature ",
"for the Daymet observation location in Imperial County.<br><br>",
"<b>Color scale:</b> yellow indicates cooler temperatures; ",
"orange and red indicate progressively hotter temperatures.<br><br>",
"<b>Important:</b> This point represents one Daymet grid location, ",
"not the temperature across the entire county.",
"</div>"
),
position = "topright"
)
heat_map
# Temperature anomaly relative to 2005
baseline <- annual_temp$mean_tmax_F[annual_temp$year == 2005]
annual_temp <- annual_temp |>
mutate(
anomaly_F = mean_tmax_F - baseline
)
ggplot(
annual_temp,
aes(x = year, y = anomaly_F)
) +
# Zero = 2005 baseline
geom_hline(
yintercept = 0,
linewidth = 1
) +
# Year-by-year anomaly
geom_col(
aes(fill = anomaly_F > 0)
) +
# Trend line
geom_smooth(
method = "lm",
se = TRUE,
color = "black",
linewidth = 1
) +
scale_fill_manual(
values = c(
"TRUE" = "#D73027",
"FALSE" = "#4575B4"
),
labels = c(
"TRUE" = "Warmer than 2005",
"FALSE" = "Cooler than 2005"
),
name = ""
) +
labs(
title = "Annual Temperature Anomaly",
subtitle = "Imperial County, California — Daymet observation location",
x = "Year",
y = "Difference from 2005 mean maximum temperature (°F)",
caption = "Positive values indicate warmer annual mean daily maximum temperatures than the 2005 baseline."
) +
theme_minimal() +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 11)
)
## `geom_smooth()` using formula = 'y ~ x'
