GIS Coursework Tutorial

Estimation of London’s Rooftop Solar Energy Potential

In this tutorial, technical solar rooftop PV panels will be estimated. Study area is London, UK. The analysis has two main stages:

QGIS program will be used to characterize which parts of rooftop area is suitable for PV installation.
R studio will be used to calculate energy output of suitable PV area installed PV panel

To conduct analysis:

LIDAR data ¹ https://environment.data.gov.uk/DefraDataDownload/?Mode=survey
Landcover data https://www.ordnancesurvey.co.uk/business-government/products/open-map-local
London Wards data https://data.london.gov.uk/dataset/statistical-gis-boundary-files-london
Energy Consumption data https://data.london.gov.uk/dataset/electricity-consumption-borough

QGIS Part

Aspect Analysis.

At first, let’s add raster LIDAR layer in order to do aspect and shadow analysis.

Layer -> Add Layer -> Add Raster Layer.

—

Process lidar asc into aspect raster

By using QGIS Toolbox: GDAL -> Aspect

—

Reclassify aspect layer to simplify values by Raster Calculator.

(“Aspect@1” < 135) * 0 + (“Aspect@1” <= 225) * 1 + (“Aspect@1” > 225) * 0

This process will also help us to reduce large file size.

—

As a result, we have binary raster which norther facing parts have 0 value and others have 1 value.

Shadow Analysis

To conduct shadow analysis, UMEP plugin is required.

-> UMEP -> Processor -> Solar Radiation -> Daily Shadow Patter -> Adding LIDAR data

> Comparing shadow layer and aspect layer

With raster calculator suitable area for rooftop will have 1 value.

“asp@1” * “Shadow_at_20210111_1200@1” * “0.25”

Because each cell is 50 cm X 50cm, we multiply it by 0.25 cm^2. For the further process, it will help us to calculate area.

Adding londonbuildings.shp vector layer

-> Using SAGA -> Raster Statistics for polygons -> Selecting londonbuildings polygon and rasters -> Number of cells and sum to calculate suitable area on each building polygon. Then we can export results and import to R studio for calculations.

—

R Studio

# required libraries
library(raster)

## Loading required package: sp

library(rgdal)

## Warning: package 'rgdal' was built under R version 4.0.3

## rgdal: version: 1.5-18, (SVN revision 1082)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/erisf/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/erisf/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

library(rgeos)

## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-4 
##  Polygon checking: TRUE

library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, union

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(here)

## here() starts at C:/Users/erisf/OneDrive - University College London/CASA/Modules/GIS/GIS Coursework/R studio/Coursework/Gis-Coursework

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

## v tibble  3.0.4     v purrr   0.3.4
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.0

## Warning: package 'tibble' was built under R version 4.0.3

## Warning: package 'readr' was built under R version 4.0.3

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()

library(spatstat)

## Warning: package 'spatstat' was built under R version 4.0.3

## Loading required package: spatstat.data

## Warning: package 'spatstat.data' was built under R version 4.0.3

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## The following object is masked from 'package:raster':
## 
##     getData

## Loading required package: rpart

## Registered S3 method overwritten by 'spatstat':
##   method     from
##   print.boxx cli

## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'

## 
## Note: spatstat version 1.64-1 is out of date by more than 8 months; we recommend upgrading to the latest version.

## 
## Attaching package: 'spatstat'

## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift

library(sp)
library(rgeos)
library(maptools)

## Checking rgeos availability: TRUE

library(GISTools)

## Warning: package 'GISTools' was built under R version 4.0.3

## Loading required package: RColorBrewer

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:spatstat':
## 
##     area

## The following object is masked from 'package:dplyr':
## 
##     select

## The following objects are masked from 'package:raster':
## 
##     area, select

library(tmap)
library(geojson)

## 
## Attaching package: 'geojson'

## The following object is masked from 'package:graphics':
## 
##     polygon

library(geojsonio)

## 
## Attaching package: 'geojsonio'

## The following object is masked from 'package:base':
## 
##     pretty

library(tmaptools)

Importing necessary data…

#london boroughs and wards datasets

londonwards <- st_read(here("data","raw", "London_Ward.shp"))

## Reading layer `London_Ward' from data source `C:\Users\erisf\OneDrive - University College London\CASA\Modules\GIS\GIS Coursework\R studio\Coursework\Gis-Coursework\data\raw\London_Ward.shp' using driver `ESRI Shapefile'
## Simple feature collection with 649 features and 7 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## projected CRS:  OSGB 1936 / British National Grid

londonboroughs <- st_read(here("data","raw", "London_Borough_Excluding_MHW.shp"))

## Reading layer `London_Borough_Excluding_MHW' from data source `C:\Users\erisf\OneDrive - University College London\CASA\Modules\GIS\GIS Coursework\R studio\Coursework\Gis-Coursework\data\raw\London_Borough_Excluding_MHW.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## projected CRS:  OSGB 1936 / British National Grid

