農地重金屬污染潛勢分析實作 by R

緣起

土地孕育著無數生命，如同人類的母親。然而當台灣工商業愈加發達，許多農地也在發展歷程中遭汙染。

為讓民眾瞭解農地汙染狀況，台灣環境資訊協會2015年10月發起「守護農地計畫」，希望敦促政府積極行政，讓國民吃得安心，使土地生生不息。

透過智庫驅動（DSP）發起的「D4SG資料英雄計畫」，環資與政治大學資訊科學系/所、新聞系的學生2016年3月下旬組成跨領域團隊，透過資料科學方法，協助環資找出可能受重金屬汙染卻未受政府管制的農地，且分析各縣市的差異，而後希望發展相關專題報導。

這份教案在保留原汁原味的前提下，以環保署2016年的列管資料重現當時的研究成果，以本著作係採用創用 CC 姓名標示 3.0 台灣授權條款授權。

記事

2015.10. 台灣環境資訊協會發起守護農地計畫

2015.11. 智庫驅動與開拓基金會舉辦公益加值資料工作坊

2016.03. 參與D4SG資料英雄計畫，透過資料解決問題

2016.06. 成果發表會、環資發布新聞稿

2016.07. 成果納入環保署土壤及地下水污染調查及整治計畫

2017.02. 完成8百多處初篩，部分案例進入第二階段細調

2017.11. 將成果製作成R語言教材對外發布

2018.03. 預計完成細調並提報告

環境設定

1.載入實作練習包

先確認您的電腦安裝了R 與 R studio
進入課程練習包教材 Github頁面
將請到本頁面中點選 Clone or download -> Download ZIP，將課程包下載到自己的電腦中
解開壓縮後點擊 FPP.Rproj 開啟實作練習教材

2.安裝需到的R套件

list.of.packages <- c("knitr", "kableExtra", "stringr", "dplyr", "plotly", "scales", "leaflet")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

3.資料匯入

使用 read.csv 匯入資料
建議把游標移到" " (雙引號) 裡面按Tab鍵選取檔案

data_EPA <- read.csv("data/環保署列管污染農地_utf8.csv", fileEncoding = "utf8")
data_TARI <- read.csv("data/農地重金屬超標未列管_utf8.csv", fileEncoding = "utf8")

PART1–環保署資料實作

環保署「土壤及地下水列管場址」資料，提供有公告列管的農地污染控制場址地號，截至2016-12-18共計5485筆。

1.使用 `dplyr` 套件製作敘述統計表

library(dplyr)
#因為沒有全部要用，擷取需要用到的資訊，並另外命名為英文
data_EPA <- data_EPA %>% 
  select(county=1, coordinate=3, area=7, control_date=8, free_date=10)

# 預覽資料
head(data_EPA)

  county            coordinate    area control_date  free_date
1 臺北市 X：301277，Y：2777431  387.00   2002/10/04 2004/11/02
2 臺北市 X：301267，Y：2777372 7738.00   2002/10/04 2004/11/02
3 臺北市 X：301310，Y：2777998  580.00   2004/01/27 2005/03/01
4 臺北市 X：299512，Y：2779782 3205.65   2004/12/10 2008/04/18
5 臺北市 X：299518，Y：2779865 2855.14   2004/12/10 2008/04/18
6 臺北市 X：299669，Y：2779961 3513.00   2004/12/10 2008/04/18

summarise_data <- data_EPA %>% 
group_by(county) %>%                                     # 以縣市做加總
  summarise(count = n(),                                 # 所有案件數
            sum_area = sum(area),                        # 場址面積總和
            control_area = sum(area[free_date=="無"]),   # 列管面積總和 (平方公尺)
            free_area = sum(area[free_date!="無"]),      # 解除面積總和 (平方公尺)
            count_control = sum(free_date=="無"),        # 列管案件
            count_free = sum(free_date!="無")) %>%       # 解除列管案件
  mutate(avg_area = sum_area/count,                      # 平均場址面積 (平方公尺)
         ratio_control = sprintf("%1.0f%%",count_control/count*100), # 列管比例
         ratio_free = sprintf("%1.0f%%",count_free/count*100))       # 解除列管比例

