CODEWITHSUNDEEP
rsurfaceauc
Cole Robertson
Cole Robertson

Reputation: 639

Method for calculating volume under a surface in R

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")))

Upvotes: 0

Views: 514

Answers (1)

the-mad-statter
the-mad-statter

Reputation: 8711

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

Upvotes: 1

Related Questions