interpolate irregular x,y data points into regular grid for contour mapping

Question

I am a Geologist needing to create few hundred consistent contour maps in a project with varying x y z data sets.

Contouring irregular x y z data points involves creation of a 'grid' of interpolated (extrapolated) z values at a uniform x-y mesh. Outside R - this step is called 'gridding'. I am relatively new to R and trying to set a strong workflow to grid a large number of irregular data points. I am struggling!

On classic contour mapping software and workflow the steps are:

1. Read x y z data
2. Determine Area of Interest (AOI) for the final map. XMIN, XMAX, YMIN, YMAX
3. Determine Grid intervals (XINT, YINT) - sets the No of Rows and No of Columns to the 'grid' (NROW, NCOL)
4. Apply one of the needed interpolators - that create 'z' at the regular mesh/ grid (Common interpolators are: inverse distance, inverse square distance, weighted average, polynomial, kriging, spline, etc.)
5. Contour the resulting 'grid'

I am trying to script R to exactly follow the above sequence of steps for flexibility and control throughout the analysis.

df is the data frame consisting of an example data set.

     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

The Area of Interest is determined from study area as below:

    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

The values are: 754700, 791500,26196000,2658600, 36800, 39000, 500 respectively.

After a lot of failed trials - got the interp() function from package - akima to do the needed interpolation. Thanks to answer under "Plotting contours on an irregular grid"

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

This did not allow me to specify the AOI controls as desired. I tried using package MBA and yet working on creating the xy.est parameter (Grid mesh) as required input.

If proper 'grid' is generated, ggplot2 and other display functions are powerful and sufficient.

Are there proper 'Gridding' packages or 'Steps'. Thanks in advance.

Answer 1

I don't think you need to use other package than akima (and a graphic package like ggplot2). You can give AOI and number of 'grid' as interp's arguments xo and yo. And you can get xy.est parameter by interp2xyz(interp.obj).

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)