Oil & Gas Coding with R (Part 2)

In this tutorial, I will run you through the aRpsDCA package in R, actually developed by a guy many of you may know, Derrick Turk of Terminus Data Science. He’s likely not aware, but when I was first getting started I lobbed him a few questions and he gratefully answered and helped me to get better at this stuff. It’s a great package, but out of the box is a little tough to work through to apply specifically to the data I usually have access to, and that’s wonky state-level monthly production data.

Follow My Blog


About Me

I’m not a computer science guy. I was a reservoir engineer that went into finance and who then decided to start picking this up. Made me hugely faster and more efficient, but the key advantage is I’m used to working with all of the data sources you are, I know general oil and gas lingo, and I think like you. And domain knowledge trumps generalist coders almost every time. So my syntax probably sucks. There are probably faster ways to do things. But I will get you from point A to point B.


Quick Notes:

Everything is an object. It’s a key premise of both Python and R. Basically, anything that you have in memory is something that can be manipulated. Lists, strings, entire dataframes (think of these as an excel sheet).

The memory part is an important component too, as working with very large data files can get us near the edge of our computers’ memory, which can either signficiantly hamper performance or sometimes just refuse to run. So we need to make sure we’re deleting objects we don’t need anymore and doing general memory cleanup.

Also I would advise working through all of this example code line by line to see what happens, and help you to learn. Don’t just copy paste and run. Legitimately try to learn at each step.


Some real basic stuff


Base functionality

I’ll just run you through some really easy stuff in base R that’s ok to do without the use of packages.


Create a data frame (table)

df <- data.frame(As = c('a', 'A', 'aa', 'Aa'),
                 Bs = c('b', 'B', 'bb', 'Bb'),
                 number = 1,
                 numbers = seq(1, 4, 1))

What did I do here? I created a Data Frame, which is basically a table. The seq function basically says go from 1 to 4 in increments of 1. If you click on df under Data in the top right box (assuming you left the default configuration), it will pop up the table in another tab.

AsBsnumbernumbers
ab11
AB12
aabb13
AaBb14
Table 1: Example Data Frame

Simple math

Now, I can apply things to this table. For example.

df$new <- df$number + df$numbers

# OR the tidyverse version
library(tidyverse)
df <- df %>% mutate(new1 = number + numbers)

I can create a new column that is just adding the number columns to create a new column. I also show the tidyverse method.

AsBsnumbernumbersnewnew1
ab1122
AB1233
aabb1344
AaBb1455
Table 2: Example Data Frame With Adds

If I go and change the number column to 2, none of the other columns will be impacted (new and new1 will stay the same).


Taking the Average

All of the typical functions are available (+ – * /). You can also do averages.

df$averageNew <- mean(df$new)
AsBsnumbernumbersnewnew1averageNew
ab11223.5
AB12333.5
aabb13443.5
AaBb14553.5
Table 3: Example Data Frame With Average

Text manipulation

Another thing you may end up doing in R with our data is text manipulation/cleanup. Say I do a lot of analysis and I realize that the multiple characters were mistakes. Well, we can fix that pretty easily.

df$As <- substr(df$As, 1, 1)

substr in this case basically takes from position 1 to position 1 and that is your new column. Simple enough, but I use this one A LOT. Other good ones are toupper (capitalize all), and tolower (lower case all).

AsBsnumbernumbersnewnew1averageNew
ab11223.5
AB12333.5
abb13443.5
ABb14553.5
Table 4: Example Data Frame With Substr

Replacement with variables

I can also define a single value within a variable. For example:

x <- 2
df$number <- x
AsBsnumbernumbersnewnew1averageNew
ab21223.5
AB22333.5
abb23443.5
ABb24553.5
Table 5: Example Data Frame With Replacement

A little fancier

Since we have the tidyverse locked and loaded, let’s get a little more advanced. We will use a few functions, like select, summarise, mutate, and group_by.


Saving

I will just be doing calculations off of this base data frame. I always like to have a dataframe sitting around holding my basic data, and I create new objects off of that in the event that I screw everything up and have to retrace my steps. Of course if it gets too big, we can compress it and store it.

saveRDS(df, 'df.rds')
rm(df)

This will store your data so you don’t have to keep it in memory, and then I can delete the original. There’s always a balance between keeping things in memory and storing it, and there are actually some much faster ways nowadays with things like the fst package, though it can only really store raw data and nothing like shapefile data which I will get to at some point.


Reloading

To work with that data, we can either call it.

df <- readRDS('df.rds')