作圖

library(plotly)
plot_ly(summarise_data, x = ~county, y = ~count_control, type = "bar") %>% 
  layout(title = "所有列管案件",
         xaxis = list(title = '全台縣市'),
         yaxis = list(title = '列管案件數'),
         width = 750, height = 400)

2.計算平均列管月份

#data_EPA的格式是Factor

data_EPA <- mutate(data_EPA,
                   control_date = as.Date(control_date), # 把列管時間轉成日期格式
                   free_date = as.Date(free_date), # 把列管時間轉成日期格式
                   # 把尚未解除列管的時間帶資料收集截止時間"2016-12-18"，以便計算列管時間
                   free_date = replace(free_date, 
                                       which(is.na(free_date)), 
                                       as.Date("2016-12-18"))
                   )

# 計算列管月份
data_EPA$totol_month <- NA
for(i in 1:nrow(data_EPA)){
  tryCatch({ # 因為資料裡有其中一筆沒有列管時間，我們加tryCatch讓他忽略錯誤，讓迴圈可以正常執行完
  data_EPA$totol_month[i] <- 
    length(seq(from=data_EPA$control_date[i], to=data_EPA$free_date[i], by='month')) 
  }, error=function(e){})
}

# 計算各縣市列管中與解除列管的平均列管月份
avg_month <- data_EPA %>% group_by(county) %>% 
  summarise(control_month = mean(totol_month[free_date=="2016-12-18"], na.rm = TRUE) %>% 
                            round(), 
  #找出列管中
            free_month = mean(totol_month[free_date!="2016-12-18"], na.rm = TRUE) %>% 
                         round())
  #找出解列的 
# 因為資料中有一筆沒有列管時間，因此把參數 na.rm 改為TRUE，才可以運算

# 把無列管中的平均月份帶0
avg_month[is.na(avg_month)] <- 0

# 合併回summarise_data
summarise_data <- summarise_data %>% 
  left_join(avg_month, by="county")

3.計算全國加總

tmp1 <- colSums(summarise_data[,2:8]) %>% t() %>% as.data.frame() %>% 
  mutate(ratio_control = sprintf("%1.0f%%",count_control/count*100),
         ratio_free = sprintf("%1.0f%%",count_free/count*100))
tmp2<-
  colMeans(summarise_data[,11:12]) %>% round() %>% t() %>% as.data.frame() %>%
  mutate(county="全國") 
#各縣市平均列管月數及平均解列月數轉換成全國平均列管月及解列月數
tmp<-cbind(tmp1,tmp2)
#將 全國相關的兩個資料做合併
summarise_data <- rbind(tmp, summarise_data)
#與原先各縣市資料做合併

作圖

Bar Chart

# x軸依據案件數做排序
plot_ly(summarise_data, 
         x = ~reorder(county, -count_control), 
         y = ~count_control, type = "bar") %>% 
  layout(title = "所有列管案件",
         xaxis = list(title = '全台縣市'),
         yaxis = list(title = '列管案件數'),
         width = 750, height = 400)

4.將欄位改成中文敘述並進行展示

tmp3 <- 
  summarise_data %>% 
  select(`全台縣市`=county,
         `案件數`=count,
         `場址面積總和`=sum_area,
         `列管面積總和`=control_area,
         `解除列管面積總和`=free_area,
         `平均場址面積`=avg_area,
         `列管案件`=count_control,
         `解除列管案件`=count_free,
         `列管案件比`=ratio_control,
         `解除列管案件比`=ratio_free,
         `列管中_平均列管月份`=control_month,
         `解除列管_平均列管月份`=free_month) %>% 
  arrange(desc(`案件數`))
tmp3

