Delhi Project Visualizations

Daily Peak Loads from 2010-2013

# Read in load data
setwd("~/Documents/SLDC_big_data/Extracted Data")

load2010 <- read.csv('loads2010-2011.csv')
load2010 <- load2010[,-c(1,6)] # remove 'X' and 'timepoint' variable

load2011 <- read.csv('loads2011-2012.csv')
load2011 <- load2011[,-c(1,6)] # remove 'X' and 'timepoint' variable

load2012 <- read.csv('loads2012-2013.csv')
load2012 <- load2012[,-c(1,6)] # remove 'X' and 'timepoint' variable

# Combine all observations from 2010 - 2013
loadAll <- rbind(load2010, load2011, load2012)
loadAll$date <- as.Date(loadAll$date)
table(loadAll$discom)
summary(loadAll$date)

# Get dataset with daily peak load over time
require(plyr)
dailyPeak <- ddply(.data=loadAll, .variables=.(date, discom), subset, mW==max(mW))

# We have an issue where if an idential reading was recorded at different timepoints,
# then we get 2 records for daily peak value
dateOccur <- as.data.frame(table(dailyPeak$date))
dateOccur[which(dateOccur$Freq==6),]

# Plot daily peak load over time by discom
require(ggplot2)

# Plot all years
# All years: Daily peak
ggplot(dailyPeak, aes(x=date, y=mW, color=discom)) + 
  geom_point() + 
  facet_grid(discom ~ ., scales="free_y") +
  scale_x_date(breaks="3 months") +
  ggtitle("Daily Peak Loads from 2010-03-29 to 2013-03-10")

# Daily peak: Plot only Year 1
ggplot(subset(dailyPeak, date<"2011-03-30"), aes(x=date, y=mW, color=discom)) + 
  geom_line(size=1.5) + 
  facet_grid(discom ~ ., scales="free_y") +
  scale_x_date(breaks="2 months") +
  ggtitle("Daily Peak Loads from 2010-03-29 to 2011-03-29")

# All recordings: Plot only Year 1
ggplot(subset(loadAll, date<"2011-03-30"), aes(x=date, y=mW, color=discom)) + 
  geom_point() + 
  facet_grid(discom ~ ., scales="free_y") +
  scale_x_date(breaks="3 months")

### Future plots
# At what times of the day do peak loads occur?
# How does a graph from the daily peak load differ from the daily minimum load?

Delhi Districts

# Mapping Delhi districts

# Install packages
install.packages('rgdal', type='source')
install.packages('sp')
install.packages('maps')
install.packages('maptools')
install.packages('rgeos', type='source')
install.packages('spatstat')
install.packages('deldir')
install.packages('devtools')
install.packages('base64enc')

require('devtools')
install_github('ramnathv/rCharts@dev')
install_github('ramnathv/rMaps')

# Load packages
require(rgdal)
require(sp)
require(maps)
require(maptools)
require(rgeos)
require(spatstat)
require(deldir)
require(rCharts)
require(rMaps)

###########
# Set directory to folder containing .shp and other files
setwd("~/Documents/SLDC_big_data/MappingFiles")

# See: 'Read OGR vector maps into Spatial objects'
ogrListLayers("Districts.shp") # this lists the layers
districts=readOGR("Districts.shp", layer="Districts") # this produces a SpatialPolygonsDataFrame

# Understanding the .shp file
class(districts)
getClass("SpatialPolygonsDataFrame")
row.names(districts)
str(districts@polygons[3],max.level=3)

# Plot shape
plot(districts)

# Create fill colors for districts
require('RColorBrewer')
distriClr <- brewer.pal(9, "Pastel1")
plot(districts, col=distriClr, lwd=3.0)
legend("right", legend=districts@data$DISTRICT, fill=distriClr, bty="n",
       y.intersp=0.6, x.intersp=0.8 )

?legend

# Read in data containing electricity grid coordinates
setwd("~/Documents/SLDC_big_data")
gridData <- read.csv('grid_coord.csv')
gridCoord <- data.frame(x=gridData$long, y=gridData$lat)

# Slightly modified Voronoi function, takes an additional 
# spatial polygons argument and extends to that box:
# Source: https://stackoverflow.com/questions/12156475/combine-voronoi-polygons-and-maps/12159863#12159863
voronoipolygons <- function(x,poly) {
  require(deldir)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  bb = bbox(poly)
  rw = as.numeric(t(bbox(poly)))
  z <- deldir(crds[,1], crds[,2],rw=rw)
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  require(sp)
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  
  voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
                                                          y=crds[,2], 
                                                          row.names=sapply(slot(SP, 'polygons'), 
                                                          function(x) slot(x, 'ID'))))
  return(voronoi)  
}

