library(sgeostat)
data(maas)

# grid <- list(x=seq(min(maas$x), max(maas$x), by=100),
#           y=seq(min(maas$y), max(maas$y), by=100))

grid <- list(x=seq(min(maas$x), max(maas$x), length=40),
           y=seq(min(maas$y), max(maas$y), length=40))

grid$xy <- data.frame(cbind(c(matrix(grid$x, length(grid$x), length(grid$x))),
             c(matrix(grid$y, length(grid$y), length(grid$y), byrow=T))))
                   
colnames(grid$xy) <- c("x", "y")

maas.lm1 <- lm(zinc ~ x + y, maas)
grid$lm1 <- predict(maas.lm1, data.frame(x=grid$xy$x, y=grid$xy$y))
image(grid$x, grid$y, matrix(grid$lm1, length(grid$x),length(grid$y)))
points(maas$x, maas$y, pch="x",col="black")
contour(grid$x, grid$y, matrix(grid$lm1, length(grid$x),
        length(grid$y)), add=T)

maas.lm2 <- lm(zinc ~ I(x^2) + I(y^2) + x + y + x:y, maas)
grid$lm2 <- predict(maas.lm2, data.frame(x=grid$xy$x, y=grid$xy$y))
image(grid$x, grid$y, matrix(grid$lm2, length(grid$x),length(grid$y)))
points(maas$x, maas$y, pch="x",col="black")
contour(grid$x, grid$y,
                matrix(grid$lm2, length(grid$x), length(grid$y)), add=T)