or we can just work off of that file so we don’t have to pull it into memory.


Quick dplyr example

For example, let’s try and do some quick table manipulation with dplyr.

df1 <- readRDS('df.rds') %>% subset(select = -c(Bs)) %>%
   mutate(new = number+numbers) %>%
   group_by(As) %>%
   summarise(averageNew = mean(new)) %>%
   ungroup()

So now we have a new dataframe. I had changed the number column to 2 earlier, so what I did was:

  1. Read in the dataframe
  2. Remove the Bs column
  3. Recalculate the new column to be number + numbers
  4. Grouped By the As column
  5. Averaged the new column to be averageNew
  6. Ungrouped

How does our new table look?

AsaverageNew
a3
A4
Table 6: Maniuplated Table

aRpsDCA

Now that the housekeeping is over (and hopefully you’re not bored to tears), let’s start doing some oil and gas specific stuff. For this, we will just be creating gas forecasts with single-segment or multi-segment decline.

First things first, install and load:

install.packages('aRpsDCA')
library(aRpsDCA)

Now, let’s create a quick monthly forecast. Why? Because I’m used to almost always having to work off of monthly data. You can do this off of daily as well, but I haven’t used it much for that. Note, that the default for the forecasts is Nominal Decline, but I almost always use effective as that’s what ARIES uses (I think? It’s been awhile), but luckily the package has a conversion for it. And also, when you try to recreate these in ARIES, use the H/ factor instead of B/ factor. Once again, I’m playing off memory here but I believe that’s correct.


Single-Segment

I will be using Modified Arps with a curtailment period at the beginning.

qiGas <- 20000 #mcfd
bGas <- 1 #B-factor
DiGas <- 0.9 #Effective Annual Decline (Initial)
DfGas <- 0.08 #Effective Annual Decline (Terminal)
curtail <- 5 #Months of curtailment
wellLife <- 30 #Years

gasFcst =  curtailed.q(arps.decline(
  qiGas*365.25, as.nominal(DiGas), bGas,
  as.nominal(DfGas)),
  curtail/12.0,seq(1/12, wellLife, by= (1/12)))/12

sum(gasFcst)

The first six variables are basic decline curve values. The gasFcst is just a vector of numbers that is your gas forecast.

  1. qiGas*365.25 -> Converts my daily rate to annual
  2. as.nominal(DiGas) -> Effective to Nominal Initial Decline
  3. as.nominal(DfGas) -> Effective to Nominal Terminal Decline
  4. curtail/12 -> Covert my flat monthly production to annual
  5. seq(1/12, wellLife, by = (1/12)) -> Perform the calculation in a sequence from month 1 to 360 in increments of 1.
  6. Divide the whole thing by 12 as it gives the values as an annual rate at that time.

You will see that this quick gas forecast is for 15,745,789 mcf, or 15.75 BCF.

Let’s check it out with highcharter.

library(highcharter)

highchart() %>% 
  hc_add_series(data.frame(months = seq(1, wellLife*12, 1), gas = gasFcst),
                type = 'spline',
                hcaes(x = months, y = as.integer(gas/30.4375)),
                color = 'red',
                name = 'Gas',
                marker = list(enabled = FALSE),
                showInLegend = FALSE) %>%
  hc_yAxis(type = 'logarithmic',
           title = list(text = '<b>Daily MCF</b>', 
                        style = list(fontSize = '18px')),
           labels = list(style = list(fontSize = '12px',
                                      fontWeight = 'bold')))%>%
  hc_xAxis(#type = 'logarithmic',
           title = list(text = '<b>Months</b>', 
                        style = list(fontSize = '18px')),
           labels = list(style = list(fontSize = '12px',
                                      fontWeight = 'bold'))) %>%
  hc_credits(enabled = TRUE, text = 'Powered by Highcharts', 
             href = "https://www.highcharts.com/") %>%
  hc_title(text = 'Modified Arps Type Curve', align = 'left') %>%
  hc_subtitle(text = 'aRpsDCA Package', align = 'left')

Figure 1: Example Type Curve – Single Segment


Multi-segment

In this example, we will pretend that the first segment is 36 months, and then new factors kick in at the end of that period (ie b-factor and Di changes).

period1 <- 36 #months Segment 1

qiGas1 <- 20000 #mcfd
bGas1 <- 1 #B-factor Segment 1
DiGas1 <- 0.9 #Effective Annual Decline (Initial) Segment 1
DfGas <- 0.08 #Effective Annual Decline (Terminal)
curtail <- 3 #Months of curtailment
wellLife <- 30 #Years