全台縣市	案件數	場址面積總和	列管面積總和	解除列管面積總和	平均場址面積	列管案件	解除列管案件	列管案件比	解除列管案件比	列管中_平均列管月份	解除列管_平均列管月份
全國	5485	9263582.55	3576260.23	5687322.33	51145.5310	2468	3017	45%	55%	40	35
彰化縣	2503	5038536.34	2282395.32	2756141.02	2012.9989	1340	1163	54%	46%	25	28
桃園市	1890	2511211.33	1216032.30	1295179.03	1328.6832	1060	830	56%	44%	48	42
臺中市	600	742225.90	10292.02	731933.88	1237.0432	11	589	2%	98%	32	30
新竹市	202	362010.73	2746.50	359264.23	1792.1323	2	200	1%	99%	7	27
臺南市	104	194793.93	34659.06	160134.87	1873.0186	37	67	36%	64%	105	30
高雄市	49	84948.44	0.00	84948.44	1733.6416	0	49	0%	100%	0	32
苗栗縣	37	43766.63	5729.23	38037.40	1182.8819	5	32	14%	86%	12	27
雲林縣	25	58230.54	4329.00	53901.54	2329.2216	2	23	8%	92%	78	31
臺北市	22	48852.15	0.00	48852.15	2220.5523	0	22	0%	100%	0	37
嘉義市	19	46035.64	12938.00	33097.64	2422.9284	5	14	26%	74%	28	49
新北市	13	37235.00	0.00	37235.00	2864.2308	0	13	0%	100%	0	22
南投縣	11	4046.00	2747.00	1299.00	367.8182	4	7	36%	64%	167	103
宜蘭縣	5	11798.15	3156.64	8641.51	2359.6300	1	4	20%	80%	98	22
屏東縣	3	75150.81	1235.16	73915.65	25050.2700	1	2	33%	67%	6	22
嘉義縣	2	4740.96	0.00	4740.96	2370.4800	0	2	0%	100%	0	16

PART2–農試所資料實作

農委會農業試驗所（農試所）之「土壤品質及生產力調查」資料。原始資料是1992年至2008年間進行的全國土壤採樣調查資料，累積約13萬筆，涵蓋78萬公頃表土，包含鎘、鉻、銅、鎳、鉛、鋅等六項重金屬有效性之調查。該資料是以250公尺*250公尺網格（6.25公頃）為單位。

由於農試所資料與環保署現今重金屬管制標準檢測方法不同，前者是透過0.1M鹽酸方法，後者是經過王水消化法，需要進行公式轉換，在此提供933筆環保數超標數據。

# 預覽資料
head(data_TARI)

         X        Y SAMELO OBJECTID  AREA PERIMETER PHTW911_ PHTW911_ID    XLO     YLO         SNO
1 120.1952 22.92694      0    37374 62500      1000    29807      33253 166625 2536625 GR85EF28521
2 120.2026 22.90666      0   110116 62500      1000       45         45 167375 2534375 GR95EF38841
3 120.2085 23.15728      0    32255 62500      1000    24435       3773 168125 2562125 GR83EG39041
4 120.2117 23.44625      0    66322 62500      1000    21553       6566 168625 2594125 GR91EH292A1
5 120.2172 22.91351      0    37673 62500      1000    30119      33631 168875 2535125 GR85EF39471
6 120.2190 23.02639      0    35537 62500      1000    27871       8128 169125 2547625 GR82EG89411
  YEAR    LNO      FE     MN     CU     ZN   CD    CR    NI     PB
1 1996 EF2852  606.90  58.19  35.40 280.98 0.86  0.38  6.05  75.95
2 2006 EF3884  617.95  50.72  26.54 243.58 1.57  1.02 11.63  40.78
3 1994 EG3904  909.50 151.41 165.91 197.95 1.01  1.22  3.93   9.88
4 2002 EH292A  449.21  42.09 164.63  13.50 0.03  0.17  1.40  17.71
5 1996 EF3947  315.54 121.34 413.09  48.83 1.27  0.45  1.94 126.43
6 1993 EG8941 1175.47  81.93 115.91  40.59 0.50 11.20 23.90  12.91
                                   add county
1    717台灣台南市仁德區中正西路1238號 台南市
2   829台灣高雄市湖內區中正路二段397巷 高雄市
3            723台灣台南市西港區金砂橋 台南市
4   Unnamed Road, 朴子市嘉義縣台灣 613 嘉義縣
5 717台灣台南市仁德區行大一街265-269號 台南市
6          710台灣台南市永康區河堤道路 台南市

