How can I predict values for a specific point usin

2019-03-01 06:13发布

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.

2条回答
劳资没心,怎么记你
2楼-- · 2019-03-01 06:32

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 
查看更多
你好瞎i
3楼-- · 2019-03-01 06:37

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
查看更多
登录 后发表回答