这篇文章是2014年10月21日~29日在台北参加“国际Open Data应用实务班”的学员Yu-Lang Chiang的报告。

项目

我们利用台湾省污染分布的开放数据和R相关的技术来可视化这些数据。下载幻灯片台灣公害資料視覺化

背景

TGOS(Taiwan Geospatial One Stop,地理資訊圖資雲服務平台) 是台湾的地理空间数据服务平台。你可以在这里查看关于该平台的详细介绍信息。因为我们来自一家空间地理信息技术的公司,上了SupStat开放数据课程之后,我们决定利用R对TGOS上面的2013年污染分布的数据做可视化分析。

数据清洗

尽管任何人都可以在TGOS的网站上非常容易地下载到数据,但是数据的质量参差不齐。例如,不同部门使用不同的坐标系给地理信息编码。数据的元数据经常是空的,没有提供关于数据模式(data schema)的足够信息。所以我们使用RGDAL包来做坐标变换的工作。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
library(rgdal)
setwd("C:\\temp\\opendataWorkDir")

# import original data
data <- system.file("data/102_data_demo.shp", package = "rgdal")

# read harzard data
harzard_s <- readOGR(dsn="data","102_data_demo",encoding="big5")
# read county line
county_s <- readOGR(dsn="data","city",encoding="big5")

# transfer to the WGS84 datum
EPSG <- make_EPSG()
EPSG[3826,3]

harzard.linglat <- spTransform(harzard_s, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
county.linglat <- spTransform(county_s, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

我们本来想联系提供污染分布数据的部门,以验证数据模式,但是没有收到回复。我们很幸运,回去跟同事们讨论之后,发现有位同事以前就在该部门的信息系统开发团队。最终我们找到数据模式,并安全地去除了没用的数据。

1
2
3
4
5
6
7
8
9
10
11
#use the over function to do spatial comparision
sc.linglat <-over(harzard.linglat, county.linglat)

# combine data
cc <- cbind(harzard.linglat@data,sc.linglat)

# remove repeated data
cc2 <- subset(cc,select=c(SerID,ALIAS,D_NAME,D_NAME1))

# merge with the original data
harzard_sc.linglat <- merge(harzard.linglat,cc2,by="SerID")

经过几次尝试之后,我们发现直接在R里用地理信息数据的原始格式绘制图形是非常慢的。我们把数据导出成rds格式就解决了这个问题。

1
2
3
# export to RDS format
saveRDS(harzard_sc.linglat, "data/harzard.rds")
saveRDS(county.linglat, "data/county.rds")

可视化

我们想用颜色的深浅来表示污染的数量。越深的颜色表示有越多的污染。我们用classIntRColorBrewer函数来完成这项工作:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# Load plyr library
library(plyr)

# Load specific data
subharz <- subset(harzard_sc,pollution1 == '水污染')
# Change to dataframe
dfsubharz <- as.data.frame(subharz)
# caculate the report number of differenct administrative distric
areacount <- ddply(dfsubharz,~D_NAME,summarise,count=length(unique(SerID)))


# color palette configuration
library(RColorBrewer)
library(classInt)
nclr <- 8
plotclr <- brewer.pal(nclr, "BuPu")

class <- classIntervals(areacount$count, nclr, style = "equal")

colcode <- findColours(class, plotclr)
tmpCounty <- merge(county,areacount,by.x="D_NAME",by.y="D_NAME.x")
# draw the map
spplot(tmpCounty,"count", col.regions  =  plotclr, at  =  round ( class $ brks ,digits  =  1 ))

例如,下面是台湾噪音污染的分布地图(点击可查看大图):

用shiny做分享

我们在这次课程学习了几种数据可视化的工具,比如shiny、D3、ECharts。最终我们选择用shiny来分享工作成果(点击可查看大图):

用户可以使用左侧的组合框来选择不同的污染类型来做可视化,比如噪音(noise)、垃圾(garbage)、振动(vibration)等。界面也可以在右侧显示详细的数据表。

感想

我们以前从没用过R,但是参加了这次课程之后,我们真心觉得R非常强大。R语言的在线社区也非常活跃。在网络上可以找到大量的工具、扩展包以及文档。如果未来我们面临统计学、大数据、可视化的问题,我们会尝试用R来解决。