Loading the CSV file for Electric Charging Stations Database from the Department of Energy. Hidden because it includes my API key. You can get the free API key at https://developer.nrel.gov/signup/ The code looks something like this, exchanging the YOUR_KEY_HERE with your key:

myurl<-“https://developer.nrel.gov/api/alt-fuel-stations/v1.csv?fuel_type=ELEC&api_key=YOUR_KEY_HERE&format=CSV” efuelstations <- read.csv(file=url(myurl), header=TRUE)

Showing the column names:

colnames(efuelstations)
##  [1] "Fuel.Type.Code"          "Station.Name"           
##  [3] "Street.Address"          "Intersection.Directions"
##  [5] "City"                    "State"                  
##  [7] "ZIP"                     "Plus4"                  
##  [9] "Station.Phone"           "Status.Code"            
## [11] "Expected.Date"           "Groups.With.Access.Code"
## [13] "Access.Days.Time"        "Cards.Accepted"         
## [15] "BD.Blends"               "NG.Fill.Type.Code"      
## [17] "NG.PSI"                  "EV.Level1.EVSE.Num"     
## [19] "EV.Level2.EVSE.Num"      "EV.DC.Fast.Count"       
## [21] "EV.Other.Info"           "EV.Network"             
## [23] "EV.Network.Web"          "Geocode.Status"         
## [25] "Latitude"                "Longitude"              
## [27] "Date.Last.Confirmed"     "ID"                     
## [29] "Updated.At"              "Owner.Type.Code"        
## [31] "Federal.Agency.ID"       "Federal.Agency.Name"    
## [33] "Open.Date"               "Hydrogen.Status.Link"   
## [35] "NG.Vehicle.Class"        "LPG.Primary"            
## [37] "E85.Blender.Pump"        "EV.Connector.Types"     
## [39] "Country"                 "Hydrogen.Is.Retail"     
## [41] "Access.Code"             "Access.Detail.Code"     
## [43] "Federal.Agency.Code"

We see the columns Longitude and Latitude - the key columns to project the stations onto a map. There is also Open.Date as one of the key variables and Groups.With.Access.Code which distinguishes between private and public stations (e.g. only for city works).

Now, let’s put this on a map.

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library(RQGIS)
## Loading required package: reticulate
library(mapview)
usco<-st_read(dsn="E:/Data/USCounties", layer = "tl_2015_us_county")
## Reading layer `tl_2015_us_county' from data source `E:\Data\USCounties' using driver `ESRI Shapefile'
## Simple feature collection with 3233 features and 17 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.44106
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
us_proj<-st_transform(usco,crs=102003)
#1020003 is the EPSG code for our Albers projection

plot(us_proj["STATEFP"])

For the poster, I am using census tracts. Unfortunally, there are some inconsistencies with Guam and Virgin Islands, etc. (has no EV charging stations listed with the Dept. of Energy) and Alaska (the Census bereau split one tract on the map to both ends, making a centroid imposible to be reliable - it’s over Eastern Canada) Hence, I will ignore everything but the continental 48 states on the map since I work with tracts and centroids for distances.

#Loading the shapefile for the US Census Tracts oft he 48 continental states, no Alaska. I created the file in QGIS by downloading all state shapefiles from the Census Burea TIGER ftp and merging them into one file. Then, I had to download a shapefile that excludes all coastal zones and the Great Lakes, cookiecutting the combined file down to land masses. Took a while but it worked just fine. The now trimmed layer is the one below. If I would not have done that, then, the centroids of the census tracts around the lakes and on the coasts would have been in the water - making the algorithm useless.
us48tracts<-st_read(dsn="E:/Data/USTracts", layer = "us48tracts")
## Reading layer `us48tracts' from data source `E:\Data\USTracts' using driver `ESRI Shapefile'
## Simple feature collection with 72373 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -124.7631 ymin: 24.5231 xmax: -66.9499 ymax: 49.38436
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
#Loading the Hub Distance. Hub distance is the centroid-to-next-station distance from the census tract centroids.
hubdistance<-st_read(dsn="E:/Data/USTracts", layer="us48HubDistance")
## Reading layer `us48HubDistance' from data source `E:\Data\USTracts' using driver `ESRI Shapefile'
## Simple feature collection with 72373 features and 14 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -124.6281 ymin: 24.54923 xmax: -67.01832 ymax: 48.98816
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
us_proj<-st_transform(us48tracts,crs=5070)
#1020003 is the EPSG code for our Albers projection

plot(us_proj["STATEFP"])

