如何在稀疏点之间插入数据以在 R & plotly 中绘制等高线图
Posted
技术标签:
【中文标题】如何在稀疏点之间插入数据以在 R & plotly 中绘制等高线图【英文标题】:How to interpolate data between sparse points to make a contour plot in R & plotly 【发布时间】:2019-07-10 15:09:18 【问题描述】:我想根据第一个图中以下彩色点的浓度数据在 xy 平面上创建等高线图。我在每个高度都没有角点,所以我需要将浓度外推到 xy 平面的边缘 (xlim=c(0,335),ylim=c(0,426))。
点的情节 html 文件可在此处获得:https://leeds365-my.sharepoint.com/:u:/r/personal/cenmk_leeds_ac_uk/Documents/Documents/HECOIRA/Chamber%20CO2%20Experiments/Sensors.html?csf=1&e=HiX8fF
dput(df)
structure(list(Sensor = structure(c(11L, 12L, 13L, 14L, 15L,
16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L,
29L, 1L, 3L, 4L, 5L, 6L, 8L, 30L, 31L, 32L, 33L, 34L, 35L), .Label = c("N1",
"N2", "N3", "N4", "N5", "N6", "N7", "N8", "N9", "Control", "A1",
"A10", "A11", "A12", "A13", "A14", "A15", "A16", "A17", "A18",
"A19", "A2", "A3", "A4", "A5", "A6", "A7", "A8", "A9", "R1",
"R2", "R3", "R4", "R5", "R6"), class = "factor"), calCO2 = c(2237,
2389.5, 2226.5, 2321, 2101.5, 1830.5, 2418, 2356.5, 435, 2345.5,
2376, 2451, 2397, 2466, 2518.5, 2087, 2463, 2256.5, 2345.5, 3506,
2950, 3386, 2511, 2385, 3441, 2473, 2357.5, 2052.5, 2318, 1893.5,
2251), x = c(83.75, 167.5, 167.5, 167.5, 251.25, 167.5, 251.25,
251.25, 0, 83.75, 251.25, 167.5, 251.25, 83.75, 83.75, 83.75,
83.75, 251.25, 167.5, 335, 0, 0, 335, 167.5, 167.5, 167.5, 0,
335, 335, 167.5, 167.5), y = c(213, 319.5, 319.5, 110, 319.5,
213, 110, 110, 356, 213, 319.5, 110, 213, 110, 319.5, 319.5,
110, 213, 213, 0, 0, 426, 426, 426, 0, 213, 213, 70, 213, 426,
0), z = c(155, 50, 155, 155, 155, 226, 50, 155, 178, 50, 50,
50, 50, 155, 50, 155, 50, 155, 50, 0, 0, 0, 0, 0, 0, 0, 130,
50, 120, 130, 130), Type = c("Airnode", "Airnode", "Airnode",
"Airnode", "Airnode", "Airnode", "Airnode", "Airnode", "Airnode",
"Airnode", "Airnode", "Airnode", "Airnode", "Airnode", "Airnode",
"Airnode", "Airnode", "Airnode", "Airnode", "Naveed", "Naveed",
"Naveed", "Naveed", "Naveed", "Naveed", "Rotronic", "Rotronic",
"Rotronic", "Rotronic", "Rotronic", "Rotronic")), .Names = c("Sensor",
"calCO2", "x", "y", "z", "Type"), row.names = c(NA, -31L), class = "data.frame")
require(plotly)
plot_ly(data = subset(df,z==0), x=~x,y=~y, z=~calCO2, type = "contour") %>%
layout(
xaxis = list(range = c(340, 0), autorange = F, autorange="reversed"),
yaxis = list(range = c(0, 430)))
我正在尝试找到类似的东西。任何帮助将不胜感激。
【问题讨论】:
而不是打印head(df)
,您应该包含dput(df)
的输出。这将使您的数据框可重现
我了解@SeGa
我同意,我即将与您分享代码,但我想与您一起尝试完整的数据。你能通过云与我们分享吗?
嗨@CésarArquero,非常感谢你!数据集是上面的 df,我现在已经使用 dput(df) 来展示它的完整构造。
【参考方案1】:
首先,您必须考虑到 +-30 分不足以获得您在示例中看到的那些干净的分离层。说了这么多,开始工作吧:
首先,您可以监督您的数据,以猜测这些层的形状将如何。在这里您可以很容易地看到较低的 z 值具有较高的 CO2 值。
require(dplyr)
require(plotly)
require(akima)
require(plotly)
require(zoo)
require(raster)
plot_ly(df, x=~x,y=~y, z=~z, color =~calCO2)
重要的是你必须定义你将拥有的层。这些层必须由整个表面的值插值制成。所以:
为每一层定义您使用的数据。 插值 z 和 calCO2。这很重要,因为这是两件不同的事情。 z 插值将制作图形,而 calCO2 将制作颜色(浓度或其他)。在您来自 (https://plot.ly/r/3d-surface-plots/) 的图像中,颜色和 z 代表相同,而在这里,我猜您想代表 z 的表面并用 calCO2 对其进行着色。这就是为什么您需要为两者插入值。插值方法是一个世界,这里我只是做了一个简单的插值,我用平均值填充了 NA。代码如下:
## Define your layers in z range (by hand or use quantiles, percentiles, etc.)
df1 <- subset(df, z >= 0 & z <= 125) #layer between 0 and 150m
df2 <- subset(df, z > 125) #layer between 150 and max
#interpolate values for each layer and for z and co2
z1 <- interp(df1$x, df1$y, df1$z, extrap = TRUE, duplicate = "mean") #interp z layer 1 with spline interp
ifelse(anyNA(z1$z) == TRUE, z1$z[is.na(z1$z)] <- mean(z1$z, na.rm = TRUE), NA) #fill na cells with mean value
z2 <- interp(df2$x, df2$y, df2$z, extrap = TRUE, duplicate = "mean") #interp z layer 2 with spline interp
ifelse(anyNA(z2$z) == TRUE, z2$z[is.na(z2$z)] <- mean(z2$z, na.rm = TRUE), NA) #fill na cells with mean value
c1 <- interp(df1$x, df1$y, df1$calCO2, extrap = F, linear = F, duplicate = "mean") #interp co2 layer 1 with spline interp
ifelse(anyNA(c1$z) == TRUE, c1$z[is.na(c1$z)] <- mean(c1$z, na.rm = TRUE), NA) #fill na cells with mean value
c2 <- interp(df2$x, df2$y, df2$calCO2, extrap = F, linear = F, duplicate = "mean") #interp co2 layer 2 with spline interp
ifelse(anyNA(c2$z) == TRUE, c2$z[is.na(c2$z)] <- mean(c2$z, na.rm = TRUE), NA) #fill na cells with mean value
#THE PLOT
p <- plot_ly(showscale = TRUE) %>%
add_surface(x = z1$x, y = z1$y, z = z1$z, cmin = min(c1$z), cmax = max(c2$z), surfacecolor = c1$z) %>%
add_surface(x = z2$x, y = z2$y, z = z2$z, cmin = min(c1$z), cmax = max(c2$z), surfacecolor = c2$z) %>%
add_trace(data = df, x = ~x, y = ~y, z = ~z, mode = "markers", type = "scatter3d",
marker = list(size = 3.5, color = "red", symbol = 10))%>%
layout(title="Stack Exchange Plot")
p
【讨论】:
如果双三次样条对于给定情况是一种合理的插值/外插方法,这很好。您在回答中注意到“插值方法是一个世界,在这里我只是做了一个简单的插值”。这一点需要强调,一种具体的方法是否合理或有意义,完全取决于基础数据的性质。我们无法从 OP 的问题中看出这一点。 非常感谢您的回答!我非常喜欢这个。我看到的问题是它不能外推到该区域的全部范围(0-335,0-426),除非传感器位于极值点。查看带有 extrap=true 的 interp 命令,我认为这会解决这个问题。有没有办法解决这个问题?【参考方案2】:正如 Cesar 所指出的,您需要定义要在此 3d 系统中插值的“层”。
在这里,我提出了一种假设一层的方法(即 - 我使用沿 z 方向的所有点)。查看您的值表将帮助您定义中断发生的位置。您可以为您定义的每个“层”重复使用下面的代码。
> table(d$z)
0 50 120 130 155 178 226
7 10 1 3 8 1 1
既然你在处理空间数据,让我们使用 R 中的空间对象来解决这个问题。
首先,我将您的数据复制/粘贴到一个名为 d
的变量中。
# make d into a SpatialPointsDataFrame object
library(sp)
coords <- d[, c("x", "y")]
s <- SpatialPointsDataFrame(coords = coords, data = d)
# interpolate with a thin plate spline
# (or another interpolation method: kriging, inverse distance weighting).
library(raster)
library(fields)
tps <- Tps(coordinates(s), as.vector(d$calCO2))
p <- raster(s)
p <- interpolate(p, tps)
# plot raster, points, and contour lines
plot(p)
plot(s, add=T)
contour(p, add=T)
您可以想象根据点的z
值将数据拆分为层,然后重新运行此代码以生成每个层的插值。请务必阅读各种插值方法,以确定哪种方法最适合您的系统。一旦你有了这些层,就不需要再像上图那样将这些数据移植到 ploty 中了。
编辑:取基 --> ggplot --> 情节很简单:
# ggplot
library(ggplot2)
p <- ggplot(as.data.frame(p, xy = TRUE), aes(x, y, fill = layer)) +
geom_tile() +
geom_contour(aes(z = layer), color = "white") +
scale_fill_viridis_c() +
theme_minimal()
Here's some instructions on adding contour labels.
把它变成一个交互式的情节对象。
library(plotly)
ggplotly(p)
第一篇文章中的代码带你进入 3d。
【讨论】:
非常感谢您在这方面的帮助。我非常喜欢这种简单性。我注意到,这将所有 CO2 数据处理在同一平面上(相同的 z),但它的 CO2 随高度而变化。你怎么看?以上是关于如何在稀疏点之间插入数据以在 R & plotly 中绘制等高线图的主要内容,如果未能解决你的问题,请参考以下文章
如何映射分类变量以在 R Plotly 中为 3D 散点图中的点轮廓着色?
如何在 Greenplum/Postgres 中使用 PL/R 反序列化模型对象?