# Create Voronoi polygons for grid areas
gridArea <-voronoipolygons(gridCoord, districts)

# Plot Voronoi polygons for Delhi
plot(gridArea, lty='dashed', add=TRUE)

# Plot points of grid locations and label them by discom
points(x=gridData$long, y=gridData$lat, cex=1.8)
text(x=gridData$long, y=gridData$lat, labels=gridData$griddis, cex=0.5, pos=4, offset = 0.5, col=gridData$Discom)


#########################################################################################################
### This is the part that creates the intersection of the Voronoi polygons and Delhi district boundaries
# Set the proj4 strings equal to one another
proj4string(gridArea) <- proj4string(districts)

# Find intersection of (1) polygons based on grid coordinates and (2) Delhi district boundaries
gridDelhi = gIntersection(districts, gridArea, byid=TRUE)

# Plot the intersection
plot(gridDelhi)
###
###########################################################################################################

# Plot Delhi outer boundary
delhi <- gUnionCascaded(districts)
plot(delhi)

# Set the proj4 strings equal to one another
proj4string(gridArea) <- proj4string(delhi)

# Find intersection of (1) polygons based on grid coordinates and (2) Delhi district boundaries
gridDelhi = gIntersection(delhi, gridArea, byid=TRUE)

# Plot the intersection
plot(gridDelhi)
polygonClr <- palette(brewer.pal(9, "Pastel1"))
plot(gridDelhi, col=polygonClr)
points(x=gridData$long, y=gridData$lat, col='black', pch=19, cex=1.8)
text(x=gridData$long, y=gridData$lat, labels=gridData$grid_name, cex=0.6, pos=2, col=gridData$Discom)

Delhi Weekly Consumption Map

# This file reads in weekly net consumption data
# from the 22 grids listed in the file (units: Mus)
# 'company wise details.xls'
# NDMC and MES were added separately

# Read in file with weekly consumption data
setwd("~/Documents/SLDC_big_data")
data2012 <- read.csv('weeklyNetConsumption2012.csv', header=TRUE)

# Read in file with grid coordinates
gridData <- read.csv('grid_coord3.csv', header=TRUE)
gridnames <- as.character(gridData$name)
gridData$state_code <- as.numeric(gridData$state_code)

# Rename columns to match names in gridData
colnames(data2012) <- c(colnames(data2012[,1:5]), gridnames)

# Melt dataframe so that 'grid' is a variable
# install.packages("reshape")
library("reshape")
longData <- melt(data2012, id=c(colnames(data2012[,1:5])))

# rename column names
colnames(longData) <- c(colnames(longData[,1:5]),"name", "consump")


############## Creating maps ###############
# First create shapefile to be turned into TopoJSON
# Set directory to folder containing Delhi .shp and other files
setwd("~/Documents/SLDC_big_data/MappingFiles")
districts=readOGR("Districts.shp", layer="Districts") # this produces a SpatialPolygonsDataFrame
plot(districts)

# Create data frame for grid coordinates
gridCoord <- data.frame(x=gridData$long, y=gridData$lat)

# Slightly modified Voronoi function, takes an additional spatial polygons argument and extends to 
# that box. Source: https://stackoverflow.com/questions/12156475/combine-voronoi-polygons-and-maps/12159863#12159863
voronoipolygons <- function(x,poly) {
  require(deldir)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  bb = bbox(poly)
  rw = as.numeric(t(bbox(poly)))
  z <- deldir(crds[,1], crds[,2],rw=rw)
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  require(sp)
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  
  voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
                                                          y=crds[,2], 
                                                          row.names=sapply(slot(SP, 'polygons'), 
                                                                           function(x) slot(x, 'ID'))))
  return(voronoi)  
}

# Create Voronoi polygons for grid areas
gridArea <-voronoipolygons(gridCoord, districts)

# Create Delhi outer boundary SpatialPolygons object
delhi <- gUnionCascaded(districts)
plot(delhi)

# Set the proj4 strings for gridArea and delhi equal to one another 
# so that we can create the intersection
proj4string(gridArea) <- proj4string(delhi)

# Find intersection of (1) polygons based on grid coordinates and (2) Delhi district boundaries
gridDelhi = gIntersection(delhi, gridArea, byid=TRUE)

