[R-sig-Geo] Inverse distance weighting in 3 dimension using gstat package

[Frede Aakmann Tøgersen]

Hi

Anyone got experiences with interpolation in 3 dimensions?

I am doing some inverse distance interpolation in 3 dimensions using the gstat package. The interpolation is based on observation at few locations (here only 2) and this is why I'm using IDW. The 3rd dimension is a scaled version of height above sea level. The parameter A physically increase by increasing height but not following any linear model.   An png file of an IDW interpolation of a parameter, A, to a 25 by 25 meter regular grid can be downloaded from https://www.dropbox.com/s/ilp4npzdhew48i6/Ainterp.png?dl=0. As you can see the interpolated values to a large extent follows the height contours of the region as wanted. But there are anomalies around the mast called "mast_1" in the upper left corner where the parameter A has higher values down-slope than at the top of the hill and further down-slope. Some sort of a ripple or wave effect.

Is this to be expected from IDW interpolation or am I doing something wrong? Your help is much appreciated.

The data can be downloaded from https://www.dropbox.com/s/ia2czqckdbjtla0/mastCFD.csv?dl=0  and https://www.dropbox.com/s/fpqt9745onlw7v0/cfd.csv?dl=0.  

Here is my R code for doing the interpolation.

### R script

library(lattice)
library(gstat)

mastCFD <- read.table("mastCFD.csv", h = TRUE, sep = ";", dec = ".")
cfd <- read.table("cfd.csv", h = TRUE, sep = ";", dec = ".")

mastCFD.gstat <- gstat(id = "A", formula = A ~ 1, locations = ~ x + y + aspZ, data = mastCFD, set = list(idp = 2))
mastCFDscale <- predict(mastCFD.gstat, cfd[,c("x","y", "aspZ")], debug.level = 0)

cfd[,"Ainterp"] <- mastCFDscale$A.pred

png("Ainterp.png", width = 1800, height = 900)

levelplot(Ainterp ~ x*y, data = cfd, aspect = myasp,
          panel = function(x,y,z,subscripts,...){
              panel.levelplot(x,y,z,subscripts,...)
              panel.contourplot(cfd$x, cfd$y, cfd$z, subscripts, region=FALSE,
                                contour=TRUE, at=pretty(cfd$z,10),
                                lwd=.3, labels = list(cex = 1))
              panel.points(mastCFD$x, mastCFD$y, pch=19,col ="black", cex=1.1*1.5,lwd=.3)
              panel.points(mastCFD$x, mastCFD$y, pch=19, col ="red",cex=0.4*1.5)
          }, col.regions = terrain.colors,
          ylab = "Northings (m)", xlab = "Eastings (m)")

dev.off()

Yours sincerely / Med venlig hilsen


Frede Aakmann Tøgersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender.