1. Research Objectives and Scope

This research aims to identify building connectivity based on hourly energy demands and potential supply by solar photovoltaic, cluster energy sharing groups of buildings, and reconfigure block boundaries which can be based on multiple microgrids. The research applied into a community located in Kyojima, Tokyo, Japan.

2. Preparation of R packages, Set working directory, and set dataframe

library(data.table)
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(mapview)
library(leaflet)
library(ggplot2)
library(ggvoronoi)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(spdep)
## Loading required package: sp
## 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(FNN)
library(ggspatial)
library(scales)
library(ggpubr)
## Loading required package: magrittr
library(igraph)
## 
## Attaching package: 'igraph'
## The following object is masked from 'package:FNN':
## 
##     knn
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
# set working directory
setwd("C:/Users/Soowon/Desktop/Data_SpatialClustering")

# read shapefile
buildings2 <- st_read("CCM22_AllBuildings_4.shp")
## Reading layer `CCM22_AllBuildings_4' from data source `C:\Users\Soowon\Desktop\Data_SpatialClustering\CCM22_AllBuildings_4.shp' using driver `ESRI Shapefile'
## Simple feature collection with 870 features and 106 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -1488.226 ymin: -31730.25 xmax: -933.4516 ymax: -31052.29
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=36 +lon_0=139.8333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
CCM <- st_read("Kyojima1_boundary.shp")
## Reading layer `Kyojima1_boundary' from data source `C:\Users\Soowon\Desktop\Data_SpatialClustering\Kyojima1_boundary.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 23 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 139.8167 ymin: 35.71393 xmax: 139.8233 ymax: 35.72027
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
st_crs(buildings2) <- 6677 # JGD2011_9 (meter scale)
buildings2_map <- st_transform(buildings2, crs=6668) # JGD2011 (degree scale)

3. Mapping existing blocks

mapview(buildings2, zcol="Num_model", map.types="CartoDB.Positron")
ggplot()+
  geom_sf(data=buildings2,aes(fill=as.factor(Num_model)))+
  theme_bw()

4. Calaulationg parts

# dataframe: just only x and y coordinates
coord <- st_set_geometry(buildings2[,c("longm","latm")],NULL) # meter scale
coords <- st_set_geometry(buildings2[,c("long","lat")],NULL) # degree scale

# number of buildings
n <- dim(buildings2)[1]

# hourly building PV Supply simulation data
Supply <- st_set_geometry(buildings2[,c(59:82)],NULL)
# hourly building Energy Demand simulation data
Demand <- st_set_geometry(buildings2[,c(83:106)],NULL) 

# Summarize supply and demand data in dataframe manner
sup <- Supply # dim(sup): 870rows * 24columns
dem <- Demand # dim(dem): 870rows * 24columns

#calculate self-efficiency
sup[sup<0] <- 0
dem[dem<0] <- 0     
pp <- as.matrix(sup-dem) #surplus +, shortage -
ng <- as.matrix(dem-sup) #shortage +, surplus -
self1   <-ifelse(pp>0, pp, 0)   #surplus +, shortage 0
loss1   <-ifelse(ng>0, ng, 0) #shortage +, surplus 0

for(k in seq(2,100,1)){     
  knnn <- get.knn(coord,k)
  knid <- cbind(1:n,knnn$nn.index)      
  Sup1 <- rowMeans(matrix(sup[knid],ncol=k+1)) #total supply 
  Dem1 <- rowMeans(matrix(dem[knid],ncol=k+1)) #total demand        
  PP <- Sup1-Dem1
  NG <- Dem1-Sup1
  Self1 <- ifelse(PP>0, PP, 0) #supply surplus (Max is best)    
  self1 <- cbind(self1, Self1)      
  Loss1 <- ifelse(NG>0, NG, 0) #demand surplus (Min is best)
  loss1 <- cbind(loss1, Loss1)      
}
# supply surplus
selfMax <- apply(self1, 1, max)
selfSel <- ifelse(self1==selfMax,1,0)
for(j in 1:length(selfSel[1,])) selfSel[,j]<-selfSel[,j]*j
selfSel[selfSel==0]<- 1000
self_sel <- apply(selfSel, 1, min)
self_val <- rep(0, n)
for(jj in 1:n) self_val[jj] <- self1[jj,self_sel[jj]]

# demand surplus 
lossMin <- apply(loss1, 1, min)
lossSel <- ifelse(loss1==lossMin,1,0)
for(j in 1:length(lossSel[1,])) lossSel[,j]<-lossSel[,j]*j
lossSel[lossSel==0] <- 1000
loss_sel <- apply(lossSel, 1, min)
loss_val <- rep(0, n)
for(jj in 1:n) loss_val[jj]   <- loss1[jj,loss_sel[jj]]

# arrenge objects
Self1<-NULL; Self_sel<-NULL; Loss1<-NULL; Loss_sel<-NULL
self <- Self <- cbind(Self1,self_val)
self_sel <- Self_sel <- cbind(Self_sel,self_sel)
loss <- Loss <- cbind(Loss1,loss_val)
loss_sel <- Loss_sel <- cbind(Loss_sel,loss_sel)

# calculate neibours and its distance
dis <-rep(0, n)     
for(k in seq(1,100,1)){     
  knnn    <-get.knn(coord,k)    
  dis     <-cbind(dis, rowSums(knnn$nn.dist))   
}
#
link <- c(1,seq(1,100,1))
n   <- length(dis[,1])
nl <- length(link)
#
self2    <-matrix(0, ncol=nl,nrow=n)
self2_sel<-matrix(0, ncol=nl,nrow=n)