1.重金屬轉換公式

農試所以 0.1M 鹽酸萃取抽出分析土壤中的重金屬濃度，其數據轉為王水消化法的推估值 (AR，以下為農試所使用的迴歸轉換公式) 後，再利用環保署現行管制標準進行篩檢。

重金屬名稱	重金屬轉換公式	環保署農地管制值
銅	CUAR = 2.035*CU + 11.884	200
鋅	ZNAR = 2.487*ZN + 89.711	600
鎘	CDAR = 1.4578*CD + 0.0323	5
鎘	CRAR = 17.35*CR+ 31.91	250
鎳	NIAR = 5.13*NI+ 14.56	200
鉛	PBAR = 2.811*PB+ 6.715	500

data_TARI <- data_TARI %>% 
  mutate(CUAR = 2.035*CU + 11.884,
         CUAR_OVER = ifelse(CUAR>200,1,0),
         ZNAR = 2.487*ZN + 89.711,
         ZNAR_OVER = ifelse(ZNAR>600,1,0),
         CDAR = 1.4578*CD + 0.0323,
         CDAR_OVER = ifelse(CDAR>5,1,0),
         CRAR = 17.35*CR+ 31.91,
         CRAR_OVER = ifelse(CRAR>250,1,0),
         NIAR = 5.13*NI+ 14.56,
         NIAR_OVER = ifelse(NIAR>200,1,0),
         PBAR = 2.811*PB+ 6.715,
         PBAR_OVER = ifelse(PBAR>500,1,0),
         OVER = ifelse(CUAR_OVER+ZNAR_OVER+CDAR_OVER+CRAR_OVER+
                       NIAR_OVER+PBAR_OVER>=1, 1, 0))

2 將農試所的超標與環保署列管案件在地圖上比較分布差異

2.1環保署資料場址座標轉經緯度

把TWD97位址資訊抓出來，分割成x軸與y軸

head(data_EPA,3) #可以看到座標欄位裡面有：和，

  county            coordinate area control_date  free_date totol_month
1 臺北市 X：301277，Y：2777431  387   2002-10-04 2004-11-02          25
2 臺北市 X：301267，Y：2777372 7738   2002-10-04 2004-11-02          25
3 臺北市 X：301310，Y：2777998  580   2004-01-27 2005-03-01          14

TWD97 <- data_EPA$coordinate %>% as.character() %>% 
#將原本是factor轉成character
  stringr::str_split(.,'[，：]',simplify = TRUE) %>% 
#將,;都切開
  as.data.frame(., stringsAsFactors=FALSE) %>% select(x=2,y=4) 
#只選第二及第四欄給欄位名x及y

將經緯度資訊從TWD97轉為WGS84
匯入轉碼程式包並用source()載入
利用轉碼函式將TWD97座標轉成經緯度座標系統

source("https://raw.githubusercontent.com/snexuz/TWD97TM2toWGS84/master/TWD97TM2toWGS84.R")
WGS84 <- TWD97TM2toWGS84(TWD97$x, TWD97$y) %>% as.data.frame()

#合併回原本的data
data_EPA$coordinate <- NULL
data_EPA <- cbind(data_EPA,TWD97,WGS84)
head(data_EPA)

  county    area control_date  free_date totol_month      x       y      lat      lon
1 臺北市  387.00   2002-10-04 2004-11-02          25 301277 2777431 25.10434 121.5084
2 臺北市 7738.00   2002-10-04 2004-11-02          25 301267 2777372 25.10381 121.5083
3 臺北市  580.00   2004-01-27 2005-03-01          14 301310 2777998 25.10946 121.5088
4 臺北市 3205.65   2004-12-10 2008-04-18          41 299512 2779782 25.12562 121.4910
5 臺北市 2855.14   2004-12-10 2008-04-18          41 299518 2779865 25.12637 121.4911
6 臺北市 3513.00   2004-12-10 2008-04-18          41 299669 2779961 25.12723 121.4926

2.2用leaflet畫地圖

地圖一：兩份資料交叉比對

library(leaflet)

