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

[Edzer Pebesma]

3D interpolation projected on irregular 2D is always somewhat tricky to
comprehend. After looking at the plot, I tried:

> library(sp)
> pts = SpatialPoints(mastCFD[c("x", "y")])
> grd = cfd
> gridded(grd) = ~x+y
> over(pts,grd)
         z   vxymag     aspZ  Ainterp
1 877.9992 11.70016 75491.32 4.700857
2 902.3500 13.44538 77585.03 4.669962
> mastCFD
       x       y      z                           mast     A     aspZ
1 586987 1566348 788.82 MASTBenakanahalli1_78_0m___Ch2 4.766 67823.60
2 591592 1562582 808.44 MASTBenakanahalli2_78_0m___Ch2 4.622 69510.55

and wondered why mastCFD$z is different from the cfd$z values at the
corresponding locations.


On 01/07/2015 07:31 AM, Frede Aakmann Tøgersen wrote:
> 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. 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma,  Co-Editor-in-Chief Computers & Geosciences
Institute for Geoinformatics (ifgi), University of Münster
Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20150107/489a5c36/attachment.bin>