R vectorize findInterval over large data sets

Question

I have two data frames that I'm using findInterval on. Wellbore data is data of the x,y, and z of a wellbore to produce oil (VSS = vertical subsea-depth, md = measured depth a.k.a. actual distance the drill bit traveled down the well). Perfs data is data where a wellbore has been perforated to allow for flow (top_perf = md, bot_perf = md).

Perfs:

Well_ID   top_perf    bot_perf    well_name   surface   ID  x   y   VSS
056-W        2808        2958        056-W     Ranger   2   0   0   0
056-W        3150        3250        056-W       Ranger 1   0   0   0
056-W        3150        3250        056-W       Ranger 2   0   0   0
056-W        3559        3664        056-W       UT 1   1   0   0   0
056-W        3559        3664        056-W       UT 2   2   0   0   0
057-W        2471        2952        057-W       Tar    1   0   0   0
057-W        2471        2952        057-W       Tar    2   0   0   0
058-W        2615        2896        058-W       Ranger 1   0   0   0
058-W        2615        2896        058-W       Ranger 2   0   0   0

Wellbore:

well_name   well_id      md      vss         x       y          
056-W        056-W       3260   -3251.46    4221436 4030454
056-W        056-W       3280   -3271.45    4221436 4030454
056-W        056-W       3300   -3291.45    4221435 4030453
056-W        056-W       3320   -3311.44    4221435 4030453
056-W        056-W       3340   -3331.44    4221434 4030453
056-W        056-W       3360   -3351.43    4221434 4030453
056-W        056-W       3380   -3371.43    4221433 4030453
056-W        056-W       3400   -3391.42    4221433 4030453

