Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Introduction
As you may know from a previous post I am interested in electric-vehicle (EV) trends and the transition to a more electrified transportation fleet. I wanted to do some mapping and spatial analysis, and I recently took the Creating Maps with R course by Charlie Joey Hadley, so I decided to use some of the skills I learned to create some maps of EV charging station data for Colorado.
< section id="goal" class="level2">Goal
- My goal in this post is to create choropleth map(s) showing the number of EV charging stations per county in Colorado.
Data & Analysis
suppressPackageStartupMessages(library(tidyverse)) library(httr) suppressPackageStartupMessages(library(jsonlite)) ggplot2::theme_set(theme_grey(base_size = 15)) library(leaflet) suppressPackageStartupMessages(library(tigris))
EV Stations data
Data on EV stations is obtained from the Alternative Fuels Data Center’s Alternative Fuel Stations database. See my previous post for more details on getting the data from the API.
# API key is stored in my .Renviron file
api_key <- Sys.getenv("AFDC_KEY")
target <- "https://developer.nrel.gov/api/alt-fuel-stations/v1"
# Return data for all electric stations in Colorado
api_path <- ".json?&fuel_type=ELEC&state=CO&limit=all"
complete_api_path <- paste0(target,api_path,'&api_key=',api_key)
response <- httr::GET(url = complete_api_path)
if (response$status_code != 200) {
print(paste('Warning, API call returned error code', response$status_code))
}
ev_dat <- jsonlite::fromJSON(httr::content(response,"text"))
ev <- ev_dat$fuel_stations
# filter out non-EV related fields
ev <- ev %>% select(-dplyr::starts_with("lng")) %>%
select(-starts_with("cng")) %>%
select(-starts_with("lpg")) %>%
select(-starts_with("hy")) %>%
select(-starts_with("ng")) %>%
select(-starts_with("e85")) %>%
select(-starts_with("bd")) %>%
select(-starts_with("rd")) %>%
filter(status_code == 'E')
County data
Next I need shape files for the Colorado counties to make the map; these are obtained from the tigris package.
options(tigris_use_cache = TRUE)
co_counties <- tigris::counties("CO",cb = TRUE)
Retrieving data for the year 2021
head(co_counties)
Simple feature collection with 6 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -109.0603 ymin: 36.99902 xmax: -104.6606 ymax: 39.92525
Geodetic CRS: NAD83
STATEFP COUNTYFP COUNTYNS AFFGEOID GEOID NAME NAMELSAD
7 08 037 00198134 0500000US08037 08037 Eagle Eagle County
150 08 059 00198145 0500000US08059 08059 Jefferson Jefferson County
151 08 067 00198148 0500000US08067 08067 La Plata La Plata County
165 08 077 00198154 0500000US08077 08077 Mesa Mesa County
174 08 035 00198133 0500000US08035 08035 Douglas Douglas County
209 08 043 00198137 0500000US08043 08043 Fremont Fremont County
STUSPS STATE_NAME LSAD ALAND AWATER geometry
7 CO Colorado 06 4362754228 18970639 MULTIPOLYGON (((-107.1137 3...
150 CO Colorado 06 1979735379 25071495 MULTIPOLYGON (((-105.0558 3...
151 CO Colorado 06 4376255277 25642579 MULTIPOLYGON (((-108.3796 3...
165 CO Colorado 06 8621348071 31991710 MULTIPOLYGON (((-109.0603 3...
174 CO Colorado 06 2176020279 6795841 MULTIPOLYGON (((-105.3274 3...
209 CO Colorado 06 3972613670 2235542 MULTIPOLYGON (((-106.0123 3...
Zip codes
I have the EV station data and the county shape files, so the next step is to join them together. However, I have a problem: the EV stations data does not contain the county name or code, so I can’t join them yet without a common column. There are probably a lot of different solutions to this problem (for example the EV data contains addresses so I could geo-code these to get the county). In this case, I decided the easiest solution was to download the zip code database from the USPS (free for personal use), which contains both zip codes and their corresponding county.
zips <- readr::read_csv("data/zip_code_database.csv",
show_col_types = FALSE) %>%
filter(state == "CO") %>%
select(zip, primary_city, county)
head(zips)
# A tibble: 6 × 3 zip primary_city county <chr> <chr> <chr> 1 80001 Arvada Jefferson County 2 80002 Arvada Jefferson County 3 80003 Arvada Jefferson County 4 80004 Arvada Jefferson County 5 80005 Arvada Jefferson County 6 80006 Arvada Jefferson County
Next I compute the number of stations per zip code in the EV data, and join to the zip code database to add the county column.
ev_county_counts <- ev %>% select(id,zip,city) %>% left_join(zips, by = "zip") %>% dplyr::count(county) %>% arrange(desc(n)) head(ev_county_counts)
county n 1 Denver County 320 2 Boulder County 303 3 Larimer County 145 4 Jefferson County 142 5 Arapahoe County 118 6 Eagle County 103
Combining data
Now we can finally join the data we want to plot (# EV stations per county) in ev_county_counts to our sf object (co_counties) with the county spatial data, and we are ready to make some maps.
co_ev_counts <- co_counties %>%
left_join(ev_county_counts, by = c("NAMELSAD" = "county"))
head(co_ev_counts)
Simple feature collection with 6 features and 13 fields Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: -109.0603 ymin: 36.99902 xmax: -104.6606 ymax: 39.92525 Geodetic CRS: NAD83 STATEFP COUNTYFP COUNTYNS AFFGEOID GEOID NAME NAMELSAD 1 08 037 00198134 0500000US08037 08037 Eagle Eagle County 2 08 059 00198145 0500000US08059 08059 Jefferson Jefferson County 3 08 067 00198148 0500000US08067 08067 La Plata La Plata County 4 08 077 00198154 0500000US08077 08077 Mesa Mesa County 5 08 035 00198133 0500000US08035 08035 Douglas Douglas County 6 08 043 00198137 0500000US08043 08043 Fremont Fremont County STUSPS STATE_NAME LSAD ALAND AWATER n geometry 1 CO Colorado 06 4362754228 18970639 103 MULTIPOLYGON (((-107.1137 3... 2 CO Colorado 06 1979735379 25071495 142 MULTIPOLYGON (((-105.0558 3... 3 CO Colorado 06 4376255277 25642579 30 MULTIPOLYGON (((-108.3796 3... 4 CO Colorado 06 8621348071 31991710 47 MULTIPOLYGON (((-109.0603 3... 5 CO Colorado 06 2176020279 6795841 57 MULTIPOLYGON (((-105.3274 3... 6 CO Colorado 06 3972613670 2235542 6 MULTIPOLYGON (((-106.0123 3...
Mapping
I’m going to make choropleth maps using two of the more popular mapping packages: ggplot2 (Wickham 2016) and leaflet (Cheng, Karambelkar, and Xie 2023). I think they both make good-looking maps; the main advantage to leaflet is that the map is interactive.
< section id="ggplot" class="level2">ggplot
ggplot2 makes it relatively easy to plot spatial data in an sf object with geom_sf
I’m also use the scales (Wickham and Seidel 2022) package to format the numbers in the legend
ggplot() +
geom_sf(data = co_ev_counts,
aes(fill = n)) +
scale_fill_viridis_c(labels = scales::number_format(big.mark = ","),
name = '# Ev Stations') +
ggtitle("Number of EV Stations by Colorado County") +
theme_void()
Leaflet
Using leaflet requires a little more code but allows you to create an interactive map that can be more useful to the reader.
In the map below I’ve set the popup to display the county name and number of stations when you click on the map.
You can also drag the map around and zoom in/out.
It’s also very easy with Leaflet to add a basemap (OpenStreetMap in this case) layer under the choropleth. I decided to add this here to give readers a better sense of context, and also because I wanted to highlight that the counties close to major highways (I-70 east-west and I-25 north-south) appear to have higher numbers of chargers.
Note I’ve also included some code using from the Creating Maps in R course to fix an issue in the legend where the NA entry overlaps with the other entries.
# create color palette
pal_ev <- leaflet::colorNumeric(palette = "viridis",
domain = co_ev_counts$n)
co_ev_map <- leaflet() %>%
addTiles() %>% # adds OpenStretMap basemap
addPolygons(data = co_ev_counts,
weight = 1,
color = "black",
popup = paste(co_ev_counts$NAME, "<br>",
" EV Stations: ", co_ev_counts$n, "<br>"),
fillColor = ~pal_ev(n),
fillOpacity = 0.6) %>%
addLegend(data = co_ev_counts,
pal = pal_ev,
values = ~n,
opacity = 1,
title = "# of EV Stations <br>
Per County"
)
Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs). Need '+proj=longlat +datum=WGS84'
# legend fix --------------------------------------------------------------
# for issue with na in legend
html_fix <- htmltools::tags$style(type = "text/css", "div.info.legend.leaflet-control br {clear: both;}")
co_ev_map %>%
htmlwidgets::prependContent(html_fix)
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.