Mapping the Number of EV Charging Stations by County in Colorado Using R

EV
R
visualization
mapping
API
Author

Andy Pickering, PhD

Published

August 1, 2023

Modified

November 6, 2024

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.

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

Load Libraries
suppressPackageStartupMessages(library(tidyverse))
library(httr)
suppressPackageStartupMessages(library(jsonlite))
ggplot2::theme_set(theme_grey(base_size = 15))
library(leaflet)
suppressPackageStartupMessages(library(tigris)) # to get county shapefiles for maps
library(ggspatial) # for adding basemaps to ggplot2 maps
library(DT) # make nice data tables

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.

Load EV stations data from 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} (Walker 2023) package.

Code
options(tigris_use_cache = TRUE)
co_counties <- tigris::counties("CO",cb = TRUE, progress_bar = FALSE)
Retrieving data for the year 2022
Code
head(co_counties)
Simple feature collection with 6 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -109.0603 ymin: 36.99891 xmax: -104.3511 ymax: 39.91418
Geodetic CRS:  NAD83
   STATEFP COUNTYFP COUNTYNS       AFFGEOID GEOID        NAME
15      08      077 00198154 0500000US08077 08077        Mesa
16      08      083 00198157 0500000US08083 08083   Montezuma
17      08      067 00198148 0500000US08067 08067    La Plata
42      08      031 00198131 0500000US08031 08031      Denver
46      08      019 00198125 0500000US08019 08019 Clear Creek
47      08      055 00198143 0500000US08055 08055    Huerfano
             NAMELSAD STUSPS STATE_NAME LSAD      ALAND   AWATER
15        Mesa County     CO   Colorado   06 8621348059 31991710
16   Montezuma County     CO   Colorado   06 5255990019 27208195
17    La Plata County     CO   Colorado   06 4376255278 25642579
42      Denver County     CO   Colorado   06  396460127  4275563
46 Clear Creek County     CO   Colorado   06 1023234059  3274900
47    Huerfano County     CO   Colorado   06 4120756290  5792101
                         geometry
