5

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.

Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
bbiasi
  • 1,549
  • 2
  • 15
  • 31

2 Answers2

7

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
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • Hello . Can you tell me why the results are different from the @www code? – bbiasi Aug 09 '18 at 16:25
  • 1
    @bbiasi 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. – Zheyuan Li Aug 09 '18 at 16:59
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
www
  • 38,575
  • 12
  • 48
  • 84
  • Why are the values obtained with the answers different? Using the code provided by @李哲源, I have that pH = 8.15849360036939. – bbiasi Aug 09 '18 at 16:23
  • 1
    @bbiasi I don't know but I think 李哲源's explanation is reasonable. You should consider increasing more samples around your prediction point. – www Aug 09 '18 at 16:45