#先把要點上圖的資料篩選出來
EPA <- data_EPA %>% filter(free_date=="2016-12-18") %>% select(lat, lon, area)
TARI <- data_TARI %>% select(lat = Y,lon = X, area = AREA)


leaflet() %>% setView(lng=120.58,lat=23.58,zoom=8) %>% 
  addProviderTiles("Esri.WorldImagery") %>%
  addCircles(data = EPA, color = "red", 
             lng = ~lon, lat = ~lat, weight = 1, radius = ~sqrt(area)/2) %>%
  addCircles(data = TARI, color = "yellow",
             lng = ~lon, lat = ~lat, weight = 1, radius = ~sqrt(area)/2)

Part 3 同學實作練習

1.使用 `plotly` 套件環保署繪製統計圖表

Grouped Bar Chart
with Hover Text and Rotated Labels

plot_ly(summarise_data, 
        x = ~reorder(county, -count_control), 
        y = ~count_control, type = 'bar',
        name = '列管案件', text = ~ratio_control) %>%
  add_trace(y = ~count_free, name = '解除列管案件', text = ~ratio_free) %>%
  layout(xaxis = list(title = "", tickangle = -45), 
         yaxis = list(title = 'count'), barmode = 'group',
         width = 800, height = 400)

Pie Chart

#列管縣市太多，很多面積很小，取場址面積大於中位數的來畫圖
summarise_data1 <- summarise_data %>% 
  arrange(sum_area %>% desc) %>%  # 降冪排序
  filter(sum_area > median(sum_area)) %>%  # 篩選大於中位數者
  slice(-1) # 移除全國總和數據

plot_ly(summarise_data1, labels = ~county, values = ~sum_area, type = 'pie',
        textinfo = 'label+percent') %>%
  layout(title = '場址面積總和',  width = 800, height = 500,
         margin=list(l = 100, r = 50, b = 100, t = 100, pad = 4))

2.計算農地樣區是否有超標情形

計算農試所資料與環保署案件資料場址距離
mindistance (與最近的環保署列管農地距離，單位公尺)
Is_monitored (境內是否有環保署列管農地。1有，0無)
minover (方圓1公里內其他超標點總數)

# 以農試所樣區資料為主體，計算每一個樣區與環保署列管農地的最短距離
for(i in 1:nrow(data_TARI)){
  data_TARI$mindistance[i] <- 
    ((data_TARI$XLO[i]-as.integer(data_EPA$x))^2+
       (data_TARI$YLO[i]-as.integer(data_EPA$y))^2) %>% 
    sqrt %>% min(na.rm = TRUE)
}


# 以農試所樣區為主體，計算每一個樣區方圓一公里內有多少個超標樣區
for(i in 1:nrow(data_TARI)){
  tmp <- sqrt((data_TARI$XLO[i]-data_TARI$XLO)^2+(data_TARI$YLO[i]-data_TARI$YLO)^2)
  
  data_TARI$minover[i] <- sum(tmp<=1000, na.rm = TRUE)
}

# 以農試所樣區為主體，計算每一個樣區境內是否有環保署列管農地
data_TARI <- data_TARI %>% mutate(Is_monitored=ifelse(mindistance<=200, 1, 0))

3.與環保署的資料合併，計算統計指標

library(scales)
# 計算各縣市統計 `group_by(county)`
stat_over <- 
  data_TARI %>% group_by(county) %>%
  summarise(`重金屬超標`=sum(OVER),
            `超標未曾列管`=sum(ifelse(OVER+Is_monitored==1,1,0))) %>%
  mutate(`超標未曾列管比例`=percent(`超標未曾列管`/`重金屬超標`),
         `超標列管比例`=percent(1-`超標未曾列管`/`重金屬超標`)) %>% 
  arrange(desc(`重金屬超標`))

