1 Objectives :

• To create interactive map of Inverse Distance Weighted (IDW) interpolation predication for solar_layouts_COSMO and wind_layouts_COSMO
• The transmission network model is based on Jensen and Pinson (2017)
• The dataset is downloaded from https://zenodo.org/record/35177#.XVZqicszabI

2 Import all static data csv file

generator_info      =read.csv(paste(infilepath,"generator_info.csv",sep=""),header=T,sep=",")
network_edges       =read.csv(paste(infilepath,"network_edges.csv",sep=""),header=T,sep=",")
network_nodes       =read.csv(paste(infilepath,"network_nodes.csv",sep=""),header=T,sep=",")
network_hvdc_links  =read.csv(paste(infilepath,"network_hvdc_links.csv",sep=""),header=T,sep=",")
solar_layouts_COSMO =read.csv(paste(infilepath,"solar_layouts_COSMO.csv",sep=""),header=T,sep=",")
solar_layouts_ECMWF =read.csv(paste(infilepath,"solar_layouts_ECMWF.csv",sep=""),header=T,sep=",")
wind_layouts_COSMO  =read.csv(paste(infilepath,"wind_layouts_COSMO.csv",sep=""),header=T,sep=",")
wind_layouts_ECMWF  =read.csv(paste(infilepath,"wind_layouts_ECMWF.csv", sep=""),header=T,sep=",")

3 Transmission Network in EU zone: igraph

network_edges_N = network_edges
network_nodes_N = network_nodes
df_from_to = data.frame(network_edges_N[, 1:2] )
network_nodes_meta_df = network_nodes_N
g = graph.data.frame(df_from_to, directed = FALSE, vertices = network_nodes_meta_df$ID)
project_layout =  layout.norm(as.matrix(network_nodes_meta_df[,6:5]))
nodeigraph = plot.igraph(g, layout = project_layout,  vertex.size=1.2,vertex.label=NA, vertex.color="red" )

4 Create spatial object

coordinates(generator_info) = ~longitude + latitude
class(generator_info)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

proj4string(generator_info) = CRS("+init=epsg:4124")
extent(generator_info)

## class      : Extent 
## xmin       : -9.0545 
## xmax       : 28.607 
## ymin       : 36.1828 
## ymax       : 57.0736

generator_info

## class       : SpatialPointsDataFrame 
## features    : 969 
## extent      : -9.0545, 28.607, 36.1828, 57.0736  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4124 +proj=longlat +ellps=bessel +towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs 
## variables   : 13
## # A tibble: 969 x 13
##       ID name  country origin status primaryfuel secondaryfuel capacity
##    <int> <fct> <fct>    <int> <fct>  <fct>       <fct>            <dbl>
##  1  2069 Alma~ Spain      162 No Da~ Nuclear     Unknown         1957  
##  2  2142 Aram~ France     486 Opera~ Fuel Oil    Unknown         1400  
##  3  2154 Arge~ France     457 No Da~ Hydro       Unknown           47.5
##  4  2175 Asco~ Spain      124 No Da~ Nuclear     Unknown         2060  
##  5  2300 Bell~ France     357 Opera~ Nuclear     Unknown         2726  
##  6  2329 Bezn~ Switze~    990 No Da~ Nuclear     Unknown          760  
##  7  2336 Bibl~ Germany    932 No Da~ Nuclear     Unknown         2525  
##  8  2366 Blay~ France     451 No Da~ Nuclear     Unknown         3640  
##  9  2368 Blé~ France     278 Opera~ Coal        Fuel Oil         750  
## 10  2385 Bohu~ Slovak~   1215 No Da~ Nuclear     Unknown          880  
## # ... with 959 more rows, and 5 more variables: lincost <dbl>,
## #   cyclecost <dbl>, minuptime <int>, mindowntime <int>,
## #   minonlinecapacity <dbl>

5 Generator Database

nodes_from = data.frame(network_edges_N)[, c(1,2)]
df2= network_nodes_meta_df[, c(1,5,6)] 
colnames(df2)[1]  = "fromNode"
df12 = merge(x = nodes_from, y = df2, by = "fromNode", all = FALSE)
colnames(df12)[3:4] = c("Begin_lat",  "Begin_lon")


