library(gstat)
library(sp)
setwd("D:/2019SpringTA/lab1") # Need to change the directory.
so4.orig <- read.table("so4.dt")
so4x <- -1*(so4.orig[,4] + so4.orig[,5]/60 + so4.orig[, 6]/3600)
so4y <- so4.orig[,1] + so4.orig[,2]/60 + so4.orig[,3]/3600
so4 <- data.frame(x=so4x, y=so4y, z=so4.orig[, 7])
so4[1:10, ]
coordinates(so4)=~x+y

par(mfrow = c(2, 2))
library(maps)
map("state")
points(so4$x, so4$y, pch = 19, cex = 0.5)

library(akima)
interp.so4 <- interp(so4$x, so4$y, so4$z)
image(interp.so4, xlim = c(-126, -66), ylim = c(23, 52), col = gray((11:2)/12),
      xlab = "longitude", ylab = "latitude")

library(lattice)
fit <- lm(z ~ poly(x, y, degree = 2), data = so4)
so4.fit <- data.frame(x = so4$x, y = so4$y, z = fit$fitted)
s <- with(so4.fit, interp(x, y, z))
trellis.par.set("axis.line", list(col = NA, lty = 1, lwd = 1))
plot(with(s, wireframe(z, row.values = x, col.values = y, xlab = "longitude", ylab = "latitude", 
                       zlab = "z", scales = list(arrows = FALSE), aspect = c(61/87, 0.5))),
     newpage = F, position = c(0, 0, 0.5, 0.6))

fit <- lm(z ~ poly(x, y, degree = 5), data = so4)
so4.fit <- data.frame(x = so4$x, y = so4$y, z = fit$fitted)
s <- with(so4.fit, interp(x, y, z))
trellis.par.set("axis.line", list(col = NA, lty = 1, lwd = 1))
plot(with(s, wireframe(z, row.values = x, col.values = y, xlab = "longitude", ylab = "latitude", 
                       zlab = "z", scales = list(arrows = FALSE), aspect = c(61/87, 0.5))), 
     newpage = F, position = c(0.5, 0, 1, 0.6))

par(mfrow = c(1, 1))