15 MULTIPOLYGON (((-109.0603 3...
16 MULTIPOLYGON (((-109.0459 3...
17 MULTIPOLYGON (((-108.3796 3...
42 MULTIPOLYGON (((-104.9341 3...
46 MULTIPOLYGON (((-105.927 39...
47 MULTIPOLYGON (((-105.5013 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 (Table 1).

Code
zips <- readr::read_csv("data/zip_code_database.csv",
                        show_col_types = FALSE) %>% 
  filter(state == "CO") %>% 
  select(zip, primary_city, county)

zips |>
  DT::datatable(options = list(pageLength = 5), rownames = FALSE)
Table 1: Zip code data from USPS

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 (Table 2).

Code
ev_county_counts <- ev %>% 
  select(id,zip,city) %>% 
  left_join(zips, by = "zip") %>% 
  dplyr::count(county) %>% 
  arrange(desc(n))

ev_county_counts |>
  DT::datatable(options = list(pageLength = 5), rownames = FALSE)
Table 2: Number of EV charging stations per county

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.

Code
co_ev_counts <- co_counties %>% 
  left_join(ev_county_counts, by = c("NAMELSAD" = "county"))

co_ev_counts <- sf::st_transform(co_ev_counts, 4326)

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.99891 xmax: -104.3511 ymax: 39.91418
Geodetic CRS:  WGS 84
  STATEFP COUNTYFP COUNTYNS       AFFGEOID GEOID        NAME           NAMELSAD
1      08      077 00198154 0500000US08077 08077        Mesa        Mesa County
2      08      083 00198157 0500000US08083 08083   Montezuma   Montezuma County
3      08      067 00198148 0500000US08067 08067    La Plata    La Plata County
4      08      031 00198131 0500000US08031 08031      Denver      Denver County
5      08      019 00198125 0500000US08019 08019 Clear Creek Clear Creek County
6      08      055 00198143 0500000US08055 08055    Huerfano    Huerfano County
  STUSPS STATE_NAME LSAD      ALAND   AWATER   n                       geometry
1     CO   Colorado   06 8621348059 31991710  61 MULTIPOLYGON (((-109.0603 3...
2     CO   Colorado   06 5255990019 27208195   5 MULTIPOLYGON (((-109.0459 3...
3     CO   Colorado   06 4376255278 25642579  39 MULTIPOLYGON (((-108.3796 3...
4     CO   Colorado   06  396460127  4275563 356 MULTIPOLYGON (((-104.9341 3...
5     CO   Colorado   06 1023234059  3274900  13 MULTIPOLYGON (((-105.927 39...
6     CO   Colorado   06 4120756290  5792101   4 MULTIPOLYGON (((-105.5013 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.

ggplot

  • {ggplot2} makes it relatively easy to plot spatial data in an sf object with the geom_sf function
  • The {scales} (Wickham and Seidel 2022) package is used to format the numbers in the legend
  • The {ggspatial} (Dunnington 2023) package is used to add a base map showing some of the major cities and roads
Code
ggplot() +
  ggspatial::annotation_map_tile(progress = "none") +
  geom_sf(data = co_ev_counts,
          aes(fill = n),
          alpha = 0.5) +
  scale_fill_viridis_c(labels = scales::number_format(big.mark = ","),
                       name = '# Ev Stations') +
  ggtitle("Number of EV Stations by Colorado County") +
  theme_void()
Figure 1: Choropleth map of number of EV charging stations by county, made with ggplot2

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 (Figure 2) 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.

Code
# 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"
            )

# 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)
Figure 2: Interactive choropleth map of number of EV charging stations by county. Click on a county polygon to display information.

Future Analysis

Now that I have the basic framework set up for mapping the EV data, there are a lot of other interesting questions I would like to investigate.

  • Look at breakdown by charger level/type/network etc..

  • Look at trends over time.

  • Look at relationship between demographics (population, income etc. ) and chargers. The tidycensus package would probably be useful for this.

  • Extend to other states or similar analysis at state level.

SessionInfo

In order to improve the reproducibility of this analysis, I include the sessionInfo below at the time this post was rendered.

SessionInfo

R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Denver
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] DT_0.33         ggspatial_1.1.9 tigris_2.1      leaflet_2.2.2  
 [5] jsonlite_1.8.8  httr_1.4.7      lubridate_1.9.3 forcats_1.0.0  
 [9] stringr_1.5.1   dplyr_1.1.4     purrr_1.0.2     readr_2.1.5    
[13] tidyr_1.3.1     tibble_3.2.1    ggplot2_3.5.1   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.5       bslib_0.8.0        xfun_0.46          htmlwidgets_1.6.4 
 [5] tzdb_0.4.0         vctrs_0.6.5        tools_4.4.1        crosstalk_1.2.1   
 [9] generics_0.1.3     curl_5.2.1         parallel_4.4.1     proxy_0.4-27      
[13] fansi_1.0.6        pkgconfig_2.0.3    KernSmooth_2.23-24 RColorBrewer_1.1-3
[17] uuid_1.2-0         lifecycle_1.0.4    farver_2.1.1       compiler_4.4.1    
[21] munsell_0.5.1      sass_0.4.9         htmltools_0.5.8.1  class_7.3-22      
[25] yaml_2.3.10        jquerylib_0.1.4    pillar_1.9.0       crayon_1.5.3      
[29] classInt_0.4-10    cachem_1.1.0       tidyselect_1.2.1   digest_0.6.36     
[33] rosm_0.3.0         stringi_1.8.4      sf_1.0-17          labeling_0.4.3    
[37] fastmap_1.2.0      grid_4.4.1         colorspace_2.1-0   cli_3.6.3         
[41] prettymapr_0.2.5   magrittr_2.0.3     utf8_1.2.4         e1071_1.7-14      
[45] withr_3.0.1        scales_1.3.0       rappdirs_0.3.3     bit64_4.0.5       
[49] timechange_0.3.0   rmarkdown_2.27     bit_4.0.5          png_0.1-8         
[53] hms_1.1.3          evaluate_0.24.0    knitr_1.48         viridisLite_0.4.2 
[57] rlang_1.1.4        Rcpp_1.0.13        glue_1.7.0         DBI_1.2.2         
[61] renv_1.0.4         rstudioapi_0.16.0  vroom_1.6.5        plyr_1.8.9        
[65] R6_2.5.1           units_0.8-5

References

Cheng, Joe, Bhaskar Karambelkar, and Yihui Xie. 2023. “Leaflet: Create Interactive Web Maps with the JavaScript ’Leaflet’ Library.” https://CRAN.R-project.org/package=leaflet.
Dunnington, Dewey. 2023. “Ggspatial: Spatial Data Framework for Ggplot2.” https://CRAN.R-project.org/package=ggspatial.
Walker, Kyle. 2023. “Tigris: Load Census TIGER/Line Shapefiles.” https://CRAN.R-project.org/package=tigris.
Wickham, Hadley. 2016. “Ggplot2: Elegant Graphics for Data Analysis.” https://ggplot2.tidyverse.org.
Wickham, Hadley, and Dana Seidel. 2022. “Scales: Scale Functions for Visualization.” https://CRAN.R-project.org/package=scales.