df2= network_nodes_meta_df[, c(1,5,6)] 
colnames(df2)[1]  = "toNode"
df123 = merge(x = df12 , y = df2, by = "toNode", all = FALSE)
colnames(df123)[5:6] = c("End_lat",  "End_lon")


begin.coord = data.frame(lon=df123[,4], lat=df123[,3])
end.coord   = data.frame(lon=df123[,6], lat=df123[,5])

line_sf = vector("list", nrow(begin.coord))
for (i in seq_along(line_sf)){line_sf[[i]] = st_linestring(as.matrix(rbind(begin.coord[i, ], end.coord[i,])))}

line_sfc = st_sfc(line_sf, crs = '+proj=longlat +ellps=bessel +towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs')

# lines_sp = as(line_sfc, "Spatial")
# lines_sf = st_sf( geometry = line_sfc)
plot(line_sfc, col = 'blue')


plot(line_sfc, col = 'blue')

colorlist = list("#ffff00","#1a0000", "#ff8000", "#00ffff", "#ff00ff",  "#ff0000"  , "#800000", "#ffe6e6","#a9a9a9" )
mapviewOptions(basemaps =  c( "OpenTopoMap", "CartoDB.Positron", "CartoDB.DarkMatter","OpenStreetMap", "Esri.WorldImagery" )  )

mapview(line_sfc, lwd = 1) +
mapview(generator_info, zcol = 'primaryfuel', legend = TRUE, col.regions= colorlist , burst = TRUE, cex = 4)

solar_layouts_COSMO_N = solar_layouts_COSMO
df2= network_nodes_meta_df[, c(1,5,6)] 
colnames(df2)[1]  = "node"
df_solar_COSMO = merge(x = solar_layouts_COSMO_N, y = df2, by = "node", all = FALSE)

