1 Preparation

dataFolder <- paste0(ws, "data/")

Load data from geopackage

charging_stationsB <- st_read(dsn=paste0(dataFolder, "chargingStations_berlin.shp"))
## Reading layer `chargingStations_berlin' from data source `/home/slautenb/Documents/lehre/HD/ws_2020_21/heikaLab/R4UrbanDataLab_2020/data/chargingStations_berlin.shp' using driver `ESRI Shapefile'
## Simple feature collection with 486 features and 20 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 373818.6 ymin: 5807999 xmax: 407135.9 ymax: 5829471
## projected CRS:  ETRS89 / UTM zone 33N

2 Explorative analysis

2.1 Kernel density map

We need to creat/define an empty raster that will store the results of the kernel density estimation

(bboxB <- st_bbox(charging_stationsB))
##      xmin      ymin      xmax      ymax 
##  373818.6 5807999.5  407135.9 5829470.9
gridsize=c(500,500)
range_x=list(c(bboxB[1]-1000, bboxB[3]+1000),c(bboxB[2]-1000, bboxB[4]+1000))
kdeMat <-  bkde2D(st_coordinates(charging_stationsB), bandwidth=c(500,500), gridsize=gridsize, range.x=range_x)
st_crs(charging_stationsB)
## Coordinate Reference System:
##   User input: ETRS89 / UTM zone 33N 
##   wkt:
## PROJCRS["ETRS89 / UTM zone 33N",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["UTM zone 33N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",15,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Europe - 12°E to 18°E and ETRS89 by country"],
##         BBOX[46.4,12,84.01,18.01]],
##     ID["EPSG",25833]]

Since bkde2D returns a matrix we need to convert it to a raster:

kdeRas <- raster(list(x=kdeMat$x1,y=kdeMat$x2,z=kdeMat$fhat), crs=CRS('+init=EPSG:25833') )
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum European_Terrestrial_Reference_System_1989 in CRS
## definition
kdeRas_scaled <- kdeRas*xres(kdeRas)*yres(kdeRas)* nrow(charging_stationsB)
cellStats(kdeRas_scaled$layer, stat=sum )
## [1] 485.8813
kdeRas_scaledContour <- rasterToContour(kdeRas_scaled)  %>%   st_as_sf()
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(kdeRas_scaled) + tm_raster(col="layer", alpha = .5, palette ="-plasma", title = "Kernel density estimate", style = "kmeans" ) + tm_shape(charging_stationsB) + tm_dots(alpha=.5) + tm_shape(kdeRas_scaledContour) + tm_lines() + tm_text(text="level", col="white", shadow=TRUE, along.lines=TRUE)

2.2 How far are the stations away from each other?

We use a for loop to calculate the distance to the clostet charging station in the data set. Therefore, we loop over all points and calculate the distance between the point i and all points including itself. This rather brute force approache returns the whole distance matrix which when can be used to calculate differences.

n <- nrow(charging_stationsB)
distMat <- matrix(data=0, nrow=n, ncol=n)
for(i in seq_len(n))
{
  allDists <- st_distance(charging_stationsB, charging_stationsB[i,])
  distMat[i,] <- allDists
  
}

In many cases several charging stations are at the same point, rendering a distance of zero.

getSecondSmallest <- function(x, pos=2)
{
   sort(x)[pos]
}
dist2Shortest <- apply(distMat, MARGI=1, FUN= getSecondSmallest, pos=2)
summary(dist2Shortest)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0   217.6   322.6   440.4  3392.7

We can force the analysis to deliver us the shortest distance to a neighbor not at the same location. Therefore, we: * sort values * get the indices of all values larger than zero * get the first element of the indices that is larger then zero

getSecondSmallestNotAtSame <- function(x, pos=1)
{
   sortedValues <- sort(x)
   idx <- which(sortedValues > 0)
   return(sortedValues[idx][pos])
}
dist2ShortestNotAtSame <- apply(distMat, MARGI=1, FUN= getSecondSmallestNotAtSame)
summary(dist2ShortestNotAtSame)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    1.303  185.636  348.441  470.288  583.351 3392.683

Or we could get the average distance to the closest ten charging stations

getMeanDistToClostest <- function(x, n=10)
{
   sortedValues <- sort(x)
   return(mean(sortedValues[1:n+1]))
}
distMean10 <- apply(distMat, MARGI=1, FUN= getMeanDistToClostest)
summary(distMean10)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    9.104  548.124  756.472 1009.831 1235.546 6542.951

If we want to show distances in the map we can simply add the vectors to the simple features object and use tmap afterwards.

charging_stationsB['meanDist10'] <- distMean10
osm_b <- read_osm(charging_stationsB, ext=1.1)
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(osm_b) + tm_rgb() + tm_shape(charging_stationsB) + tm_dots(size="meanDist10", legend.size.show = TRUE, title.size = "Mean distance to closest 10 stations" ) + tm_basemap(server= "OpenStreetMap.DE") + tm_layout(legend.outside=TRUE) + tm_scale_bar()

We see that distances are (not unexpectedly increaing from the city center).

A convenient and fast imlementations is kNNdist from the dbscan package:

nn4Mat <- kNNdist(st_coordinates(charging_stationsB), k=4, all=TRUE)
head(nn4Mat)
##          1         2        3        4
## 1   0.0000  94.39053 422.7904 427.1573
## 2   0.0000  94.39053 422.7904 427.1573
## 3   0.0000 266.16057 330.6558 766.6068
## 4 209.5521 565.63630 656.6839 852.1462
## 5 362.3375 565.63630 750.6605 808.8525
## 6 330.6558 330.65578 340.0982 586.6176

It returns the distance matrix for the k nearest neighbors.

2.3 Identify spatial clusters using dbscan

DBSCAN estimates the density around each data point by counting the number of points in a user-specified eps-neighborhood and applies a used-specified minPts thresholds to identify core, border and noise points. In a second step, core points are joined into a cluster if they are density-reachable (i.e., there is a chain of core points where one falls inside the eps-neighborhood of the next). Finally, border points are assigned to clusters. The algorithm only needs parameters eps and minPts.

db_B <- dbscan(st_coordinates(charging_stationsB), eps = 650, minPts = 7)
print(db_B)
## DBSCAN clustering for 486 objects.
## Parameters: eps = 650, minPts = 7
## The clustering contains 12 cluster(s) and 268 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9  10  11  12 
## 268  17  25  11   7  10  77  11   8  16  17   7  12 
## 
## Available fields: cluster, eps, minPts
charging_stationsB['dbCluster'] <- db_B$cluster
nClust <- length(unique(db_B$cluster))

Define our own palette based on the Set3 palette but uses white for the noise points.

mypal <- RColorBrewer::brewer.pal(nClust-1, "Set3")
mypal <- c("#FFFFFF", mypal)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(charging_stationsB) + 
  tm_dots(col = "dbCluster", size=0.25, alpha=.5, palette=mypal, n=nClust) +
  tm_basemap(server= c("OpenStreetMap.DE", "Esri.WorldImagery")) + 
  tm_scale_bar()