【问题标题】:plotting and coloring data on irregular grid在不规则网格上绘制和着色数据
【发布时间】:2011-08-02 06:39:11
【问题描述】:

我有 (x, y, z) 形式的数据,其中 x 和 y 不在常规网格上。我希望显示这些数据的 2D 颜色图,其中强度(例如灰度)映射到 z 变量。一个明显的解决方案是在规则网格上进行插值(见下文),

d <- data.frame(x=runif(1e3, 0, 30), y=runif(1e3, 0, 30))
d$z = (d$x - 15)^2 + (d$y - 15)^2


library(akima)
d2 <- with(d, interp(x, y, z, xo=seq(0, 30, length = 30),
                     yo=seq(0, 30, length = 50), duplicate="mean"))

pal1 <- grey(seq(0,1,leng=500))
with(d2, image(sort(x), sort(y), z, useRaster=TRUE, col = pal1))
points(d$x, d$y, col="white", bg=grey(d$z/max(d$z)), pch=21, cex=1,lwd=0.1)

但是,这会丢失初始网格的信息(具有实际数据的点的位置),这些信息在某些位置可能非常精细或非常粗糙。我更喜欢使用三角形的 delaunay 平铺,它准确地代表了原始数据点的实际位置和密度。

理想的解决方案是

  • 在绘图函数之外计算曲面细分,以便可以使用ggplot2lattice 或基本图形绘制生成的多边形

  • 要快。在我的实际示例中(约 1e5 分),通过 deldir 计算镶嵌可能非常缓慢。

“镶嵌”是指 Delaunay 三角形或 Voronoi 图,尽管我更喜欢前者。然而,它带来了额外的复杂性,即根据原始数据点对每个三角形的颜色进行插值。

【问题讨论】:

    标签: r delaunay voronoi spatstat tesselation


    【解决方案1】:

    这是基于maptools 包中的dirichlet 的解决方案,

    d <- data.frame(x=runif(1e3, 0, 30), y=runif(1e3, 0, 30))
    d$z = (d$x - 15)^2 + (d$y - 15)^2
    
    library(spatstat) 
    library(maptools)
    
    W <- ripras(df, shape="rectangle") 
    W <- owin(c(0, 30), c(0, 30)) 
    X <- as.ppp(d, W=W) 
    Y <- dirichlet(X) 
    Z <- as(Y, "SpatialPolygons") 
    plot(Z, col=grey(d$z/max(d$z)))
    

    我仍然不确定从这个 SpatialPolygons 类中提取多边形的方法。

    另外,如果有一种简单的方法可以为相关的 delaunay 镶嵌生成“正确”颜色,我想听听。

    【讨论】:

    • “正确的”颜色是什么意思?您是否仍在尝试将正确的灰色阴影附加到正确的瓷砖上?如果是这样,Z 的@plotOrder 插槽有帮助吗?我的意思是plot(Z, col=grey(d$z/max(d$z))[Z@plotOrder])
    • 确实,我在意识到真正的问题是 delaunay 瓦片比 z 值更多之前尝试了这个。试试上面的例子 delaunay(), Y
    【解决方案2】:

    这是一个使用deldir的格子解决方案

    d <- data.frame(x=runif(1e3, 0, 30), y=runif(1e3, 0, 30))
    d$z = (d$x - 15)^2 + (d$y - 15)^2
    
    pal1 <- grey(seq(0,1,leng=500))
    library(latticeExtra)
    
     levelplot(z~x*y, data=d,
               panel = function(...) panel.voronoi(..., points=FALSE),
               interpolate=TRUE,
               col.regions = colorRampPalette(pal1)(1e3), cut=1e3)
    

    【讨论】:

      猜你喜欢
      • 2018-08-03
      • 2019-09-24
      • 1970-01-01
      • 1970-01-01
      • 2018-08-09
      • 1970-01-01
      • 2014-09-06
      • 2019-04-17
      相关资源
      最近更新 更多