在 R 中绘制网格曲面图 [英] Plot a gridded surface plot in R

问题描述

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

I'm a beginner to R and I am trying to plot a surface plot on a specific grid. Basically I have a data-set of points from across the UK containing the longitude, latitude and amount of rainfall for a particular day. Using the following code I can plot this data onto a UK map:

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):

So far so good. I can also perform a thin plate spline interpolation (TPS) using:

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

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

and then I can do a surface plot at a grid scale of my choice e.g.:

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

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

This effectively gives me a gridded data plot at the resolution of 100*100.

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

Following help from another user I can now extract this data in a grid by using:

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 相同)再次绘制曲面图?你知道我该怎么做吗?

What I would like to do now is plot the surface plot again using the same grid as the predict function (i.e. same as xvals and yvals) as above? Do you know how I can do this?

感谢您的帮助

推荐答案

一旦您在 griddg 中预测了新值，您就可以在技术上使用 Tps 重新插值，并且然后像以前一样继续绘制表面图和地图:

Once you have predicted your new values in griddg, you can technically re-interpolate with Tps and then proceed with the surface plot and map as before:

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 两次.此方法需要一些外部函数，但您的映射将具有更大的灵活性.

For more control over your maps, you could also plot your new grid directly - This is probably more correct in that the above method essentially fits your interpolation Tps twice. This method requires some external functions, but you will have more flexibility in your mapping.

#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()

这篇关于在 R 中绘制网格曲面图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！