Method for calculating volume under a surface in R

Question

I'm trying to calculate the volume under a 3d surface in R.

My data, dat, looks like:

          0.003     0.019     0.083      0.25       0.5         1
0     1.0000000 0.8884265 0.8603268 0.7719994 0.7443621 0.6571405
0.111 0.6909722 0.6775000 0.6443750 0.6243750 0.5914730 0.5698242
0.25  0.5847205 0.6022367 0.5572917 0.5432991 0.5170673 0.4835819
0.429 0.5210938 0.5139063 0.4995312 0.4864062 0.4648636 0.4163698
0.667 0.4363103 0.4526562 0.4321859 0.4027519 0.4046011 0.3661616
1     0.3958333 0.4167468 0.3964428 0.3810459 0.3486328 0.3487930

where x = rownames(dat), y = colnames(dat) and z = dat.

I've looked here, here, and here, but can't seem to figure out how to apply those answers to my use case.

Here's a reproducible version of my data:

dat = structure(c(1,0.690972222222222,0.584720477386935,0.52109375,0.436310279187817,0.395833333333333,0.888426507537688,0.6775,0.602236675126904,0.51390625,0.45265625,0.416746794871795,0.860326776649746, 0.644375, 0.557291666666667,0.49953125,0.432185913705584,0.396442819148936,0.771999378109453,0.624375,0.543299129353234,0.48640625,0.402751865671642,0.381045854271357,0.744362113402062,0.591472989949749,0.517067307692308,0.464863578680203,0.404601130653266,0.3486328125,0.657140544041451,0.56982421875,0.483581852791878,0.41636981865285,0.366161616161616,0.348792989417989),.Dim = c(6L, 6L), .Dimnames = list(c("0","0.111","0.25","0.429","0.667","1"),c("0.003","0.019","0.083","0.25","0.5","1")))

Answer 1

You could use the getVolume() function provided in the answer you linked here provided that your matrix were in the requisite dataframe format.

Here is some code to make that dataframe:

df <- expand.grid(x = as.numeric(rownames(dat)), y = as.numeric(colnames(dat)))
df$z = as.vector(dat)

Then define the function and apply:

library(geometry)

getVolume=function(df) {
  #find triangular tesselation of (x,y) grid
  res=delaunayn(as.matrix(df[,-3]),full=TRUE,options="Qz")
  #calulates sum of truncated prism volumes
  sum(mapply(function(triPoints,A) A/3*sum(df[triPoints,"z"]),
             split.data.frame(res$tri,seq_along(res$areas)),
             res$areas))
}

getVolume(df)

[1] 0.4714882