import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

!git clone https://www.github.com/yohanesnuwara/pyreservoir

import sys
sys.path.append('/content/pyreservoir/volumetrics')

from volumetrics import get_contours, compute_area, trapezoid, pyramidal, simpson

# load data
welldat = pd.read_csv('/content/pyreservoir/data/volumetric_data.csv')
x, y, z, h, poro, sw = welldat.x, welldat.y, welldat.depth, welldat.h, welldat.poro, welldat.sw


# plot contour map of OIP from the resulted N data"
from scipy.interpolate import griddata
import matplotlib.cm as cm

# define the z values for plotting
z = np.array([[z, h],
              [poro, sw]])

# grid x and y coordinate data
xi = np.linspace(min(x), max(x), 100)
yi = np.linspace(min(y), max(y), 100)
xi, yi = np.meshgrid(xi, yi)

# define subplots and colormap
figs, ax = plt.subplots(nrows=2, ncols=2, figsize=(12,10))

# define each subplot title
title = np.array([['Depth Map', 'Thickness Map'],
                  ['Porosity Map', 'Water Saturation Map']])

for i in range(2):
  for j in range(2):
    zi = griddata((x,y), z[i,j], (xi,yi), method='cubic')
    levels = np.linspace(min(z[i,j])-0.1, max(z[i,j])+0.1, 10)
    ax[i,j].contourf(xi, yi, zi, levels=levels, cmap='rainbow')  
    ax[i,j].plot(x, y, 'ko', ms=3) # plot the well points
    ax[i,j].set_title(title[i,j], pad=10, size=15)
    

# known
Bo = 1.435 # in RB/STB
owc = -4960 # oil water contact, in ft

N = welldat.h * welldat.poro * (1 - welldat.sw) / Bo
N = N * (1 / 5.61458) # convert the result from ft3/ft2 to standard STB/ft2
df = pd.DataFrame({"well_id": welldat.well_id, "N (STB/ft2)": N})
df = pd.merge(welldat, df, on='well_id')

# display result of OIP per ft2 calculation
df


# plot contour map of OIP from the resulted N data"
from scipy.interpolate import griddata

# define z value for plotting as: N (OIP per area)
z = N 

# grid N data with cubic interpolation method
zi = griddata((x,y),z,(xi,yi),method='cubic')

levels = np.arange(0, 1.1, 0.1)

# plot OIP per area
plt.figure(figsize=(10, 10))
figs = plt.contour(xi, yi, zi, levels=levels, colors=['0.25', '0.5', '0.25', '0.5', '0.25'], linewidths=[1.0, 0.5, 1, 0.5, 1])
plt.contourf(xi,yi,zi,levels=levels, cmap="rainbow")
plt.plot(x, y, 'ko', ms=3) # plot the well points
plt.clabel(figs, figs.levels[::2], inline=1, fontsize=10) # give labels for contours
plt.title("OIP contour map", pad=10, size=20)
plt.colorbar()
plt.show()


# get the individual contour lines and plot
contour_all = get_contours(figs, xi, yi, plot='Yes')

contour_area = compute_area(contour_all)

areas = pd.DataFrame({"OIP Contour(STB)": levels, "Area(ft2)": contour_area})
areas

oip_trapezo = trapezoid(contour_area, 0.1)
oip_pyramidal = pyramidal(contour_area, 0.1)
oip_simpson = simpson(contour_area, 0.1)


print("Oil in Place using Trapezoidal Method:", np.round((oip_trapezo / 1E+06), 3), "million STB")
print("Oil in Place using Pyramidal Method:", np.round((oip_pyramidal / 1E+06), 3), "million STB")
print("Oil in Place using Simpson 1/3 Method:", np.round((oip_simpson / 1E+06), 3), "million STB")