# 加入全國訊息
stat_over <-
  data_TARI %>% 
  summarise(`重金屬超標`=sum(OVER),
            `超標未曾列管`=sum(ifelse(OVER+Is_monitored==1,1,0))) %>%
  mutate(`超標未曾列管比例`=percent(`超標未曾列管`/`重金屬超標`),
         `超標列管比例`=percent(1-`超標未曾列管`/`重金屬超標`)) %>% 
  mutate(county="全國") %>% 
  select(5, 1:4) %>%    # 使用mutate新增的欄位會排在最後面，利用select把county的順序往前移 
  full_join(stat_over)
  
# stat_over

county	重金屬超標	超標未曾列管	超標未曾列管比例	超標列管比例
全國	993	875	88.1%	11.9%
桃園市	411	361	87.8%	12.2%
彰化縣	159	98	61.6%	38.4%
屏東縣	86	86	100.0%	0.0%
高雄市	85	85	100.0%	0.0%
台中市	77	71	92.2%	7.8%
台南市	32	32	100.0%	0.0%
新竹縣	30	30	100.0%	0.0%
嘉義縣	23	23	100.0%	0.0%
宜蘭縣	22	22	100.0%	0.0%
苗栗縣	18	17	94.4%	5.6%
南投縣	13	13	100.0%	0.0%
新北市	11	11	100.0%	0.0%
雲林縣	10	10	100.0%	0.0%
花蓮縣	8	8	100.0%	0.0%
台東縣	5	5	100.0%	0.0%
嘉義市	2	2	100.0%	0.0%
台北市	1	1	100.0%	0.0%

4.用log(mindistance)畫地圖

TARI <- data_TARI %>% 
  select(lat = Y,lon = X, area = AREA ,mindistance) %>%
  mutate(log_mindistance=log(mindistance))

leaflet() %>% 
  addCircles(data=TARI,lng = ~lon, lat = ~lat, 
             radius = ~sqrt(area)/2, weight = 1,
             fill=TRUE, fillOpacity = 0.8,
             color=~colorNumeric(c("#CD0000", "#FFFFFF","#0B752F"),
                                 log_mindistance)(log_mindistance)) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addLegend(position = 'topleft',
            pal =colorNumeric(c("#CD0000", "#FFFFFF","#0B752F"),
                              domain=TARI$log_mindistance),
            values=TARI$log_mindistance,
            title = 'log-mindistance')

農地重金屬污染潛勢分析實作 by R

Ting-An Jiang, Yu Ying Han, Johnson Hsieh, Ning Chen

2017年11月25日

緣起

記事

環境設定

1.載入實作練習包

2.安裝需到的R套件

3.資料匯入

PART1–環保署資料實作

1.使用 `dplyr` 套件製作敘述統計表

作圖

2.計算平均列管月份

3.計算全國加總

作圖

4.將欄位改成中文敘述並進行展示

PART2–農試所資料實作

1.重金屬轉換公式

2 將農試所的超標與環保署列管案件在地圖上比較分布差異

2.1環保署資料場址座標轉經緯度

2.2用leaflet畫地圖

地圖一：兩份資料交叉比對

Part 3 同學實作練習

1.使用 `plotly` 套件環保署繪製統計圖表

2.計算農地樣區是否有超標情形

3.與環保署的資料合併，計算統計指標

4.用log(mindistance)畫地圖

農地重金屬污染潛勢分析實作 by R

Ting-An Jiang, Yu Ying Han, Johnson Hsieh, Ning Chen

2017年11月25日

緣起

記事

環境設定

1.載入實作練習包

2.安裝需到的R套件

3.資料匯入

PART1–環保署資料實作

1.使用 dplyr 套件製作敘述統計表

作圖

2.計算平均列管月份

3.計算全國加總

作圖

4.將欄位改成中文敘述並進行展示

PART2–農試所資料實作

1.重金屬轉換公式

2 將農試所的超標與環保署列管案件在地圖上比較分布差異

2.1環保署資料場址座標轉經緯度

2.2用leaflet畫地圖

地圖一：兩份資料交叉比對

Part 3 同學實作練習

1.使用 plotly 套件環保署繪製統計圖表

2.計算農地樣區是否有超標情形

3.與環保署的資料合併，計算統計指標

4.用log(mindistance)畫地圖

1.使用 `dplyr` 套件製作敘述統計表

1.使用 `plotly` 套件環保署繪製統計圖表