コーセラのDeveloping Data Productで知りたい場所の緯度経度を入力すると過去の観測データからオゾン濃度とPM25汚染濃度を予測するコードを紹介していた。Video Lectureのみだったので、コーディングしてコメントを追加しておく。
サイトのURLはこちら https://class.coursera.org/devdataprod-002/
setwd("C:/r/yhat")
## 大気汚染のデータをRに取り込む。
d <- read.csv("annual_all_2013.csv")
str(d)
## PM2.5およびオゾンデータのみとりこむ。さらに必要な変数のみにする。
sub <- subset(d, (Parameter.Name %in% c("PM2.5 - Local Conditions", "Ozone")
& Pollutant.Standard %in% c("Ozone 8-Hour 2008", "PM25 Annual 2006")),
c(Longitude, Latitude, Parameter.Name, Arithmetic.Mean))
head(sub)
## 時系列で重複しているので地点で集計する。
pollavg <- aggregate(sub[,"Arithmetic.Mean"],
sub[,c("Longitude", "Latitude", "Parameter.Name")],
mean, na.rm = TRUE)
names(pollavg)[4] <- "level"
pollavg <- transform(pollavg, Parameter.Name = factor(Parameter.Name))
pollavg
rm(d, sub)
## 監視地点の行列を作成する。
monitors <- data.matrix(pollavg[,c("Longitude", "Latitude")])
monitors
library(fields)
## input parameter : data frame
## lon: longitude
## lat: latitude
## radius: Radius in miles for findng monitors.
pollutant <- function(df){
x <- data.matrix(df[, c("lon", "lat")])
r <- df$radius
## パラメータで引き渡された地点と監視地点の距離を計算する
d <- rdist.earth(monitors, x)
# dはx地点からのmonitors地点への距離が含まれる。
# dはリストで戻される、行数=観察対象の地点数、列数=モニタしている地点数
## モニタしている地点すべてについて 距離 < rとなる地点を探す
use <- lapply(seq_len(ncol(d)), function(i) {
which (d[,i] < r[i])
})
## オゾンとPM25について集計する
levels <- sapply(use, function(idx){
with(pollavg[idx,], tapply(level,Parameter.Name, mean))
})
dlevels <- as.data.frame(t(levels))
data.frame(df, dlevels)
}
## テスト
pollutant(data.frame(lon = -76.61, lat=39.2, radius=40))