#reading london rooftop data processed in QGIS

rooftoplondon <- st_read(here("data","qgisprocessed", "londonpvstatistics.gpkg"))

## Reading layer `londonpvstatistics' from data source `C:\Users\erisf\OneDrive - University College London\CASA\Modules\GIS\GIS Coursework\R studio\Coursework\Gis-Coursework\data\qgisprocessed\londonpvstatistics.gpkg' using driver `GPKG'
## Simple feature collection with 669310 features and 31 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZM
## bbox:           xmin: 503970.2 ymin: 156581.2 xmax: 561188.1 ymax: 200024.3
## z_range:        zmin: 0 zmax: 0
## m_range:        mmin: 0 mmax: 0
## projected CRS:  OSGB 1936 / British National Grid

#Replacing null values with 0 
rooftoplondon <- rooftoplondon %>%
  replace(is.na(.), 0)

Because the LIDAR file of London as total is very large size, approxiamately 8.5 gb. We divided London into 10 parts to calculate fast. Therefore, in the following step, sutaible PV area is calculated

#adding total potential pv area column to dataframe 

rooftoplondon$pvarea <- rooftoplondon$X1..SUM.+ rooftoplondon$X2..SUM.+rooftoplondon$X3..SUM+rooftoplondon$X4..SUM.+rooftoplondon$X5..SUM.+
  rooftoplondon$X6..SUM.+rooftoplondon$X7..SUM.+rooftoplondon$X8..SUM.+rooftoplondon$X9..SUM.+rooftoplondon$X11..SUM.

Next step is calculating each building footprint area.

#calculating total rooftop area of each building 

rooftoplondon$area <- st_area(rooftoplondon)
#Transforming to integer 
rooftoplondon$area <- gsub("C([0-9]+)_.*", "\\1", rooftoplondon$area)%>%
  as.integer(rooftoplondon$area)

There are missing values in LIDAR data resulted in building rooftop area as 0 Therefore, I estimate them in a traditional way. Multiply rooftop which are 0 by general ratio and removing decimals because we have 1X1

#calculating total rooftop area of each building 
ratio <- sum(rooftoplondon$pvarea)/sum(rooftoplondon$area[rooftoplondon$pvarea!=0])
#General ratio is 0.70
rooftoplondon <- rooftoplondon %>%
  mutate(pvareanew = ifelse(pvarea == 0,lapply(area*ratio, as.integer),pvarea))

#Calculating total pv area in London 
rooftoplondon$pvareanew <- unlist(rooftoplondon$pvareanew)
totalpvarea <- sum(rooftoplondon$pvareanew) # we have 175.826.108 m^2 available roofplane area

getting center points of each building in order to find which wards they are from

centroids <- rooftoplondon %>%
  st_zm()%>%
  st_centroid()

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

#ward intersection for each building 
centroids <- st_set_crs(centroids, 27700)
londonwards <- st_transform(londonwards, 27700)
rooftoplondonwards <- st_intersection(centroids, londonwards)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

#Removing unnecessary columns 

rooftoplondonwards <- subset(rooftoplondonwards, select = c("ID","GSS_CODE.1","NAME.1", "BOROUGH", "area", "pvareanew"))

#calculating total potential rooftop area grouping by each wards

rooftoplondonwards <- rooftoplondonwards %>% 
  group_by(GSS_CODE.1) %>%
  summarise(sumpvarea = sum(pvareanew))

## `summarise()` ungrouping output (override with `.groups` argument)

#creating final dataframe 
finallondonwards <- londonwards %>% arrange(GSS_CODE)
finallondonwards$pvarea <- rooftoplondonwards$sumpvarea

Calculating Solar Electricty Output and CO^2 reduce output

finallondonwards$kWp <- (finallondonwards$pvarea/1.7)*200/1000
finallondonwards$Ep <- finallondonwards$kWp*0.75*1108*0.70
finallondonwards$EpAsGWh <- finallondonwards$Ep/10^6

sum(finallondonwards$EpAsGWh) # sum of amount of electricty provided by solar panels = 17603.38 GWh

## [1] 12032.71

#calculating greenhouse emission reduced
finallondonwards$EpAsMWh <- finallondonwards$EpAsGWh*10^6
finallondonwards$MtCO2 <- finallondonwards$EpAsMWh*0.2977/10^9

sum(finallondonwards$MtCO2) #773341.4 tCo2 is reduced

## [1] 3.582138

consumptiondata <- read.csv(here("data","raw", "energyconsumption.csv"))

#Turning to borough scale
finallondonboroughcalculated <- finallondonwards %>% 
  group_by(LB_GSS_CD) %>%
  summarise(sumpvarea = sum(pvarea), sumMtCO2 = sum(MtCO2), sumGWH = sum(EpAsGWh))

## `summarise()` ungrouping output (override with `.groups` argument)

finallondonboroughcalculated$consumption <- consumptiondata$consumption