period2 <- (wellLife*12 - period1)/12 #Life of Segment 2

#Segment 1
gasFcst1 =  curtailed.q(arps.decline(
  qiGas*365.25, as.nominal(DiGas1), bGas1,
  as.nominal(0)),
  curtail/12.0,seq(1/12, period1/12, by= (1/12)))/12


qiGas2 <- gasFcst1[length(gasFcst1)]/30.4375 #Calculate New qiGas as last period of Segment 1
bGas2 <- 0.8 #New B Factor
DiGas2 <- 0.7 #New Initial Decline

#Segment 2
gasFcst2 =  curtailed.q(arps.decline(
  qiGas2*365.25, as.nominal(DiGas2), bGas2,
  as.nominal(DfGas)),
  0/12.0,seq(1/12, period2, by= (1/12)))/12

#Combine together
gasFcst <- append(gasFcst1, gasFcst2)

And in highcharter.

Figure 2: Example Type Curve – Multi-Segment


And a little Shiny to top it off

One of the best things about R is the ability to develop apps with Shiny, and it’s largely done so typical R folks (engineers or scientists) can actually display their results. To get started:

install.packages('shiny')

The rest of it is quite difficult, so I will just share the entire code block. Create a new R script, paste it in and save it as app.R.

library(shiny)
library(aRpsDCA)
library(highcharter)

ui <- fluidPage(
  
  # App title ----
  titlePanel("Type Curve Example"),
  
  # Sidebar layout with input and output definitions ----
  sidebarLayout(
    
    # Sidebar panel for inputs ----
    sidebarPanel(
      
      numericInput('qiGas', 'Gas IP, MCFD', value = 20000, min= 1, max = 100000),
      numericInput('bGas', 'B-Factor', value = 1),
      numericInput('DiGas', 'Effective Annual Di', value = 0.9),
      numericInput('DfGas', 'Effective Annual Terminal Decline', value = 0.08),
      numericInput('curtailGas', 'Flat Months', value = 4),
      numericInput('wellLife', 'Well Life, Yrs', value = 30)
      
    ),
    
    # Main panel for displaying outputs ----
    mainPanel(
      
      # Output: Decline Curve
      highchartOutput(outputId = "tcPlot"),
      textOutput('gasEUR')
      
    )
  )
)

server <- function(input, output, session) {
  
  output$tcPlot <- renderHighchart({
    
    gasFcst =  (curtailed.q(arps.decline(
      input$qiGas*365.25, as.nominal(input$DiGas), input$bGas,
      as.nominal(input$DfGas)),
      input$curtailGas/12.0,seq(1/12, input$wellLife, by= (1/12)))/12)
    
    output$gasEUR <- renderText(paste0('Gas EUR: ', as.integer(sum(gasFcst)/1000), ' MMCF'))
    

    
    highchart() %>% 
      hc_add_series(data.frame(months = seq(1, wellLife*12, 1), gas = gasFcst),
                    type = 'spline',
                    hcaes(x = months, y = as.integer(gas/30.4375)),
                    color = 'red',
                    name = 'Gas',
                    marker = list(enabled = FALSE),
                    showInLegend = FALSE) %>%
      hc_yAxis(type = 'logarithmic',
               title = list(text = '<b>Daily MCF</b>', 
                            style = list(fontSize = '18px')),
               labels = list(style = list(fontSize = '12px',
                                          fontWeight = 'bold')))%>%
      hc_xAxis(title = list(text = '<b>Months</b>', 
                     style = list(fontSize = '18px')),
        labels = list(style = list(fontSize = '12px',
                                   fontWeight = 'bold'))) %>%
      hc_credits(enabled = TRUE, text = 'Powered by Highcharts', 
                 href = "https://www.highcharts.com/") %>%
      hc_title(text = 'Modified Arps Type Curve', align = 'left') %>%
      hc_subtitle(text = 'aRpsDCA Package', align = 'left') 
    
  })
  
}

shinyApp(ui, server)

Once you save it, it will say Run App on the top right portion of the first R window. Run it and have fun!

Summary

Thanks for running through this quick analysis. In the next one, we will dig a bit deeper into functions and how to apply a quick economic model.

3 thoughts on “Oil & Gas Coding with R (Part 2)”

  1. Pingback: Oil & Gas Coding in Python (Part 2) - Shale Insights

  2. Pingback: Explaining the Shiny App (Oil & Gas Coding Series) - Shale Insights

  3. Pingback: EnergyNet Override Valuation (Oil & Gas Coding with R) - Shale Insights

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: