WaterHackerthon

台灣超抽地下水潛在分析

資料集使用：

區域用戶用水量：http://www.water.gov.tw/ct.aspx?xItem=153656&CtNode=3394&mp=1

地層下陷磁層分層監測井監測資訊：http://data.gov.tw/node/41572

地層下陷磁層分層監測井基本資訊：http://data.gov.tw/node/41571

台灣電力公司_表燈用電不及底度用戶數：http://data.gov.tw/node/33463

雲林一級發布區：http://data.gov.tw/node/20597

地層下陷

magneticbasic <- as.data.frame(magneticbasic$WaterResourcesAgencyContinuousLandSubsidenceMagneticRingStationProfile_OPENDATA)
magnetic <- left_join(magneticnumber,magneticbasic,by = "LandSubsidenceMonitoringWellIdentifier")
magnetic <- magnetic[,c(-5,-8,-11,-12,-14,-15)]
magnetic <- magnetic[complete.cases(magnetic),]
magnetic$ObservationtTime <- strsplit(magnetic$ObservationtTime,
                                      split = "T") %>%
    lapply("[",1) %>% unlist() %>%
    as.Date(format = "%Y-%m-%d")
magnetic$SetTime <- strsplit(magnetic$SetTime,split = "T") %>%
    lapply("[", 1) %>% unlist() %>%
    as.Date(format = "%Y-%m-%d")
magnetic <- magnetic[magnetic$ObservationtTime >= as.Date(c("2011-01-04")),] %>%
    filter(StratumCompressionValue != "")
magnetic[,c(2,4,9)] <- sapply(magnetic[,c(2,4,9)], as.numeric)
yunlinmagnetic <- filter(magnetic,CountyName == "雲林縣")
yunlinmagneticmax <- group_by(yunlinmagnetic,LandSubsidenceMonitoringWellIdentifier) %>% 
    summarise(MagneticRingNumber = max(MagneticRingNumber)) %>% 
    left_join(y = yunlinmagnetic,by = c("MagneticRingNumber","LandSubsidenceMonitoringWellIdentifier"))
yunlinmagneticmax %<>% arrange(desc(StratumCompressionValue))
x <- filter(yunlinmagneticmax, TownName == "土庫鎮")

雲林

這裡從雲林挑出地層壓縮量最大的地區－土庫鎮

地層壓縮量在2014年底至2015年中有明顯上升的趨勢

造成地層壓縮量如此巨大變動的原因：

1.火山爆發：無

2.地震：2014年底至2015年中並無明顯的大地震

3.地層嚴重下陷：在2014年底至2015年中台灣發生了大乾旱，造成大量的農作停灌、工業用水減供，而尋求其他的取水來源？地下水！

推測：超抽地下水

說到地下水一定要鑿的地下水井，地下水井分為淺水井與深水井，目前研究都趨向於深水井是造成地層下陷的主因，淺水井造價便宜，而深水井動輒數百萬元，一般農民不可能會有錢去鑿深水井，當然選造價便宜的淺水井，能夠鑿的起深水井的，只剩下工廠、工業等等大型機關了。

g1 <- ggplot(x[x$ObservationtTime >= as.Date(c("2014-04-01")),][1:21,],
       aes(x = as.character(ObservationtTime),y = StratumCompressionValue)) +
    geom_point() + 
    geom_line(aes(group = 1)) +
    labs(x = "Year(2014~2015)") +
    ggtitle(label = "MW-HLES") +
    theme(plot.title = element_text(size = 25),
          axis.text = element_blank())
g2 <- ggplot(x[x$ObservationtTime >= as.Date(c("2014-04-01")),][22:41,],
       aes(x = as.character(ObservationtTime),y = StratumCompressionValue)) +
    geom_point() + 
    geom_line(aes(group = 1)) +
    labs(x = "Year(2011~2015)") +
    ggtitle(label = "MW-SETS") +
    theme(plot.title = element_text(size = 25),
          axis.text = element_blank())
grid.arrange(g1,g2,ncol = 2)

