How can I predict values for a specific point using the idw() function in R?

Question

Using this answer from Ege Rubak as an example, how can I predict pH values for a specific point, say lat = -23.49184 and long = 152.07185, using the idw() function in R?

The closest answer I found was through this document in RPubs, but I could not extract just a specific value.

library(gstat)
library(sp)

lat <-  c(-23.49174, -23.49179, -23.49182, -23.49183, -23.49185, -23.49187)
long <- c(152.0718, 152.0718, 152.0717, 152.0717, 152.0717, 152.0717)
pH <- c(8.222411, 8.19931, 8.140428, 8.100752, 8.068141, 8.048852)
sample <- data.frame(lat, long, pH)

x.range <- range(sample$long)
y.range <- range(sample$lat)

x<-seq(x.range[1], x.range[2], length.out=20)
y<-seq(y.range[1], y.range[2], length.out=20)
grd<-expand.grid(x,y)

coordinates(sample) = ~long+lat
coordinates(grd) <- ~ Var1+Var2
gridded(grd) <- TRUE

proj4string(sample) <- CRS("+proj=longlat +datum=WGS84")
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84")

dat.idw <- idw(formula=pH ~ 1, locations = sample, newdata = grd, idp = 2.0)
#> [inverse distance weighted interpolation]

I did not specifically ask Ege Rubak in a comment because I do not yet have 50 reputations.

Absorb

Answer 1

You don't need a grid. Provide your new locations in a consistent way that your observed locations are represented.

library(gstat)
library(sp)

lat <-  c(-23.49174, -23.49179, -23.49182, -23.49183, -23.49185, -23.49187)
long <- c(152.0718, 152.0718, 152.0717, 152.0717, 152.0717, 152.0717)
pH <- c(8.222411, 8.19931, 8.140428, 8.100752, 8.068141, 8.048852)

sample <- data.frame(lat, long, pH)
coordinates(sample) = ~long+lat
proj4string(sample) <- CRS("+proj=longlat +datum=WGS84")

loc <- data.frame(long = 152.07185, lat = -23.49184)
coordinates(loc)  <- ~ long + lat
proj4string(loc) <- CRS("+proj=longlat +datum=WGS84")

oo <- idw(formula=pH ~ 1, locations = sample, newdata = loc, idp = 2.0)
oo@data$var1.pred
#[1] 8.158494

Counterpoise

Hello . Can you tell me why the results are different from the @www code? – Absorb

@Absorb I should be more precise. It is because your location does not lie exactly on that grid. Although the grid points are predicted with inverse distance method, locations elsewhere are probably interpolated with linear methods by projection. Your linked answer used a grid, because he want to use plot to make plot. – Counterpoise

Answer 2

You can use the extract function from the raster package. Notice that your point is outside of your original grid, so I increase by 1.5 to cover the point.

library(gstat)
library(sp)

lat <-  c(-23.49174, -23.49179, -23.49182, -23.49183, -23.49185, -23.49187)
long <- c(152.0718, 152.0718, 152.0717, 152.0717, 152.0717, 152.0717)
pH <- c(8.222411, 8.19931, 8.140428, 8.100752, 8.068141, 8.048852)
sample <- data.frame(lat, long, pH)

x.range <- range(sample$long)
y.range <- range(sample$lat)

x<-seq(x.range[1], x.range[2] * 1.5, length.out=20)
y<-seq(y.range[1], y.range[2] * 1.5, length.out=20)
grd<-expand.grid(x,y)

coordinates(sample) = ~long+lat
coordinates(grd) <- ~ Var1+Var2
gridded(grd) <- TRUE

proj4string(sample) <- CRS("+proj=longlat +datum=WGS84")
proj4string(grd) <- CRS("+proj=longlat +datum=WGS84")

dat.idw <- idw(formula=pH ~ 1, locations = sample, newdata = grd, idp = 2.0)


library(raster)

# Convert to raster
dat.r <- raster(dat.idw)

# Create a matrix showing the coordinate of interest
p <- SpatialPoints(matrix(c(152.07185, -23.49184), ncol = 2))
proj4string(p) <- projection(dat.r)

# Extract the values
extract(dat.r, p)
# 8.048852 

Hambrick

Why are the values obtained with the answers different? Using the code provided by @李哲源, I have that pH = 8.15849360036939. – Absorb

@Absorb I don't know but I think 李哲源's explanation is reasonable. You should consider increasing more samples around your prediction point. – Hambrick