finallondonboroughcalculated$rate <- finallondonboroughcalculated$sumGWH/finallondonboroughcalculated$consumption

sum(finallondonboroughcalculated$consumption)

## [1] 37800.2

finallondonboroughcalculated$GWHasinteger <- as.integer(finallondonboroughcalculated$sumGWH)
finallondonboroughcalculated$percent <- as.integer(finallondonboroughcalculated$rate*100)

Plotting

tm1 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("sumGWH",
              palette="seq",
              style="jenks",
              midpoint=NA,
              border.col = "white", 
              border.alpha = 0.2,
              popup.vars=c("LB_GSS_CD", "sumGWH"),
              title="Percentage") +
  tm_legend(show=FALSE)

legend1 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("sumGWH",
              palette="seq",
              style="jenks",
              midpoint=NA,
              popup.vars=c("LB_GSS_CD", "sumGWH"),
              title="Percentage") +
  tm_layout(legend.only = TRUE, legend.position=c(0.2,0.25),asp=0.1)+
  tm_compass(north=0, position=c(0.65,0.6))

t1=tmap_arrange(tm1, legend1)
t1

tmap_save(t1, 'GWHjenks.png')

## Map saved to C:\Users\erisf\OneDrive - University College London\CASA\Modules\GIS\GIS Coursework\R studio\Coursework\Gis-Coursework\GWHjenks.png

## Resolution: 2100 by 2100 pixels

## Size: 7 by 7 inches (300 dpi)

tm2 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("percent",
              palette="PuBu",
              midpoint=NA,
              border.col = "white", 
              border.alpha = 0.2,
              popup.vars=c("LB_GSS_CD", "percent"),
              title="Percentage") +
  tm_legend(show=FALSE)

legend2 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("percent",
              palette="PuBu",
              midpoint=NA,
              popup.vars=c("LB_GSS_CD", "percent"),
              title="Percentage") +
  tm_layout(legend.only = TRUE, legend.position=c(0.2,0.25),asp=0.1)+
  tm_compass(north=0, position=c(0.65,0.6))
  
t2=tmap_arrange(tm2, legend2)
t2

tmap_save(t2, 'percentage.png')

## Map saved to C:\Users\erisf\OneDrive - University College London\CASA\Modules\GIS\GIS Coursework\R studio\Coursework\Gis-Coursework\percentage.png

## Resolution: 2100 by 2100 pixels

## Size: 7 by 7 inches (300 dpi)

tm3 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("sumMtCO2",
              palette="BuPu",
              midpoint=NA,
              border.col = "white", 
              border.alpha = 0.2,
              popup.vars=c("LB_GSS_CD", "sumMtCO2"),
              title="Percentage") +
  tm_legend(show=FALSE)

legend3 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("sumMtCO2",
              palette="BuPu",
              midpoint=NA,
              popup.vars=c("LB_GSS_CD", "sumMtCO2"),
              title="Percentage") +
  tm_layout(legend.only = TRUE, legend.position=c(0.2,0.25),asp=0.1)+
  tm_compass(north=0, position=c(0.65,0.6))

t3=tmap_arrange(tm3, legend3)
t3

tmap_save(t3, 'CO2NEW.png')

## Map saved to C:\Users\erisf\OneDrive - University College London\CASA\Modules\GIS\GIS Coursework\R studio\Coursework\Gis-Coursework\CO2NEW.png

## Resolution: 2100 by 2100 pixels

## Size: 7 by 7 inches (300 dpi)

tm4 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("sumMtCO2",
              palette="PuBu",
              midpoint=NA,
              border.col = "white", 
              border.alpha = 0.2,
              popup.vars=c("LB_GSS_CD", "sumMtCO2"),
              title="Percentage") +
  tm_legend(show=FALSE)

legend4 <- tm_shape(finallondonboroughcalculated) +
  tm_polygons("sumMtCO2",
              palette="PuBu",
              midpoint=NA,
              popup.vars=c("LB_GSS_CD", "sumMtCO2"),
              title="Percentage") +
  tm_layout(legend.only = TRUE, legend.position=c(0.2,0.25),asp=0.1)+
  tm_compass(north=0, position=c(0.65,0.6))

t4=tmap_arrange(tm4, legend4)
t4

tmap_save(t4, 'CO21.png')

## Map saved to C:\Users\erisf\OneDrive - University College London\CASA\Modules\GIS\GIS Coursework\R studio\Coursework\Gis-Coursework\CO21.png

## Resolution: 2100 by 2100 pixels

## Size: 7 by 7 inches (300 dpi)

tmap_mode("plot")

## tmap mode set to plotting

—

LIDAR Data divided into cells for London that contains almost 120 file. In the following part of analysis, processed data will be shared.↩︎

GIS Coursework Tutorial

Ömer Faruk Eriş

11/01/2021

Estimation of London’s Rooftop Solar Energy Potential

QGIS Part

R Studio