The goal is to find Perfs$top_perf and Perfs$bot_perf that are closest in value to Wellbore$md where Perfs$Well_ID = Wellbore$well_id and then extract the vss, x, and y from Wellbore and add it to Perfs. (I don't care about interpolating if it's in between, just need something that's close).

Here is my code to do this:

for(i in 1:dim(Perfs)[1]){
  if(Perfs$ID[i] == 1){
    Wellbore_temp <- Wellbore[which(Wellbore$well_id == Perfs[i,"Well_ID"]),]
    interval <- findInterval(Perfs[i,"top_perf"], Wellbore_temp$md)
    Perfs[i,c("x","y","VSS")] <- Wellbore_temp[interval, c("x","y","vss")]
  }else{
    Wellbore_temp <- Wellbore[which(Wellbore$well_id == Perfs[i,"Well_ID"]),]
    interval <- findInterval(Perfs[i,"bot_perf"], Wellbore_temp$md)
    Perfs[i,c("x","y","VSS")] <- Wellbore_temp[interval, c("x","y","vss")]
  }
}

This code does work, it's just far too slow for the application this will be used in. How can I get rid of the loop and do this in a more vectorized manner to speed things up? Also open to suggestions outside of findInterval.

Answer 1

Below I present a data.table solution. I've only tested it on the small subset of data that you've shown, and on that small a dataset it works out slower than your solution but I think it might scale better. If not, consider parallelizing.

If you've not used data.table before I think it's often pretty quick but syntax can be a bit convoluted. .SD refers to the subset of the wellbore data that joins to row i of the perfs data (iterating through .EACHI). This saves a mammoth join of everything to everything. Rather than using the findInterval function, I calculate an error (top_perf - md or bot_perf - md) and minimize absolute error. The advantage of this approach over a rolling join ('nearest') is that you can see what the error is, and filter if necessary.

library(data.table)

Perfs <- fread(input = 'Well_ID   top_perf    bot_perf    well_name   surface   ID  x   y   VSS
056-W        2808        2958        056-W     Ranger   2   0   0   0
056-W        3150        3250        056-W       Ranger 1   0   0   0
056-W        3150        3250        056-W       Ranger 2   0   0   0
056-W        3559        3664        056-W       UT_1   1   0   0   0
056-W        3559        3664        056-W       UT_2   2   0   0   0
057-W        2471        2952        057-W       Tar    1   0   0   0
057-W        2471        2952        057-W       Tar    2   0   0   0
058-W        2615        2896        058-W       Ranger 1   0   0   0
058-W        2615        2896        058-W       Ranger 2   0   0   0')

Wellbore <- fread(input = 'well_name   well_id      md      vss         x       y          
056-W        056-W       3260   -3251.46    4221436 4030454
056-W        056-W       3280   -3271.45    4221436 4030454
056-W        056-W       3300   -3291.45    4221435 4030453
056-W        056-W       3320   -3311.44    4221435 4030453
056-W        056-W       3340   -3331.44    4221434 4030453
056-W        056-W       3360   -3351.43    4221434 4030453
056-W        056-W       3380   -3371.43    4221433 4030453
056-W        056-W       3400   -3391.42    4221433 4030453')


#top
setkey(Wellbore, 'well_id')
setkey(Perfs, 'Well_ID', 'top_perf')
top_matched <- Wellbore[unique(Perfs), .SD[which.min(abs(top_perf-md)),.(md, top_perf, err=top_perf-md, x,y,vss)],nomatch=0, by=.EACHI]
setkey(top_matched, 'well_id', 'top_perf')
top_joined <- top_matched[Perfs]
top_joined[,`:=`(i.x=NULL, i.y=NULL,VSS=NULL)]
setnames(top_joined, old=c('err', 'x', 'y', 'vss'), new=paste0('top_', c('err', 'x', 'y', 'vss')))

#bottom
setkey(Perfs, 'Well_ID', 'bot_perf')
bot_matched <- Wellbore[unique(Perfs), .SD[which.min(abs(bot_perf-md)),.(md, bot_perf, err=bot_perf-md, x,y,vss)],nomatch=0, by=.EACHI]
setkey(bot_matched, 'well_id', 'bot_perf')
bot_joined <- bot_matched[Perfs]
bot_joined[,`:=`(i.x=NULL, i.y=NULL,VSS=NULL)]
setnames(bot_joined, old=c('err', 'x', 'y', 'vss'), new=paste0('bot_', c('err', 'x', 'y', 'vss')))


answer <- cbind(top_joined[,c(1:2,9:11,3:7), with=F], bot_joined[,3:7,with=F])

# well_id   md well_name surface ID top_perf top_err   top_x   top_y  top_vss bot_perf bot_err
# 1:   056-W 3260     056-W  Ranger  2     2808    -452 4221436 4030454 -3251.46     2958    -302
# 2:   056-W 3260     056-W  Ranger  1     3150    -110 4221436 4030454 -3251.46     3250     -10
# 3:   056-W 3260     056-W  Ranger  2     3150    -110 4221436 4030454 -3251.46     3250     -10
# 4:   056-W 3400     056-W    UT_1  1     3559     159 4221433 4030453 -3391.42     3664     264
# 5:   056-W 3400     056-W    UT_2  2     3559     159 4221433 4030453 -3391.42     3664     264
# 6:   057-W   NA     057-W     Tar  1     2471      NA      NA      NA       NA     2952      NA
# 7:   057-W   NA     057-W     Tar  2     2471      NA      NA      NA       NA     2952      NA
# 8:   058-W   NA     058-W  Ranger  1     2615      NA      NA      NA       NA     2896      NA
# 9:   058-W   NA     058-W  Ranger  2     2615      NA      NA      NA       NA     2896      NA
# bot_x   bot_y  bot_vss
# 1: 4221436 4030454 -3251.46
# 2: 4221436 4030454 -3251.46
# 3: 4221436 4030454 -3251.46
# 4: 4221433 4030453 -3391.42
# 5: 4221433 4030453 -3391.42
# 6:      NA      NA       NA
# 7:      NA      NA       NA
# 8:      NA      NA       NA
# 9:      NA      NA       NA

Answer 2

Found the answer to the question here: Join R data.tables where key values are not exactly equal--combine rows with closest times

Based on the thoughts of a data.table provided by @ds440

Here is the code I used and it runs very fast:

Perf.Data <- Perfs


Wellbore.Perfs <- data.table(Wellbore[,c("well_id","md","vss")])
Spotfire.Top.Perf <- data.table(Perf.Data[,c("Well_ID","top_perf", "bot_perf")])
Spotfire.Bot.Perf <- data.table(Perf.Data[,c("Well_ID","bot_perf", "top_perf")])

#Change the column names to match up with Wellbore.Perfs
#Add in the bot_perf to .top.perf and the top_perf to the .bot.perf is done to make these unique and ensure everything is captured from the perfs table
colnames(Spotfire.Top.Perf) <- c("well_id","md", "bot_perf")
colnames(Spotfire.Bot.Perf) <- c("well_id","md","top_perf")

#set key to join on
setkey(Wellbore.Perfs, "well_id","md")

#roll = "nearest" will take the nearest value of md in .top.perf or .bot.perf and match it to the md in wellbore.perfs where Well_ID = Well_ID
Perfs.Wellbore.Top <- Wellbore.Perfs[Spotfire.Top.Perf, roll = "nearest"]
Perfs.Wellbore.Bot <- Wellbore.Perfs[Spotfire.Bot.Perf, roll = "nearest"]