Do the delanay triangulation
library(deldir)
## deldir 0.1-15
# Generate points
x <- rnorm(500, 0, 1.5)
y <- rnorm(500, 0, 1)
vtess <- deldir(x, y)
##
## PLEASE NOTE: The components "delsgs" and "summary" of the
## object returned by deldir() are now DATA FRAMES rather than
## matrices (as they were prior to release 0.0-18).
## See help("deldir").
##
## PLEASE NOTE: The process that deldir() uses for determining
## duplicated points has changed from that used in version
## 0.0-9 of this package (and previously). See help("deldir").
# Delaunay triangulation
plot(x, y, type="n", asp=1)
plot(vtess, wlines="triang", wpoints="none", number=FALSE, add=TRUE, lty=1)
points(x, y, pch=20, col="black", cex=0.5)
# Get the things
#vtess$delsgs # x1, y1, x2, y2, index 1, index 2
estimate.point <- function(x, y, t) {
}
line.has.points <- function(i, j, line) {
(line[5]==i && line[6]==j) || (line[5]==j && line[6]==i)
}
# Finding which triangle a point is in
# (from https://stackoverflow.com/questions/2049582/how-to-determine-if-a-point-is-in-a-2d-triangle)
which.side.of.line <- function(p, p1, p2) {
(p$x - p2$x) * (p1$y - p2$y) - (p1$x - p2$x) * (p$y - p2$y);
}
is.inside.triangle <- function(p, p1, p2, p3) {
b1 <- which.side.of.line(p, p1, p2) < 0;
b2 <- which.side.of.line(p, p2, p3) < 0;
b3 <- which.side.of.line(p, p3, p1) < 0;
((b1 == b2) && (b2 == b3))
}
get.cell <- function(p, points, lines) {
for (i in 1:length(points)) {
for (j in (i+1):length(points)) {
# Check if points i and j are connected
skip <- TRUE
for (l in 1:length(lines)) {
if (line.has.points(i, j, lines[l, ])) {
skip <- FALSE
break
}
}
if (skip) next
for (k in (j+1):length(points)) {
# Check if points i, j, k forms a triangle
skip.1 <- TRUE
skip.2 <- TRUE
for (l in 1:length(lines)) {
if (line.has.points(i, k, lines[l, ]))
skip.1 <- FALSE
if (line.has.points(j, k, lines[l, ]))
skip.2 <- FALSE
}
if (skip.1 || skip.2) next
if (is.inside.triangle(p, points[i], points[j], points[k]))
return( c(i, j, k) )
}
}
}
}
# Mesh interpolation
get.triangle.area <- function(p1, p2, p3) {
abs(det(matrix(c(p1$x, p2$x, p3$x, p1$y, p2$y, p3$y, 1, 1, 1), nrow=3, ncol=3)))/2
}
mesh.interpolate <- function(ps) {
num <- 0
den <- 0
# ps <- c(p1, p2, p3, p4)
for (i in 1:3) {
V.i <- do.call(get.triangle.area, as.list(ps[-i]) )
num <- num + V.i * ps[[i]]$w
den <- den + V.i
}
num / den
}
# p1 <- list(x=0, y=0, w=1)
# p2 <- list(x=2, y=0, w=2)
# p3 <- list(x=0, y=4, w=3)
# p4 <- list(x=1, y=1)
# ps <- list(p1, p2, p3, p4)
# mesh.interpolate(ps)
library(ggplot2)
load("with_ecan.RData")
#ids <- c(105, 115, 114, 109) # last one is in middle # replaced 104 with 114
#sites <- c(1, 2, 8, 18) # changed 13 to 8
ids <- c(102, 105, 113, 109)
sites <- c(10, 1, 9, 18)
data <- data.frame(date = with.ecan$date, wind.speed = with.ecan$wind.speed)
for (i in 1:4) {
id <- ids[i]
site <- sites[i]
train.data <- data.frame(x = with.ecan[, paste0('pm2.5.odin.', id, '.site.18')],
y = with.ecan[, 'pm2.5'])
linear.model <- lm(y ~ x, train.data, na.rm=TRUE)
deployed.data <- data.frame(date = with.ecan$date,
x = with.ecan[, paste0('pm2.5.odin.', id, '.site.', site)])
results <- predict(linear.model, deployed.data)
deployed.data[, paste0('pm2.5.odin.', id)] <- predict(linear.model, deployed.data)
deployed.data$x <- NULL
data <- merge(data, deployed.data, by='date')
}
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
## extra argument 'na.rm' will be disregarded
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
## extra argument 'na.rm' will be disregarded
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
## extra argument 'na.rm' will be disregarded
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
## extra argument 'na.rm' will be disregarded
# Remove rows with NAs
data <- na.omit(data)
experiment <- function(data) {
ps <- list(NA,NA,NA,NA)
i <- 1
for (i in 1:4) {
p <- list()
p.id <- ids[i]
p.site <- sites[i]
p.x <- locs[which(locs$id==p.site), 'lon']
p.y <- locs[which(locs$id==p.site), 'lat']
p.w <- data[, paste0('pm2.5.odin.', p.id)]
p <- list(x=p.x, y=p.y, w=p.w)
ps[[i]] <- p
}
y.hat <- mesh.interpolate(ps)
y <- data[, paste0('pm2.5.odin.', p.id)]
mse <- mean( (y.hat - y)^2 )
mrae <- mean( abs(y.hat - y)/y )
print('date')
print(data$date[1])
print('mse')
print(mse)
print('mrae')
print(mrae)
}
experiment(data)
## [1] "date"
## [1] "2016-07-26 00:30:00 +12"
## [1] "mse"
## [1] 59.40089
## [1] "mrae"
## [1] 0.6131045
date.lt <- as.POSIXlt(data$date)
# 2016-08-09 gives na
start.date <- as.POSIXct('2016-08-28 12:00', format='%Y-%m-%d %H:%M', tz='Etc/GMT-12')
end.date <- as.POSIXct('2016-08-29 12:00', format='%Y-%m-%d %H:%M', tz='Etc/GMT-12')
low.wind.1 <- data[which(start.date < data$date & data$date < end.date), ]
experiment(low.wind.1)
## [1] "date"
## [1] "2016-08-28 12:10:00 +12"
## [1] "mse"
## [1] 289.8246
## [1] "mrae"
## [1] 0.4857813
start.date <- as.POSIXct('2016-08-14 12:00', format='%Y-%m-%d %H:%M', tz='Etc/GMT-12')
end.date <- as.POSIXct('2016-08-15 12:00', format='%Y-%m-%d %H:%M', tz='Etc/GMT-12')
low.wind.2<- data[which(start.date < data$date & data$date < end.date), ]
experiment(low.wind.2)
## [1] "date"
## [1] "2016-08-14 12:10:00 +12"
## [1] "mse"
## [1] 266.7016
## [1] "mrae"
## [1] 0.4229557
start.date <- as.POSIXct('2016-09-11 12:00', format='%Y-%m-%d %H:%M', tz='Etc/GMT-12')
end.date <- as.POSIXct('2016-09-12 12:00', format='%Y-%m-%d %H:%M', tz='Etc/GMT-12')
low.wind.3 <- data[which(start.date < data$date & data$date < end.date), ]
experiment(low.wind.3)
## [1] "date"
## [1] "2016-09-11 12:10:00 +12"
## [1] "mse"
## [1] 24.49705
## [1] "mrae"
## [1] 0.3911135
experiment.2 <- function(data) {
ps <- list(NA,NA,NA,NA)
i <- 1
for (i in 1:4) {
p <- list()
p.id <- ids[i]
p.site <- sites[i]
p.x <- locs[which(locs$id==p.site), 'lon']
p.y <- locs[which(locs$id==p.site), 'lat']
p.w <- data[, paste0('pm2.5.odin.', p.id)]
p <- list(x=p.x, y=p.y, w=p.w)
ps[[i]] <- p
}
y.hat <- mesh.interpolate(ps)
y <- data[, paste0('pm2.5.odin.', p.id)]
mse <- mean( (y.hat - y)^2 )
rae <- abs((y.hat - y) / y)
df.1 <- data.frame(date=data$date, data=rae, lab=rep('rae',length(data$date)))
df.2 <- data.frame(date=data$date, data=data$wind.speed, lab=rep('wind',length(data$date)))
df <- rbind(df.1, df.2)
plt <- ggplot(df) +
geom_line(aes(date, data,colour=lab))
print(plt)
mrae <- mean( rae )
}
experiment.2(low.wind.1)
experiment.2(low.wind.2)
experiment.2(low.wind.3)
