This work has been published by Mehdi Rezvandehy and Clayton V. Deutsch at Journal of Petroleum Geoscience https://doi.org/10.1144/petgeo2016-161. Python implementation for calculating uncertainty in spatial correlation (variogram) and how to incorporate this uncertainty in geostatistical modeling is presented.
Summary
Inferring a stable variogram model (spatial correlation) from widely-spaced data is a longstanding challenge. Correct uncertainty in the calculated variogram should be quantified and incorporated in geostatistical modeling. A new approach of variogram uncertainty is presented by computing the degrees of freedom. Variogram realizations are drawn from the uncertainty interval of lag distances honoring the correlation between lags. The realizations, respecting the correlation between lag distances of different directions, are transferred to geostatistical simulation to incorporate variogram uncertainty in the final numerical models. A realistic example is presented for variaogram uncertainty.
Python functions and data files needed to run this notebook are available via this link.
Required Libraries for this work
In [2]:
import pandas as pd
import numpy as np
import matplotlib
import pylab as plt
from matplotlib import gridspec
from matplotlib.offsetbox import AnchoredText
import geostatspy.GSLIB as GSLIB
import geostatspy.geostats as geostats
from IPython.display import clear_output
import os
os.environ['path']+=";.\\Codes" # All required codes are in this folder
import warnings
warnings.filterwarnings('ignore')
from IPython.display import HTML
A realistic example is presented for variogram uncertainty. The variable of interest is averaged permeability (in milliDarcies) over the main reservoir layer. Realizations of variogram for major and minor directions of continuity are calculated and used for geostatistical simulation one at the time: each variogram realization for each realization. This leads to incorporate uncertainty of the calculated variogram in final model.
In [3]:
# Looking at the data
Data=GSLIB_View('Data.out')
display(Data[0:10])
In [4]:
Data.describe()
Out[4]:
Figure below shows location map of porosity data with 62 samples (left) and histogram of data (right).
In [5]:
font = {'size' : 9}
matplotlib.rc('font', **font)
fig, axs = plt.subplots(figsize=(10, 5), dpi= 200, facecolor='w', edgecolor='k')
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1], wspace=0.12)
ax1=plt.subplot(gs[0])
plt.scatter(Data['X'], Data['Y'],c=Data['Per'],s=55,linewidths=1,edgecolors='k',
cmap='jet',vmin=0, vmax=14)
plt.colorbar(shrink=0.8,label='Permeability')
plt.xlabel('Easting (ft)',fontsize='12')
plt.ylabel('Northing (ft)',fontsize='12')
plt.title('Location Map of Permeability',fontsize=12)
ax1=plt.subplot(gs[1])
plt.xlabel('Permeability',fontsize='10')
plt.ylabel('Frequency',fontsize='10')
plt.hist(Data['Per'], bins=15,ec='black')
n=len(Data['Per'])
Mean=np.mean(Data['Per'])
SD=np.sqrt(np.var(Data['Per']))
Max=np.amax(Data['Per'])
Min=np.amin(Data['Per'])
#plt.xlim(0,20)
plt.title('Histogram of Permeability',fontsize=10)
plt.gca().set_aspect('0.7')
txt='n=%.0f \n $\mu$=%.3f \n $\sigma$=%.3f \n Max=%.1f \n Min=%.1f'
anchored_text = AnchoredText(txt %(n,Mean,SD,Max,Min), loc=1)
ax1.add_artist(anchored_text)
plt.show()
Variogram¶
In [6]:
#Experimental Variogram with Nromal Score and then fit a theoretical model
#Azimith 0: Experimental Variogram
tmin=-100;tmax=10000;xlag=1000;xltol=750;nlag=10;azm=0;atol=22.5;bandwh=99999;isill=1
Lag1, Var1, npp_az0=geostats.gamv(Data,'X', 'Y','NS:Per',tmin,tmax,xlag,xltol,nlag,azm,atol,bandwh,isill)
#Azimith 0: Fitted variogram Model Parameters
nlag=100; xlag=100
vario={'nst':2,'nug':0.0, 'cc1':0.2,'it1':2,'azi1':0,'hmaj1':1000,'hmin1':1000,
'cc2':0.8,'it2':1,'azi2':0,'hmaj2':15000,'hmin2':15000}
index_az0,Lag1_fit,Var1_fit,cov_az0,ro_az0=geostats.vmodel(nlag, xlag, azm, vario)
#Azimith 90: Experimental Variogram
tmin=-100;tmax=10000;xlag=1000;xltol=750;nlag=10;azm=90;atol=22.5;bandwh=99999;isill=1
Lag2, Var2, npp_az90=geostats.gamv(Data,'X', 'Y','NS:Por',tmin,tmax,xlag,xltol,nlag,azm,atol,bandwh,isill)
#Fitted Model: Fitted variogram Model Parameters
nlag=100; xlag=100
vario={'nst':2,'nug':0.0, 'cc1':0.2,'it1':2,'azi1':0,'hmaj1':1000,'hmin1':1000,
'cc2':0.8,'it2':1,'azi2':0,'hmaj2':4500,'hmin2':4500}
index_az90,Lag2_fit,Var2_fit,cov_az90,ro_az90=geostats.vmodel(nlag, xlag, azm, vario)
clear_output()
In [7]:
# Plot variogram models
font = {'size' : 12}
matplotlib.rc('font', **font)
fig, axs = plt.subplots(figsize=(12,12), dpi= 200, facecolor='w', edgecolor='k')
fig.subplots_adjust(hspace = .05, wspace=0.25)
#
#Label and Title of Variogram Plot
label1='Azimuth 0 ' ;symbol1='--or';label1_fit='Fitted Model ';symbol1_fit='-r'
label2='Azimuth 90 ';symbol2='--sb';label2_fit='Fitted Model ';symbol2_fit='-b'
xlabl='Distance (ft)'; xlim=[0,7000]; vlim=[0,1.38]; sill=1;loc=4
ax1=plt.subplot(1,2,1)
#Plot Variograms for Azimith 0
Lag1=np.delete(Lag1, np.where(Var1==0))
Var1=np.delete(Var1, np.where(Var1==0))
plt.plot(Lag1, Var1,symbol1,linewidth=2,markersize=7,label=label1)
plt.plot(Lag1_fit, Var1_fit,symbol1_fit,linewidth=2.5,label=label1_fit)
plt.legend(fontsize=13,loc=loc,frameon=False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.gca().set_aspect('4000')
plt.axhline(y=sill,linewidth=2,color='k',linestyle='--')
title='Variogram and Fitted Model_Azimuth 0'
plt.title(title,fontsize=15)
plt.xlabel(xlabl,fontsize=15)
ylabl='$\gamma (\mathbf{h}) $'
plt.ylabel(ylabl,rotation=0,fontsize=17,labelpad=15)
plt.xlim(float(xlim[0]),float(xlim[1]))
plt.ylim(float(vlim[0]),float(vlim[1]))
ax1=plt.subplot(1,2,2)
#Plot Variogramz for Azimith 90
Lag2=np.delete(Lag2, np.where(Var2==0))
Var2=np.delete(Var2, np.where(Var2==0))
plt.plot(Lag2, Var2,symbol2,linewidth=2,markersize=7,label=label2)
plt.plot(Lag2_fit, Var2_fit,symbol2_fit,linewidth=2.5,label=label2_fit)
plt.legend(fontsize=13,loc=loc,frameon=False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.gca().set_aspect('4000')
plt.axhline(y=sill,linewidth=2,color='k',linestyle='--')
title='Variogram and Fitted Model_Azimuth 90'
plt.title(title,fontsize=15)
plt.xlabel(xlabl,fontsize=15)
ylabl='$\gamma (\mathbf{h}) $'
plt.ylabel(ylabl,rotation=0,fontsize=17,labelpad=15)
plt.xlim(float(xlim[0]),float(xlim[1]))
plt.ylim(float(vlim[0]),float(vlim[1]))
clear_output()
plt.show()
Run VarUn Code for Calculating Variogram Uncertainty¶
Call VarUn Fortran code in python for calculating variogram uncertainty for Azimuth 0 and 90 degrees. 100 variogram realizations are simulated.
In [8]:
# Call VarUn Fortran code in Python for variogram declustering
txt=" Parameters for VarUn \n"\
+" *************************** \n"\
+"START OF PARAMETERS:\n"\
+"Data.out - file with data \n"\
+"2 3 0 - columns for X, Y, Z coordinates \n"\
+"10 - col number for variable \n"\
+"-100 1.0e21 - trimming limits \n"\
+"0 1 - global kriging variance yes(GK=1)/no(GK=0), simple(1)/ordinary(0) \n"\
+"50 100 200 - nx, xmn, xsiz (ignore if GK=0) \n"\
+"50 100 200 - ny, ymn, ysiz (ignore if GK=0) \n"\
+"1 0.5 1.0 - nz, zmn, zsiz (ignore if GK=0) \n"\
+"3 - grid increment, 1 considering all blocks (ignore if GK=0) \n"\
+"7 - number of lags \n"\
+"1000 - lag separation distance \n"\
+"750 - lag tolerance \n"\
+"2 - number of directions \n"\
+"0 22.5 99999 0.0 0.5 5.0 - azm,atol,bandh,dip,dtol,bandv \n"\
+"90 22.5 99999 0.0 0.5 5.0 - azm,atol,bandh,dip,dtol,bandv \n"\
+"1 - standardize experimental variogram \n"\
+"1 - standardize variogram realizations \n"\
+"1251 - random number seed \n"\
+"100 - number of variogram realizations \n"\
+"1.4 - range scale: enforce variogram range of first structure (<1 not apply) \n"\
+"0 -1.2 - variogram realizations limits for vriogram plot (from data if max<min) \n"\
+"./VariogramUncertainty/corr.out - correlation matrix between lag distance \n"\
+"./VariogramUncertainty/FOM.out - Variogram realizations of FOM approach \n"\
+"./VariogramUncertainty/FOM.ps - output file for variogam ploting of FOM approach \n"\
+"./VariogramUncertainty/DoF.out - Variogram realizations of DoF approach \n"\
+"./VariogramUncertainty/DoF.ps - output file for variogam ploting of DoF approach \n"\
+"./VariogramUncertainty/GK.out - Variogram realizations of global kriging approach \n"\
+"./VariogramUncertainty/GK.ps - output file for variogam ploting of global kriging approach \n"\
+"2 0.000 - nst, nugget effect \n"\
+"2 0.2 0.0 0.0 0.0 - it,cc,ang1,ang2,ang3 \n"\
+" 1000 1000 1 - a_hmax,a_hmin,a_vert \n"\
+"1 0.8 0.0 0.0 0.0 - it,cc,ang1,ang2,ang3 \n"\
+" 15000.000 4500.000 1 - a_hmax,a_hmin,a_vert \n"
f1 = open('temp', "w");
f1.write(txt)
f1.close()
! echo temp | VarUn
! rm temp
clear_output()
In [9]:
# Plot covariance matrix between lag distances
font = {'size' : 7}
matplotlib.rc('font', **font)
fig, ax=plt.subplots(figsize=(5, 5), dpi= 200, facecolor='w', edgecolor='k')
corr=GSLIB_View('./VariogramUncertainty/corr.out')
corr.drop('lag',axis=1, inplace=True)
corr=corr.values
colu_Az0=['Lag ' +str(i+1)+', Azimuth 0' for i in range(8)]
colu_Az90=['Lag ' +str(i+1)+', Azimuth 90' for i in range(8)]
columns=colu_Az0+colu_Az90
cax =ax.matshow(corr, cmap='jet', interpolation='nearest',vmin=-0.6, vmax=0.6)
fig.colorbar(cax,shrink=0.5,label='Correlation Coefficient')
plt.title('Correlation Matrix of Variogram Lags', fontsize=8,y=1.22)
for i in range(16):
for j in range(16):
c = corr[j,i]
ax.text(i, j, str(round(c,2)), va='center', ha='center',fontsize=5)
plt.axvline(x=7.5,linewidth =2,color='k')
plt.axhline(y=7.5,linewidth =2,color='k')
ax.set_xticks(np.arange(len(columns)))
ax.set_xticklabels(columns,fontsize=4.5, rotation='vertical')
ax.set_yticks(np.arange(len(columns)))
ax.set_yticklabels(columns,fontsize=4.5)
plt.show()
Correlation matrix between the sampled variogram realizations for eight lag distances of azimuth 0â—¦ and eight lag distances of azimuth 90â—¦.
Read Variogram Realization Calculated by VarUn¶
In [10]:
# Variogram uncertainty calculated by degree of freedom (DoF) approach works better than other approaches.
# Varigram uncertinaty calculated by DoF approach is read in python for vissualization.
DoF='./VariogramUncertainty/DoF.out'
Lag_f_0=[]; Gamma_f_0=[]
Lag_f_90=[]; Gamma_f_90=[]
nlag=100; xlag=100; azm=0
nsim=100
nst=np.zeros((nsim)); nug=np.zeros((nsim)); cc1=np.zeros((nsim)); it1=np.zeros((nsim)); azi1=np.zeros((nsim));
hmaj1_0=np.zeros((nsim)); hmin1_0=np.zeros((nsim)); hmaj1_90=np.zeros((nsim));hmin1_90=np.zeros((nsim))
#
cc2=np.zeros((nsim)); it2=np.zeros((nsim)); azi2=np.zeros((nsim));
hmaj2_0=np.zeros((nsim)); hmin2_0=np.zeros((nsim));hmaj2_90=np.zeros((nsim)); hmin2_90=np.zeros((nsim))
#
insim=0
with open(DoF) as f:
i=0
j=0
for line in f:
i=i+1
p = line.split()
if(i==2): nst[insim]=int(float(p[0])); nug[insim]=float(p[1])
if(i==3): cc1[insim]=float(p[1]); it1[insim]=int(float(p[0])); azi1[insim]=float(p[2])
if(i==4): hmaj1_0[insim]=float(p[0]); hmin1_0[insim]=float(p[0]); hmaj1_90[insim]=float(p[1]);hmin1_90[insim]=float(p[1])
if(i==5): cc2[insim]=float(p[1]); it2[insim]=int(float(p[0])); azi2[insim]=float(p[2])
if(i==6):
hmaj2_0[insim]=float(p[0]); hmin2_0[insim]=float(p[0]);hmaj2_90[insim]=float(p[1]); hmin2_90[insim]=float(p[1])
#Azimith 0: Variogram Realizations
vario_0={'nst':int(nst[insim]),'nug':nug[insim], 'cc1':cc1[insim],'it1':int(it1[insim]),
'azi1':azi1[insim],'hmaj1':hmaj1_0[insim],'hmin1':hmin1_0[insim],
'cc2':cc2[insim],'it2':int(it2[insim]),'azi2':azi2[insim],'hmaj2':hmaj2_0[insim],'hmin2':hmin2_0[insim]}
index_az0,Lag1_fit,Var1_fit,cov_az0,ro_az0=geostats.vmodel(nlag, xlag, azm, vario_0)
Lag_f_0.append(Lag1_fit); Gamma_f_0.append(Var1_fit)
#Azimith 90: Variogram Realizations
vario_90={'nst':int(nst[insim]),'nug':nug[insim], 'cc1':cc1[insim],'it1':int(it1[insim]),
'azi1':azi1[insim],'hmaj1':hmaj1_90[insim],'hmin1':hmin1_90[insim],
'cc2':cc2[insim],'it2':int(it2[insim]),'azi2':azi2[insim],'hmaj2':hmaj2_90[insim],'hmin2':hmin2_90[insim]}
index_az90,Lag2_fit,Var2_fit,cov_az90,ro_az90=geostats.vmodel(nlag, xlag, azm, vario_90)
Lag_f_90.append(Lag2_fit); Gamma_f_90.append(Var2_fit)
i=0
insim=insim+1
clear_output()
In [11]:
# Plot variogram realizations for Azimuth 0 and 90
fig=plt.figure(figsize=(6, 2), dpi= 200, facecolor='w', edgecolor='k')
font = {'size' : 5}
matplotlib.rc('font', **font)
ax1=plt.subplot(1,2,1)
xlim=[0,7000]; vlim=[0,1.2]; sill=1; fsize=5; loc=4; outfile='vargplt_Declustered.png'
plt.axhline(y=sill,linewidth=0.6,color='k',linestyle='--')
plt.title('Azimuth 0',fontsize=fsize*1.5)
plt.xlabel(xlabl,fontsize=fsize*1.5)
ylabl='$\gamma (\mathbf{h}) $'
plt.ylabel(ylabl,rotation=0,fontsize=fsize*1.9,labelpad=5)
plt.xlim(float(xlim[0]),float(xlim[1]))
plt.ylim(float(vlim[0]),float(vlim[1]))
#
symbol1= '-ob'; symbol2= '-g'; symbol3= '--'; label1='Experimental Variogram'
label2='Variogram Realization ';label3='Fitted Variogram (Declustered)'; xlabl='Distance (km)'
for i in range(nsim):
if(i==0):
plt.plot(Lag_f_0[i], Gamma_f_0[i],symbol2,color='gray',linewidth=0.5,markersize=1.5,label=label2)
else:
plt.plot(Lag_f_0[i], Gamma_f_0[i],symbol2,color='gray',linewidth=0.5,markersize=1.5)
plt.plot(Lag1, Var1,symbol1,linewidth=0.8,markersize=2,label=label1)
fsizet=fsize*1.2
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.legend(fontsize=str(fsizet),loc=loc,frameon=False)
ax2=plt.subplot(1,2,2)
xlim=[0,5000]; vlim=[0,1.3]; sill=1; fsize=5; loc=4; outfile='vargplt_Declustered.png'
plt.axhline(y=sill,linewidth=0.6,color='k',linestyle='--')
plt.title('Azimuth 90',fontsize=fsize*1.5)
plt.xlabel(xlabl,fontsize=fsize*1.5)
ylabl='$\gamma (\mathbf{h}) $'
plt.ylabel(ylabl,rotation=0,fontsize=fsize*1.9,labelpad=5)
plt.xlim(float(xlim[0]),float(xlim[1]))
plt.ylim(float(vlim[0]),float(vlim[1]))
#
symbol1= '-or'; symbol2= '-g'; symbol3= '--'; label1='Experimental Variogram'
label2='Variogram Realization ';label3='Fitted Variogram (Declustered)'; xlabl='Distance (km)'
for i in range(nsim):
if(i==0):
plt.plot(Lag_f_90[i], Gamma_f_90[i],symbol2,color='gray',linewidth=0.5,markersize=1.5,label=label2)
else:
plt.plot(Lag_f_90[i], Gamma_f_90[i],symbol2,color='gray',linewidth=0.5,markersize=1.5)
plt.plot(Lag2, Var2,symbol1,linewidth=0.8,markersize=2,label=label1)
fsizet=fsize*1.2
plt.legend(fontsize=str(fsizet),loc=loc,frameon=False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
fig.savefig(outfile, bbox_inches='tight')
100 realizations of variogram uncertainty for Azimuth 0 and 90. Realzations should the correlation between variables
Transfer Variogram Uncertainty in Geostatistical Modeling (SGS)¶
Variogram uncertainty is more influential in flow simulation because it affects the connectivity of rock properties.
In order to incorporate the uncertainty in variogram into our model, each calculated variogram realization is applied in conditional simulation one at the time; for exmaple to generate 100 realizations of geostatistical simulation, 100 variogram realizations are needed. In conventional geostatistical modeling, only one variogram model is applied that variogram uncertainty is not incorporated in the model.
Below shows 6 realization of sequential Gaussian simulation (SGS) with 6 diffent generated variogram realization.
In [12]:
nsim_=100
# Parameter file of SGS
for isim in range(nsim_):
txt=" Parameters for SGSIM \n"\
+" ********************** \n"\
+"START OF PARAMETERS:\n"\
+"Data.out -file with data\n"\
+"2 3 0 5 0 0 - columns for X,Y,Z,vr,wt,sec.var.\n"\
+"-998.0 1.0e21 - trimming limits\n"\
+"1 -transform the data (0=no, 1=yes)\n"\
+"./Simulation/sgsim.trn - file for output trans table\n"\
+"0 - consider ref. dist (0=no, 1=yes)\n"\
+"file.out - file with ref. dist distribution\n"\
+"1 0 - columns for vr and wt\n"\
+"0 27 - zmin,zmax (for tail extrapolation)\n"\
+"1 4.0 - lower tail option (1=linear), parameter\n"\
+"1 11.5 - upper tail option (1=linear), parameter\n"\
+"1 -debugging level: 0,1,2,3\n"\
+"./Simulation/sgsim.dbg -file for debugging output\n"\
+"./Simulation/sgsim_"+str(isim+1)+" -file for simulation output\n"\
+"1 -number of realizations to generate\n"\
+"100 50 100 -nx,xmn,xsiz\n"\
+"110 50 100 -ny,ymn,ysiz\n"\
+"1 1 1 -nz,zmn,zsiz\n"\
+str(10520+isim)+" -random number seed\n"\
+"4 24 -min and max original data for sim\n"\
+"24 -number of simulated nodes to use\n"\
+"1 -assign data to nodes (0=no, 1=yes)\n"\
+"1 3 -multiple grid search (0=no, 1=yes),num\n"\
+'0 -maximum data per octant (0=not used)\n'\
+str(hmaj2_0[isim])+" "+str(hmaj2_90[isim])+" 1.0 -maximum search radii\n"\
+"0 0.0 0.0 -angles for search ellipsoid\n"\
+"0 0.60 1.0 -ktype: 0=SK,1=OK,2=LVM,3=EXDR,4=COLC, corr and VRF\n"\
+"lvmfl.dat - file with LVM, EXDR, or COLC variable\n"\
+"1 - column for secondary variable\n"\
+str(nst[isim])+" "+str(nug[isim])+" -nst, nugget effect\n"\
+str(it1[isim])+" "+str(cc1[isim])+ " 0 0.0 0.0 -it,cc,ang1,ang2,ang3 \n"\
+" "+str(hmaj1_0[isim])+" "+str(hmin1_90[isim])+ " 1.0 -a_hmax, a_hmin, a_vert \n"\
+str(it2[isim])+" "+str(cc2[isim])+ " 0 0.0 0.0 -it,cc,ang1,ang2,ang3\n"\
+" "+str(hmaj2_0[isim])+" "+str(hmin2_90[isim])+ " 1.0 -a_hmax, a_hmin, a_vert \n"
#
f1 = open('temp', "w");
f1.write(txt)
f1.close()
! echo temp | sgsim
! rm temp
clear_output()
Variogram Realization 1 for SGS Simulation¶
In [13]:
isim=1
font = {'size' : 7}
matplotlib.rc('font', **font)
fig, axs = plt.subplots(figsize=(8, 5), dpi= 200, facecolor='w', edgecolor='k')
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2], wspace=0.25)
ax1=plt.subplot(gs[0])
vargunshow(isim)
ax2=plt.subplot(gs[1])
show_sim(isim)
Variogram Realization 2 for SGS Simulation¶
In [14]:
isim=2
font = {'size' : 7}
matplotlib.rc('font', **font)
fig, axs = plt.subplots(figsize=(8, 5), dpi= 200, facecolor='w', edgecolor='k')
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2], wspace=0.25)
ax1=plt.subplot(gs[0])
vargunshow(isim)
ax2=plt.subplot(gs[1])
show_sim(isim)
Variogram Realization 3 for SGS Simulation¶
In [15]:
isim=3
font = {'size' : 7}
matplotlib.rc('font', **font)
fig, axs = plt.subplots(figsize=(8, 5), dpi= 200, facecolor='w', edgecolor='k')
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2], wspace=0.25)
ax1=plt.subplot(gs[0])
vargunshow(isim)
ax2=plt.subplot(gs[1])
show_sim(isim)
Variogram Realization 4 for SGS Simulation¶
In [16]:
isim=4
font = {'size' : 7}
matplotlib.rc('font', **font)
fig, axs = plt.subplots(figsize=(8, 5), dpi= 200, facecolor='w', edgecolor='k')
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2], wspace=0.25)
ax1=plt.subplot(gs[0])
vargunshow(isim)
ax2=plt.subplot(gs[1])
show_sim(isim)
Variogram Realization 5 for SGS Simulation¶
In [17]:
isim=5
font = {'size' : 7}
matplotlib.rc('font', **font)
fig, axs = plt.subplots(figsize=(8, 5), dpi= 200, facecolor='w', edgecolor='k')
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2], wspace=0.25)
ax1=plt.subplot(gs[0])
vargunshow(isim)
ax2=plt.subplot(gs[1])
show_sim(isim)
Uncertainty in Sum of all Simulated Values with Variogram Uncertainty¶
After calculating 100 variogram realizations and incorporate them in geostatistical simulation, final distribution of uncertainty in final respond can be achieved. This could uncertainty in sum of all cell (uncertainty in resource estimation).
In [18]:
# Plot histogram of uncertainty for sum of values (resource estimation)
font = {'size' :10 }
#
matplotlib.rc('font', **font)
fig=plt.figure(figsize=(7, 4), dpi= 200, facecolor='w', edgecolor='k')
sum_=[]
for isim in range(nsim):
txt= './Simulation/sgsim_'+str(isim+1)
with open(txt) as f:
vt1=[]
for i in range(3):
next(f)
for line in f:
p = line.split()
vt1.append(float(p[0]))
sum_.append(sum(vt1))
histplt (sum_,25,title='Uncertainty in Sum with Variogram Uncertainty'
,xlabl='Sum',xlimt=None,ylimt=None,axt=None)