yunlinemap <- readOGR(dsn = "/Users/kangwenwei/Desktop/黑客入寶山-台灣自來水公司/1-10009/",
                 layer = "G97_10009_U0201_2015",encoding = "BIG-5",stringsAsFactors = FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/kangwenwei/Desktop/黑客入寶山-台灣自來水公司/1-10009/", layer: "G97_10009_U0201_2015"
## with 2160 features
## It has 11 fields
## Integer64 fields read as strings:  U_ID

## Warning in readOGR(dsn = "/Users/kangwenwei/Desktop/黑客入寶山-台灣自來水公
## 司/1-10009/", : Z-dimension discarded

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

de <- de[,-1]
colnames(de) <- c("郵遞區號","行政區",201601:201612)
z <- filter(waterconsumption,一級發布區代碼 %in% yunlinemap[yunlinemap$TOWN == "土庫鎮",]$CODE1 &
                用水種別 == 3)
tocoowater3dsum <- cbind(tFrame(filter(de,行政區 == "土庫")[-1:-2]),z)
tocoowater3dsum <- tocoowater3dsum[,-3:-4]
for (i in c(1,3)) {
    tocoowater3dsum[,i] <- tocoowater3dsum[,i]/mean(tocoowater3dsum[,i])
}

抽絲剝繭

當知道了雲林有可能有工業等大型機關違法鑿井的可能，如何盡量縮小範圍呢？

線索：對工業工廠來說，用水與用電長年來應該是成正比的。若用電與用水沒有呈現共同的趨勢，就有可能是潛在超抽地下水的兇手！

範圍：從雲林的一級發布區與供水用戶用水量資料比對，可以抽出土庫鎮的一級發布區代碼，然後就能從用水量挑出土庫鎮的一級發布區用水量，然後再從台灣電力公司表燈用電去挑出土庫鎮的用電，再進行正規化！

ggplot(tocoowater3dsum,aes(x = as.character(收費年月),y = V1)) +
    geom_bar(data = tocoowater3dsum,aes(x = as.character(收費年月), y = 用水量),stat = "identity",alpha = 0.5,fill = "blue") +
    geom_point(colour = "yellow") + geom_line(colour = "yellow",group = 1) +
    labs(x = NULL,y = NULL) +
    ggtitle("土庫鎮(肉類處理及其製品)") +
    theme(plot.title = element_text(size = 30,face = "plain"))

whoweiwater3dsum <- cbind(tFrame(filter(de,行政區 == "虎尾")[-1:-2]),z)
whoweiwater3dsum <- whoweiwater3dsum[,-3:-4]
for (i in c(1,3)) {
    whoweiwater3dsum[,i] <- whoweiwater3dsum[,i]/mean(whoweiwater3dsum[,i])
}
ggplot(whoweiwater3dsum,aes(x = as.character(收費年月),y = V1)) +
    geom_bar(data = whoweiwater3dsum,aes(x = as.character(收費年月), y = 用水量),stat = "identity",alpha = 0.5,fill = "blue") +
    geom_point(colour = "yellow") + geom_line(colour = "yellow",group = 1) +
    labs(x = NULL,y = NULL) +
    ggtitle("虎尾鎮(紡織、食品製造、化學製品)") +
    theme(plot.title = element_text(size = 30,face = "plain"))

雲林土庫鎮

yunlinelecwaterN3 <- read.csv("/Users/kangwenwei/Desktop/黑客入寶山-台灣自來水公司/雲林非工業用電 總用水(105上半年).csv",
                       fileEncoding = "BIG-5",stringsAsFactors = FALSE)
tocoo <- yunlinemap[yunlinemap$TOWN == "土庫鎮",]
tocoo1 <- yunlinemap[yunlinemap$CODE1 == filter(yunlinelecwaterN3,area == "土庫鎮")$code1,]
yunlinemap %<>% fortify(region == "TOWN")

## Regions defined for each Polygons

tocoo %<>% fortify(region == "TOWN")

## Regions defined for each Polygons

tocoo1 %<>% fortify(region == "TOWN")

## Regions defined for each Polygons

雲林縣

ggplot(yunlinemap,aes(x = long,y = lat,group = group)) +
    geom_polygon(fill = NA,colour = "black") +
    geom_polygon(data = tocoo,aes(x = long,y = lat,group = group),
                 fill = NA,colour = "orange") +
    labs(x = NULL,y = NULL) +
    theme_bw() +
    theme(panel.grid = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          panel.border = element_blank(),
          plot.title = element_text(size = 20)) +
    geom_polygon(data = tocoo1,aes(x = long,y = lat,group = group),
                 fill = "red") +
    geom_point(aes(186638.940,2620511.630),size = 7,shape = 21,colour = "red")