#
for(i in 1:n){
  data_sub <-data.frame(table(as.matrix(self_sel[i,])))
  data_sub2<-data_sub[order(data_sub[,2],decreasing=T),]
  ids   <-as.numeric(as.character(data_sub2[,1]))
  av_sel<-rep(0, nl)
  num   <-1
  for(id in ids){
    av_sel[num]<-sum(self[i,self_sel[i,]==id])
    num <-num + 1
  }
  self2[i,]<-av_sel
  self2_sel[i,1:length(ids)]<-ids
}

5A. Visualization - Building Connectivity

#vis in degree scale
coord2 <- st_set_geometry(buildings2_map[,c("long","lat")],NULL)
coord2_sf <- st_as_sf(coord2, coords = c("long", "lat"), crs=6668)

# calculate the dataframe of neibours sharing energy
A <- NULL
knns <- get.knn(coords, max(self2_sel[,1]))$nn.index
for(i in 1:length(coords[,1])){
  if(self2_sel[i,1]>1){
    knns_sub<-knns[i,1:(self2_sel[i,1])]
    for(ii in 1:length(knns_sub)){
      a <- c(coords[i,1], coords[i,2], coords[knns_sub[ii],1], coords[knns_sub[ii],2])
      A <- as.data.frame(rbind(A,a))
    }
  }
}

colnames(A) <- c("Sup_lon","Sup_lat","Dem_lon","Dem_lat")

# connectivity for energy sharing
g1 <- ggplot()+
  geom_sf(data=buildings2_map)+
  #annotation_map_tile(type="cartodark", zoomin=14)+
  geom_sf(data=coord2_sf, size=0.05)+
  geom_segment(data=A, 
               aes(x=Sup_lon, y=Sup_lat, xend=Dem_lon, yend=Dem_lat),
               color="blue",size=0.8, alpha=0.2)+
  labs(x="Longitude", y="Latitude")+
  theme_bw()
g1

# Number of Sharing Neighbours
buildings3 <- cbind(buildings2, opt=self2[,1], No._Sharing_Neighbours=self2_sel[,1])
g2 <- ggplot()+
  geom_sf(data=buildings3, aes(fill=No._Sharing_Neighbours), lwd=0.1)+ 
  scale_fill_gradient(low="black", high="white")+
  labs(x="Longitude", y="Latitude")+
  theme_bw()
g2

5B. Visualization - Community Clustering

# graph network
B <- NULL
for(i in 1:length(coords[,1])){
  if(self2_sel[i,1]>1){
    knns_sub<-knns[i,1:(self2_sel[i,1])]
    for(ii in 1:length(knns_sub)){
      b <- c(i, knns_sub[ii])
      B <- as.data.frame(rbind(B,b))
    }
  }
}
B <- as.matrix(B)
SC <- make_graph(B)
# community clustering: based on betweenness 
SC_e <- graph_from_edgelist(B, directed=F)
SC_eb <- cluster_edge_betweenness(SC_e, directed=F) # Better but being burden: SC_opt <- cluster_optimal(SC) # max(SC_opt$modularity)

#number of clusters
nc <- length(SC_eb)

# 
C <- NULL
for(i in 1:nc){
  for(ii in 1:length(SC_eb[[i]])){
    c <- c(buildID=SC_eb[[i]][ii],group=i)
    C <- as.data.frame(rbind(C,c))
  }
}
C <- C[order(C$buildID),]

clusterbuidings <- bind_cols(buildings2,C)

# Clustered buildings
g3 <- ggplot()+
  geom_sf(data=clusterbuidings, aes(fill=as.factor(group)), lwd=0.1)+
  labs(x="Longitude", y="Latitude")+
  theme_bw()
g3

5C. Visualization - Boundaries using Voronoi Diagram

coord_g_CENTER <- NULL
for(i in 1:nc){
  temp <- filter(clusterbuidings,group==i)
  coord_g <- st_set_geometry(temp[,c("long","lat")],NULL) # meter scale
  coord_g_center <- c(i, apply(coord_g,2,mean))
  coord_g_CENTER <- as.data.frame(rbind(coord_g_CENTER,coord_g_center))
}

temp <- coord_g_CENTER
coord_g_CENTER <- st_as_sf(coord_g_CENTER, coords=c("long","lat"), crs = 6677)
boundbuildings <- st_transform(clusterbuidings, crs=6677)
coord_g_CENTER <- st_as_sf(left_join(temp,coord_g_CENTER,by="V1"))
boundbuildings <- st_transform(boundbuildings, crs=6668)
coord_g_CENTER <- st_transform(coord_g_CENTER, crs=6668)
coord_g_CENTER$PX <- st_coordinates(coord_g_CENTER)[,1]
coord_g_CENTER$PY <- st_coordinates(coord_g_CENTER)[,2]

outline <- as.data.frame(matrix(c(
  min(coords[1]), min(coords[2]),
  min(coords[1]), max(coords[2]),
  max(coords[1]), max(coords[2]),
  max(coords[1]), min(coords[2])),nrow=4,ncol=2,byrow=T))

g4 <- ggplot() +
  geom_sf(data=boundbuildings, aes(fill=as.factor(group)), lwd=0.1)+
  stat_voronoi(data=coord_g_CENTER,mapping=aes(long,lat), geom="path",outline=outline, lwd=1, lty="dashed") +
  geom_point(data=coord_g_CENTER,mapping=aes(long,lat)) +
  labs(x="Longitude", y="Latitude")+
  theme_bw()
g4

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.