将不规则的x、y数据点插值为规则网格,以便进行等高线绘制。

5

我是一名地质学家,需要在一个包含不同x、y、z数据集的项目中创建几百个一致的等高线图。

对于不规则的x、y、z数据点,需要创建插值(外推)后的z值网格,并在统一的x-y网格上进行。在R之外,这一步被称为“网格化”。我相对来说是R的新手,正在努力建立强大的工作流程,以网格化大量的不规则数据点。但我遇到了困难!

经典等高线制图软件和工作流程的步骤如下:

  1. 读取x、y、z数据
  2. 确定最终地图的兴趣区域(AOI)。XMIN、XMAX、YMIN、YMAX
  3. 确定网格间隔(XINT、YINT),设置'grid'的行数和列数(NROW、NCOL)
  4. 应用所需的插值器之一,创建正常网格/网格(常见的插值器有:反距离、反平方距离、加权平均、多项式、克里金、样条等)
  5. 等高线化所得的“grid”

我正在尝试使用R脚本完全按照以上步骤进行操作,以便在分析过程中获得灵活性和控制。

df是由示例数据集组成的数据框。

     wellid property           z        x       y
    060010        1 0.008849558 756994.5 2637732
    009410        1 0.260162602 760190.9 2622262
    009910        1 0.115044248 760898.7 2637466
    051110        1 0.109243697 761690.2 2630985
    065610        1 0.066666667 763064.1 2620929
    011010        1 0.000000000 763089.3 2630888
    035210        1 0.022556391 765942.4 2625944
    052510        1 0.157894737 767058.1 2650034
    006610        1 0.045045045 768265.0 2645318
    009010        1 0.378151261 768471.8 2636731
    011210        1 0.028776978 771393.8 2629001
    064810        1 0.428571429 771394.1 2650776
    009110        1 0.064220183 775332.6 2648531
    011410        1 0.148760331 778324.8 2633905
    065010        1 0.514851485 780480.9 2654874
    052410        1 0.173913043 780961.0 2637571
    064110        1 0.019417476 781001.5 2650994
    009310        1 0.037383178 783904.7 2641130
    010810        1 0.041237113 786200.6 2652417
    052610        1 0.150537634 788007.5 2654005

感兴趣的区域是根据以下研究区域确定的:

    xmin <- signif(min(wellcoords$x),4) - 1000
    xmax <- signif(max(wellcoords$x),4) +1000
    ymin <- signif(min(wellcoords$y),5) - 1000
    ymax <- signif(max(wellcoords$y),5) +1000
    xrange <- xmax-xmin
    yrange <- ymax-ymin
    gridint <- 500     # grid interval is set same for xint and yint

数值分别为:754700、791500、26196000、2658600、36800、39000、500。

经历了许多失败的尝试后,我们从akima包中得到了必要的插值函数interp()。感谢在“在不规则网格上绘制轮廓”的回答下提供帮助。

    fld<- with(df, interp(x=df$x, y=df$y, z=df$z, xo=xcoord, yo=ycoord, linear = FALSE, extrap = TRUE))

这使我无法按照预期指定AOI控件。我尝试使用MBA包,但是需要创建所需输入的xy.est参数(网格网格)。

如果生成适当的“网格”,则ggplot2和其他显示功能非常强大且充足。

是否有适当的“分网”包或步骤。

谢谢!
1个回答

8

我认为除了使用akima(和像ggplot2这样的图形包)之外,您不需要使用其他软件包。您可以将AOI和“网格”数量作为interp的参数xoyo。您可以通过interp2xyz(interp.obj)获得xy.est参数。

df <- "your example data set"
# I didn't know What wellcoords were, so I treated df as wellcoords. These values are different from what you said.
xmin <- signif(min(df$x),4) - 1000  # 756000
xmax <- signif(max(df$x),4) + 1000  # 789000
ymin <- signif(min(df$y),5) - 1000  # 2619900
ymax <- signif(max(df$y),5) + 1000  # 2655900
gridint <- 500

library(akima)
fld<- with(df, interp(x = x, y = y, z = z, linear = FALSE, extrap = TRUE,
                      xo=seq(xmin, xmax, length=gridint), 
                      yo=seq(ymin, ymax, length=gridint)))  # give AOI and number of 'grid'
# check whether the conditions are met.
length(fld$x); length(fld$y); length(fld$z); range(fld$x); range(fld$y)
  # 500, 500, 250000 (=500^2), 756000 789000, 2619900 2655900,   # all OK

contour(fld)   # Left graph (most basic graphic output)

fld2 <- as.data.frame(interp2xyz(fld))  # the xy.est parameter (data.frame)
library(ggplot2)
ggplot(fld2, aes(x=x, y=y, z=z)) + geom_contour()  # Right graph (simple example)

enter image description here


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接