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

EV
R
visualization
mapping
API
Author

Andy Pickering

Published

August 1, 2023

Modified

May 13, 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 2021
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.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 (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.99902 xmax: -104.6606 ymax: 39.92525
Geodetic CRS:  WGS 84
  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 106 MULTIPOLYGON (((-107.1137 3...
2     CO   Colorado   06 1979735379 25071495 172 MULTIPOLYGON (((-105.0558 3...
3     CO   Colorado   06 4376255277 25642579  39 MULTIPOLYGON (((-108.3796 3...
4     CO   Colorado   06 8621348071 31991710  53 MULTIPOLYGON (((-109.0603 3...
5     CO   Colorado   06 2176020279  6795841  62 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.

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

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.

Code 
R version 4.3.1 (2023-06-16)
Platform: x86_64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.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.30         ggspatial_1.1.9 tigris_2.0.4    leaflet_2.2.0  
 [5] jsonlite_1.8.7  httr_1.4.7      lubridate_1.9.3 forcats_1.0.0  
 [9] stringr_1.5.0   dplyr_1.1.3     purrr_1.0.2     readr_2.1.4    
[13] tidyr_1.3.0     tibble_3.2.1    ggplot2_3.4.4   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.4       xfun_0.40          bslib_0.5.1        htmlwidgets_1.6.2 
 [5] tzdb_0.4.0         vctrs_0.6.4        tools_4.3.1        crosstalk_1.2.0   
 [9] generics_0.1.3     curl_5.1.0         parallel_4.3.1     proxy_0.4-27      
[13] fansi_1.0.5        pkgconfig_2.0.3    KernSmooth_2.23-22 RColorBrewer_1.1-3
[17] uuid_1.1-1         lifecycle_1.0.3    farver_2.1.1       compiler_4.3.1    
[21] munsell_0.5.0      sass_0.4.7         htmltools_0.5.6.1  class_7.3-22      
[25] yaml_2.3.7         jquerylib_0.1.4    pillar_1.9.0       crayon_1.5.2      
[29] ellipsis_0.3.2     classInt_0.4-10    cachem_1.0.8       viridis_0.6.4     
[33] tidyselect_1.2.0   digest_0.6.33      rosm_0.3.0         stringi_1.7.12    
[37] sf_1.0-14          labeling_0.4.3     fastmap_1.1.1      grid_4.3.1        
[41] colorspace_2.1-0   cli_3.6.1          prettymapr_0.2.4   magrittr_2.0.3    
[45] utf8_1.2.4         e1071_1.7-13       withr_2.5.1        scales_1.2.1      
[49] rappdirs_0.3.3     bit64_4.0.5        timechange_0.2.0   rmarkdown_2.25    
[53] bit_4.0.5          gridExtra_2.3      png_0.1-8          hms_1.1.3         
[57] evaluate_0.22      knitr_1.44         viridisLite_0.4.2  rlang_1.1.1       
[61] Rcpp_1.0.11        glue_1.6.2         DBI_1.1.3          renv_1.0.3        
[65] rstudioapi_0.15.0  vroom_1.6.4        plyr_1.8.9         R6_2.5.1          
[69] units_0.8-4

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.