testmap2<-st_read(dsn="E:/Data/USTracts", layer="us48LowStates")
## Reading layer `us48LowStates' from data source `E:\Data\USTracts' using driver `ESRI Shapefile'
## Simple feature collection with 49 features and 14 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -124.849 ymin: 24.39631 xmax: -66.88544 ymax: 49.38436
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
testmap3<-st_read(dsn="E:/Data/USTracts", layer="us48EVCS")
## Reading layer `us48EVCS' from data source `E:\Data\USTracts' using driver `ESRI Shapefile'
## Simple feature collection with 18106 features and 45 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -124.6631 ymin: 24.55055 xmax: -67.39852 ymax: 48.99049
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
#Now, I have to color the census tracts with with a Jenks break color layer according to the Hub Distance.


library(classInt)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
testmap<-us48tracts %>%
  mutate(Tract_to_EV_distance=cut(hubdistance$HubDist,breaks=data.frame(classIntervals(var=hubdistance$HubDist, n=7, style="jenks")[2])[,1], include.lowest = T),include.lowest = T)
library(ggsn)
## Loading required package: ggplot2
## Loading required package: grid
p1<-ggplot(testmap, aes(fill = Tract_to_EV_distance, color = Tract_to_EV_distance)) + 
  geom_sf() + 
  geom_sf(data = testmap2 , fill = NA , color = "black") +
  #geom_sf(data = testmap3, ) + #commented out because the would overlay the smaller census tracts
  ggtitle("Distance Census Tract Centroid to Closest EV Charging Station", 
          subtitle = "USA, Lower 48 States, June 2018, Census Bureau, US Dept. of Energy")+
    scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
  theme(plot.title = element_text(size=30,face = "bold"),plot.subtitle = element_text(size=28),axis.text.x = element_blank(), axis.text.y = element_blank(),legend.title = element_text(size=28,face="bold"),legend.text = element_text(size=28),legend.position = c(0.9, 0.3))+
  north(testmap)+
  scalebar(testmap, dist = 200,  dd2km =T, model="GRS80", st.size = 8)
ggsave("p1.png",width=30,height=18,units="in")
p1

#Writing the files into local folders
sf::st_write(testmap,dsn="E:/Data/USTracts",layer="us48TractDistance", driver="ESRI Shapefile", delete_layer=T, update=T)
## Deleting layer `us48TractDistance' using driver `ESRI Shapefile'
## Updating layer `us48TractDistance' to data source `E:\Data\USTracts' using driver `ESRI Shapefile'
## features:       72373
## fields:         15
## geometry type:  Multi Polygon

This one is a good way to fill up your RAM until R crashes - mine crashed after occupying 11GB. So, this is a headache for another day. #{r} #library(mapview) #library(RColorBrewer) #pal <- colorRampPalette(brewer.pal(6, "Blues")) #set colors #mapview(testmap["Tract_to_EV_distance"], col.regions=pal, legend=T,map.types="OpenStreetMap", #layer.name="Tract Centroid Distances to Next EVCS") #

Now, let’s try to explore the match/mismatch of area of a polygon and Hub Distance as a inverted correlation: The area for the census tracts is ALAND in square meters. One square kilometer is (ALAND/1000000) since a square kilometer is 1000x1000 meters. The population is somewhat even in census tracts except for really lowly populated areas. us48HubDistance was created in QGIS by creating centroids and calculating the HubDist (distance from centroid to next EV station) and it has the ALAND variable as well inherited from the calcuation the layer comes from. Et voila! Let’s plot them!

hubdistance2<-st_read(dsn="E:/Data/USTracts", layer="us48DistCalc")
## Reading layer `us48DistCalc' from data source `E:\Data\USTracts' using driver `ESRI Shapefile'
## Simple feature collection with 72373 features and 15 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -124.6281 ymin: 24.54923 xmax: -67.01832 ymax: 48.98816
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
p2<-ggplot(hubdistance2, aes(x = HubDist, y = Area)) +
  geom_point(size=4,shape=16) + 
  ggtitle("Scatterplot EV Charging Station Distance Relation to Census Tract Area", 
          subtitle = "USA, Lower 48 States, June 2018, Census Bureau, US Dept. of Energy") +
  theme(plot.title = element_text(size=30,face = "bold"),plot.subtitle = element_text(size=28),axis.text.x = element_text(size = 28), axis.text.y = element_text(size = 28),axis.title.x = element_text(size=28),axis.title.y = element_text(size=28)) +
  labs(x="Distance between census tract centroid and next EVCS") +
  labs(y="Area of census tract")
ggsave("p2.png",width=20,height=20,units="in")
## Warning: Removed 819 rows containing missing values (geom_point).
p2
## Warning: Removed 819 rows containing missing values (geom_point).

#It looks correct in the file. Mind that the large font etc. is for a 20x20 inches output file.