快速上手：停電預測挑戰賽.Rmd

---
title: "快速上手：停電預測挑戰賽"
author: "Li Cheng En"
date: "10/6/2017"
output: html_document
---

```{r setup, include=FALSE}
library(tidyverse)
library(tidyr)
#library(reshape)
#library(ggthemr)
#ggthemr("light")
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
#load("/Volumes/SD/DSP/Stray_Dog/Inform_data.rdata")
```

## 載入套件
```{r} 
library(tidyverse)
```

## 載入颱風停電戶資料
```{r} 
train <- read.csv("./data/train.csv")
submit <- read.csv("./data/submit.csv")
```

看一下資料長什麼樣子
```{r } 
head(train, 20)
```

```{r} 
head(submit, 20)
```

## 加入電桿資料
資料來源為[政府資料開放平台](https://data.gov.tw/dataset/33305)
在讀取電桿資料前，先做出電桿資料的路徑。
```{r }
poleString <- c("北北區處pole.csv", "嘉義區處pole.csv", "澎湖區處pole.csv", "北南區處pole.csv", "基隆區處pole.csv", "花蓮區處pole.csv", "北市區處pole.csv", "宜蘭區處pole.csv", "苗栗區處pole.csv", "北西區處pole.csv", "屏東區處pole.csv", "金門區處pole.csv", "南投區處pole.csv", "彰化區處pole.csv", "雲林區處pole.csv", "台中區處pole.csv", "新營區處pole.csv", "馬祖區處pole.csv", "台南區處pole.csv", "新竹區處pole.csv", "高雄區處pole.csv", "台東區處pole.csv", "桃園區處pole.csv", "鳳山區處pole.csv")

# 先宣告一條vector，名字是pole_wd
pole_wd <- c() 
# 以字串處理，將檔名接成完整的檔案路徑，存到pole_wd中
for(i in 1:length(poleString)){
  pole_wd[i] <- paste0("./data/poledata/", poleString[i])
}
```

讀取電桿檔案
```{r}
# 宣告一張list，預計將多個電桿資料檔一一讀取、存放到裡面
pole <- list()
for(i in 1:length(pole_wd)){
  pole[[i]] <- read.csv(pole_wd[i], # 檔案讀取路徑
                        header = TRUE, # 保留header
                        stringsAsFactors = FALSE, # 字串不要轉成類別型資料
                        sep = "\t" # 在csv檔裡是以"\t"分開
  )
  pole[[i]] <- pole[[i]][1:5] # csv檔中，只有第1~5欄是我們需要的，選取出來
}
pole <- Reduce(x = pole, f = rbind) # 將list合併為一張dataframe
```

清理電桿資料：先對「縣市」、「行政區」兩個欄位的類別清理乾淨。
```{r}
pole$縣市 <- as.factor(pole$縣市) # 轉成類別型資料
pole$行政區 <- as.factor(pole$行政區) 
levels(pole$縣市)[30] <- c("台南市") # 「臺」改成「台」
levels(pole$縣市)[28] <- c("台中市")
levels(pole$縣市)[29] <- c("台北市")
levels(pole$行政區)[332] <- c("頭份市") # 「頭份鎮」改成「頭份市」
```

新增一個「縣市行政區」的欄位
```{r}
pole$縣市 <- as.character(pole$縣市)
pole$縣市行政區 <- paste0(pole$縣市, pole$行政區)
pole$縣市行政區 <- pole$縣市行政區 %>% as.factor()
```

清理電桿型式欄位
```{r}
pole$型式 <- pole$型式 %>% as.factor()
# 將「縣市行政區」欄位為不合理值的樣本挑出來
dropSet <- filter(pole, 縣市行政區 == "分#6"|縣市行政區 == "低2"|
                    縣市行政區 == "低1"|縣市行政區 == "9)"|
                    縣市行政區 == "5"|縣市行政區 == "2左低1 "|
                    縣市行政區 == "2左低1"|縣市行政區 == "1右低1"|
                    縣市行政區 == "26"|縣市行政區 == "2)"|
                    縣市行政區 == "2"|縣市行政區 == "1右低1 "|
                    縣市行政區 == "11右2"|縣市行政區 == "11"|
                    縣市行政區 == "1"|縣市行政區 == ")"|
                    縣市行政區 == " ")
pole <- setdiff(pole, dropSet) # 將剛剛選出來的樣本濾掉
```

計算每個縣市鄉鎮區各有哪些種類、各幾支的電桿。
這邊我只有做到「鄉鎮區」，建議可以再進一步細做「村里」的部分
```{r}
pole_type <- group_by(pole, 縣市行政區, 型式) %>% #對「縣市行政區」、「型式」分類
             summarise(n = n()) %>% # 計算每個縣市鄉鎮區各有哪些種類、各幾支的電桿
             ungroup() 
```

觀察一下，第一列可以刪除
```{r}
head(pole_type, 10)
pole_type <- pole_type[-1,]
```

另外，「型式」為遺失值的狀況似乎會一直出現，可以將長型資料翻轉成寬型資料觀察一下趨勢。
```{r}
head(pole_type, 20)
pole_type <- spread(pole_type, # 將長型資料翻轉成寬型資料，方便清理
                    key = 縣市行政區, # 將「縣市行政區」的level作為欄位
                    value = n, # 資料
                    fill = 0) # NA值補0
pole_type <- t(pole_type[,-1]) %>% 
               as.data.frame()
# 將欄位名改成英文，較不會出錯
names(pole_type) <- c("型式", "pole1", "pole2", "pole3", "pole4", "pole5",
                      "pole6", "pole7", "pole8", "pole9", "pole10")
pole_type$縣市行政區 <- rownames(pole_type)
pole_type <- pole_type[, -1]
```

將電桿資料與train set資料集結合，就可以知道各鄉鎮區有哪些種類的電桿、各有幾支，這些地區在train set的各颱風中，又分別有幾戶停電。
```{r}
train$縣市行政區 <- paste0(train$CityName, train$TownName) 
train <- left_join(train, pole_type, by = "縣市行政區")
for(i in 14:23){
  train[,i][is.na(train[,i]) == TRUE] <- 0
}
```

## 加入人口戶數資料
資料來源為[政府資料開放平臺](https://data.gov.tw/dataset/32973#r0)

讀取資料
```{r}
family <- read.csv("./data/opendata10603M030.csv")
```

清理資料
```{r}
family <- family[-1, c(2,4)] # 只會用到第2~4欄的資料
family$site_id <- gsub(x = family$site_id, # 把「臺」改成「台」
                       pattern = "臺", 
                       replacement = "台")
family$site_id <- gsub(x = family$site_id, # 把「台東」改成「臺東」
                       pattern = "台東", 
                       replacement = "臺東")
family$site_id <- gsub(x = family$site_id, # 清除空格
                       pattern = "　", 
                       replacement = "")
names(family)[1] <- "縣市行政區"

family$household_no <- as.character(family$household_no) %>%
                       as.numeric() # factor型資料要轉成numeric型資料時，記得要先轉成character型資料
family <- group_by(family, 縣市行政區) %>% 
          summarise(household = sum(household_no)) # 計算縣市行政區的平均戶數
```

將地區戶數資料和train set、submit set結合
```{r}
train <- left_join(train, family, by = "縣市行政區")
```

檢查遺失值
```{r}
train$縣市行政區[is.na(train$household) == TRUE] %>% 
  as.factor %>%  
  table()
```

補遺失值
```{r}
# 屏東縣霧臺鄉有1049戶
train$household[train$縣市行政區 == "屏東縣霧臺鄉"] <- rep(1049, 6)
# 雲林縣臺西鄉有8727戶
train$household[train$縣市行政區 == "雲林縣臺西鄉"] <- rep(8727, 15)
# 高雄市三民區有134958戶
train$household[train$縣市行政區 == "高雄市三民區"] <- rep(134958, 86)
# 高雄市鳳山區有134958戶
train$household[train$縣市行政區 == "高雄市鳳山區"] <- rep(138016, 76)
```

將戶數資料和電桿資料併到submit set
```{r}
submit <- left_join(submit, train[, c(3, 14:24)], by = "VilCode")
```

## 加入颱風風力資料
讀取颱風風力資料，資料來源為[颱風資料庫](http://rdc28.cwb.gov.tw/)
```{r}
gust_csv <- read.csv("./data/gust.csv")
names(gust)[1] <- "CityName"
gust$CityName <- as.factor(gust$CityName)
```

將train set和submit set中的颱風資料選出來
```{r}
soudelor <- select(train, c(1:4, 13:24, 8))
meranti <- select(train, c(1:4, 13:24, 12))
megi <- select(submit, -c(5:6))
nesatAndHaitang <- select(submit, -c(5:6))
```

將颱風風力資料和train set和submit set中的颱風資料結合
```{r}
soudelor <- left_join(soudelor, 
                      gust[,c(1:3)], 
                      by = "CityName")
megi <- left_join(megi, 
                  gust[,c(1, 6:7)], 
                  by = "CityName")
meranti <- left_join(meranti, 
                     gust[,c(1, 12:13)], 
                     by = "CityName")
nesatAndHaitang <- left_join(nesatAndHaitang, 
                             gust[,c(1, 14:15)], 
                             by = "CityName")
```

## 建立隨機森林模型
用蘇迪勒颱風預測梅姬颱風
```{r}
library(randomForest)
names(soudelor)[18:19] <- c("maxWind", "gust")
names(megi)[16:17] <- c("maxWind", "gust")
soudelor_rf <- randomForest(Soudelor~., data = soudelor[, -c(1:5)])
soudelor_pred <- predict(soudelor_rf, newdata = megi[5:17])
megi_pred <- 1.45*soudelor_pred
```

用 莫蘭蒂+馬勒卡颱風 預測 尼莎+海棠颱風
```{r}
names(meranti)[18:19] <- c("maxWind", "gust")
names(nesatAndHaitang)[16:17] <- c("maxWind", "gust")
meranti_rf <- randomForest(MerantiAndMalakas~., data = meranti[, -c(1:5)])
meranti_pred <- predict(meranti_rf, newdata = nesatAndHaitang[5:17])
nesatAndHaitang_pred <- 1.53*meranti_pred
```

組合上傳檔案
```{r}
submit_dc <- cbind(submit[1:4], nesatAndHaitang_pred) %>% 
             cbind(megi_pred)
names(submit_dc)[5:6] <- c("NesatAndHaitang", "Megi")
write.csv(submit_dc, file = "submit_dc.csv", row.names = FALSE)
```