df_solar_COSMO_sp = df_solar_COSMO
coordinates(df_solar_COSMO_sp) = ~longitude + latitude
class(df_solar_COSMO_sp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

proj4string(df_solar_COSMO_sp) = CRS("+init=epsg:4124")
extent(df_solar_COSMO_sp)

## class      : Extent 
## xmin       : -9.243608 
## xmax       : 28.96797 
## ymin       : 35.57556 
## ymax       : 57.2472

df_solar_COSMO_sp

## class       : SpatialPointsDataFrame 
## features    : 1494 
## extent      : -9.243608, 28.96797, 35.57556, 57.2472  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4124 +proj=longlat +ellps=bessel +towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs 
## variables   : 3
## # A tibble: 1,494 x 3
##     node Proportional Uniform
##    <int>        <dbl>   <dbl>
##  1     1        1262.   1183.
##  2     2        1463.   1387.
##  3     3         889.    776.
##  4     4         569.    489.
##  5     5        1145.   1007.
##  6     6        1154.   1062.
##  7     7        1170.   1016.
##  8     8         598.    525.
##  9     9        1630.   1388.
## 10    10         958.    868.
## # ... with 1,484 more rows

mapview(line_sfc, lwd = 1) + mapview(df_solar_COSMO_sp, zcol = 'Proportional', cex = 3)

6 IDW interpolation for solar_layouts_COSMO

max_Proportional = max(df_solar_COSMO$Proportional)
min_Proportional = min(df_solar_COSMO$Proportional)

max_latitude = max(df_solar_COSMO$latitude)
min_latitude = min(df_solar_COSMO$latitude)

max_longitude = max(df_solar_COSMO$longitude)
min_longitude = min(df_solar_COSMO$longitude)

latitude_range  = as.numeric(c(min_latitude, max_latitude))
longitude_range = as.numeric(c(min_longitude, max_longitude))

7 Inverse Distance Weighted (IDW) interpolation for solar_layouts_COSMO

x_range = longitude_range
y_range = latitude_range
grd = expand.grid(x = seq(from = x_range[1],to = x_range[2],  by = 1/10),
                  y = seq(from = y_range[1],to = y_range[2],  by = 1/10))
#class(grd)
coordinates(grd) = ~x + y
gridded(grd) =  TRUE

proj4string(grd) = CRS("+init=epsg:4124")
idw_solar_COSMO_Proportional = idw(formula = Proportional ~ 1, locations = df_solar_COSMO_sp, newdata = grd, idp = 1)

## [inverse distance weighted interpolation]

idw_solar_COSMO_Proportional_raster = raster(idw_solar_COSMO_Proportional)

mapview(idw_solar_COSMO_Proportional_raster, col = rainbow(5, rev = TRUE), alpha.regions = 0.3  ) +
mapview(line_sfc, lwd = 1) + mapview(df_solar_COSMO_sp, zcol = 'Proportional', legend = FALSE, cex = 3, col.regions = rainbow(5, rev = TRUE))

wind_layouts_COSMO_N = wind_layouts_COSMO
df2= network_nodes_meta_df[, c(1,5,6)] 
colnames(df2)[1]  = "node"
df_wind_COSMO = merge(x = wind_layouts_COSMO_N, y = df2, by = "node", all = FALSE)

df_wind_COSMO_sp = df_wind_COSMO
coordinates(df_wind_COSMO_sp) = ~longitude + latitude
class(df_wind_COSMO_sp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

proj4string(df_wind_COSMO_sp) = CRS("+init=epsg:4124")
extent(df_wind_COSMO_sp)

## class      : Extent 
## xmin       : -9.243608 
## xmax       : 28.96797 
## ymin       : 35.57556 
## ymax       : 57.2472

df_wind_COSMO_sp

## class       : SpatialPointsDataFrame 
## features    : 1494 
## extent      : -9.243608, 28.96797, 35.57556, 57.2472  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4124 +proj=longlat +ellps=bessel +towgs84=414.1,41.3,603.1,-0.855,2.141,-7.023,0 +no_defs 
## variables   : 3
## # A tibble: 1,494 x 3
##     node Proportional Uniform
##    <int>        <dbl>   <dbl>
##  1     1         812.   1260.
##  2     2        1190.   1917.
##  3     3         501.    827.
##  4     4         408.    521.
##  5     5         507.   1073.
##  6     6        1241.   1625.
##  7     7         662.   1083.
##  8     8         267.    559.
##  9     9        1020.   1479.
## 10    10         575.    979.
## # ... with 1,484 more rows

8 IDW interpolation for wind_layouts_COSMO

max_Proportional = max(df_wind_COSMO$Proportional)
min_Proportional = min(df_wind_COSMO$Proportional)

max_latitude = max(df_wind_COSMO$latitude)
min_latitude = min(df_wind_COSMO$latitude)

max_longitude = max(df_wind_COSMO$longitude)
min_longitude = min(df_wind_COSMO$longitude)

latitude_range  = as.numeric(c(min_latitude, max_latitude))
longitude_range = as.numeric(c(min_longitude, max_longitude))

9 Inverse Distance Weighted (IDW) interpolation wind_layouts_COSMO

x_range = longitude_range
y_range = latitude_range
grd = expand.grid(x = seq(from = x_range[1],to = x_range[2],  by = 1/10),
                  y = seq(from = y_range[1],to = y_range[2],  by = 1/10))
# class(grd)
coordinates(grd) = ~x + y
gridded(grd) =  TRUE

proj4string(grd) = CRS("+init=epsg:4124")
idw_wind_COSMO_Proportional = idw(formula = Proportional ~ 1, locations = df_wind_COSMO_sp, newdata = grd, idp = 1)

## [inverse distance weighted interpolation]

idw_wind_COSMO_Proportional_raster = raster(idw_wind_COSMO_Proportional)

mapview(idw_wind_COSMO_Proportional_raster, col = heat.colors(5, rev = TRUE), alpha.regions = 0.4  ) +
mapview(line_sfc, lwd = 1) + mapview(df_wind_COSMO_sp, zcol = 'Proportional', legend = FALSE, cex = 3, col.regions = heat.colors(5, rev = TRUE))