# Plot the intersection
plot(gridDelhi)

# Turn SpatialPolygons object into SpatialPolygonsDataFrame, so that it can
# be converted to a shapefile
row.names(gridDelhi) <- row.names(gridData) 
gridDelhiDF <- SpatialPolygonsDataFrame(gridDelhi, data=gridData)
plot(gridDelhiDF)

# Create IDs that will be used by 'ichoropleth()' function to determine fill color
gridDelhiDF <- spChFIDs(gridDelhiDF, paste0("p",gridData$state_code))
gridDelhiDF$id <- gridData$id <- paste0('p',gridData$state_code)
str(gridDelhiDF @ data)

# Add ID to dataset containing peak consumption data
longData$id <- rep(gridData$id, each=53)

# Create shapefile and convert that to TopoJSON using terminal commands
setwd("~/Documents/SLDC_big_data/")
writeOGR(obj=gridDelhiDF, dsn = 'DelhiWeekly', layer = 'myPolygons', driver = "ESRI Shapefile")

# Here is the terminal command to convert the .shp file into a TopoJSON file
# cd ~/Documents/SLDC_big_data/DelhiWeekly
# topojson -o delhigrid.json  myPolygons.shp -p id=id,name=name  --id-property=id


#########This part creates the choropleth and saves and produces an HTML file##########
# Add week number to data
longData$week <- rep(1:53, times=22)

setwd("~/Documents/SLDC_big_data/DelhiWeekly")

delhiMap <- ichoropleth(consump ~ id, data=longData, ncuts=9, 
                        pal='YlOrRd', animate = 'week', 
                        map = 'myPolygons')

delhiMap$set(
  geographyConfig = list(
    dataUrl = "delhigrid.json"
  ),
  scope = 'myPolygons',
  setProjection = '#! function( element, options ) {
   var projection, path;
   projection = d3.geo.mercator()
    .center([77.3, 28.3]).scale(30000)
    .translate([element.offsetWidth / 2, element.offsetHeight]);

   path = d3.geo.path().projection( projection );
   return {path: path, projection: projection};
  } !#'
)

delhiMap$save('delhi.html', cdn = TRUE)

Temperature Load Curves

# Read in temperature data
setwd("~/Documents/SLDC_big_data/Extracted Data")
weather <- read.csv('2009_2013_Delhi_Temp_Wunderground.csv')
weather$IST <- as.Date(weather$IST, "%m/%d/%Y")
summary(weather$IST)

# Keep only temp
temp <- weather[,1:4]
colnames(temp) <- c("date","max", "mean", "min")

# Read in load data
setwd("~/Documents/SLDC_big_data/Extracted Data")
load2011 <- read.csv('loads2011-2012.csv')
load2011$date <- as.Date(load2011$date)
summary(load2011$date)

load2012 <- read.csv('loads2012-2013.csv')

# Select daily maximum load reading
max2011 <- aggregate(mW ~ discom + date, data = load2011, FUN = "max")
max2012 <- aggregate(mW ~ discom + date, data = load2012, FUN = "max")

# Merge 2 files together
loads <- rbind(max2011, max2012)
loads$date <- as.Date(loads$date)
summary(loads$date)

# Graph data
ggplot(loads, aes(x=date, y=mW, color=discom)) + 
  geom_line(size=1.1) + 
  facet_grid(discom ~ ., scales="free_y") +
  scale_x_date(breaks="1 month")

# Subset and just get 2012 loads
only2012 <- subset(loads, date > "2011-12-31" & date < "2013-01-01")

# Keep only 2012 temperatures
tempOnly <-temp[temp$date > "2011-12-31" & temp$date < "2013-01-01",]

# Create dataframe that repeats temperature values 5 times for 5 discoms
temp12 <- data.frame(1:1830)
temp12$date <- rep(tempOnly$date, each=5)
temp12$max <- rep(tempOnly$max, each=5)
temp12$min <- rep(tempOnly$min, each=5)
temp12$mean <- rep(tempOnly$mean, each=5)

# Add temperature variable to 'loads' data
tempload12 <- cbind(only2012, temp12)

# Create temperature-load curves based on maximum daily loads
# and maximum daily temperatures from 
ggplot(tempload12, aes(x=max, y=mW, color=discom)) + 
  geom_point() +
  facet_wrap(~discom, ncol=5, scales="free_y") +
  stat_smooth(method="loess",se=FALSE, size=2)