r - 在 R 中绘制网格曲面图

Question

我是 R 的初学者，我正在尝试在特定网格上绘制曲面图。基本上，我有一个来自英国各地的点数据集，其中包含特定日期的经度、纬度和降雨量。使用以下代码，我可以将此数据绘制到英国地图上：

dat <- read.table("~jan1.csv", header=T, sep=",")
names(dat) <- c("gauge", "date", "station", "mm", "lat", "lon", "location", "county",    "days")
library(fields)
quilt.plot(cbind(dat$lon,dat$lat),dat$mm)
world(add=TRUE)

到目前为止，一切都很好。我还可以使用以下方法执行薄板样条插值 (TPS)：

fit <- Tps(cbind(dat$lon, dat$lat), dat$mm, scale.type="unscaled")

然后我可以在我选择的网格比例上绘制曲面图，例如：

surface (fit, nx=100, ny=100)

这有效地为我提供了分辨率为 100*100 的网格数据图。

在另一个用户的帮助下，我现在可以使用以下方法在网格中提取这些数据：

xvals <- seq(-10, 4, len=20)
yvals <- seq(49, 63, len=20)
griddf <- expand.grid(xvals, yvals)
griddg <- predict(fit, x=as.matrix(griddf) )

我现在想做的是使用与上面的预测函数相同的网格（即与xvals和yvals相同）再次绘制曲面图？你知道我该怎么做吗？

谢谢你的帮助

Answer 1

一旦您预测了 中的新值griddg，您可以在技术上重新插值，Tps然后像以前一样继续进行曲面图和地图：

例子：

xvals <- seq(-10, 4, len=20)
yvals <- seq(49, 63, len=20)
griddf <- expand.grid(lon=xvals, lat=yvals)
griddg <- predict(fit, x=as.matrix(griddf) )

dat2 <- cbind(griddf, mm=griddg)
head(dat2)
fit <- Tps(cbind(dat2$lon, dat2$lat), dat2$mm, scale.type="unscaled")
surface (fit, nx=100, ny=100)
world(add=TRUE)

为了更好地控制您的地图，您还可以直接绘制新网格 - 这可能更正确，因为上述方法基本上适合您的插值Tps两次。此方法需要一些外部函数，但您的映射将具有更大的灵活性。

#option 2
source("matrix.poly.r") #http://menugget.blogspot.de/2012/04/create-polygons-from-matrix.html
source("val2col.R") # http://menugget.blogspot.de/2011/09/converting-values-to-color-levels.html
source("image.scale.R") # http://menugget.blogspot.de/2011/08/adding-scale-to-image-plot.html

#new grid and predition
xvals <- seq(-10, 4, len=100)
yvals <- seq(49, 63, len=100)
griddf <- expand.grid(lon=xvals, lat=yvals)
griddg <- predict(fit, x=as.matrix(griddf) )

#make polygons for new grid, calculate color levels
mat <- matrix(griddg, nrow=length(xvals), ncol=length(yvals))
poly <- matrix.poly(xvals, yvals, z=mat, n=seq(mat))
pal <- colorRampPalette(c("blue", "cyan", "yellow", "red"))
COL <- val2col(mat, col=pal(100))

#required packages
library(maps)
library(mapproj)

#plot
png("tmp.png", width=5, height=4, res=400, units="in")
layout(matrix(1:2, nrow=1, ncol=2), widths=c(4,1), heights=4)
par(mar=c(1,1,1,1))
map("world", proj="stereographic", orient=c(mean(yvals),mean(xvals),0), par=NULL, t="n", xlim=range(xvals), ylim=range(yvals))
for(i in seq(poly)){
 polygon(mapproject(poly[[i]]), col=COL[i], border=COL[i], lwd=0.3)
}
map("world", proj="stereographic", orient=c(mean(yvals),mean(xvals),0), par=NULL, add=T)
map.grid(col=rgb(0,0,0,0.5), labels=F)
box()

par(mar=c(5,0,5,4))
image.scale(mat, col=pal(100), horiz=FALSE, axes=FALSE, xlab="", ylab="")
axis(4)
mtext("mm", side=4, line=2.5)
box()

dev.off()