## GeostatsPy: Variogram Calculation for Subsurface Data Analytics in Python 


### Michael Pyrcz, Associate Professor, University of Texas at Austin 

##### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)


### PGE 383 Exercise: Basic Univariate Summary Statistics and Data Distribution Representativity Plotting in Python with GeostatsPy

Here's a simple workflow with some basic variogram calcuation with irregularly sampled data. This should help you get started with variogram caculation in subsurface modeling.

#### Spatial Continuity 

**Spatial Continuity** is the correlation between values over distance.

* No spatial continuity – no correlation between values over distance, random values at each location in space regardless of separation distance.

* Homogenous phenomenon have perfect spatial continuity, since all values as the same (or very similar) they are correlated. 

We need a statistic to quantify spatial continuity! A convenient method is the Semivariogram.

#### The Semivariogram

Function of difference over distance.

* The expected (average) squared difference between values separated by a lag distance vector (distance and direction), $h$:

\begin{equation}
\gamma(\bf{h}) = \frac{1}{2 N(\bf{h})} \sum^{N(\bf{h})}_{\alpha=1} (z(\bf{u}_\alpha) - z(\bf{u}_\alpha + \bf{h}))^2  
\end{equation}

where $z(\bf{u}_\alpha)$ and $z(\bf{u}_\alpha + \bf{h})$ are the spatial sample values at tail and head locations of the lag vector respectively.

* Calculated over a suite of lag distances to obtain a continuous function.

* the $\frac{1}{2}$ term converts a variogram into a semivariogram, but in practice the term variogram is used instead of semivariogram.
* We prefer the semivariogram because it relates directly to the covariance function, $C_x(\bf{h})$ and univariate variance, $\sigma^2_x$:

\begin{equation}
C_x(\bf{h}) = \sigma^2_x - \gamma(\bf{h})
\end{equation}

Note the correlogram is related to the covariance function as:

\begin{equation}
\rho_x(\bf{h}) = \frac{C_x(\bf{h})}{\sigma^2_x}
\end{equation}

The correlogram provides of function of the $\bf{h}-\bf{h}$ scatter plot correlation vs. lag offset $\bf{h}$.  

\begin{equation}
-1.0 \le \rho_x(\bf{h}) \le 1.0
\end{equation}

#### Variogram Observations

The following are common observations for variograms that should assist with their practical use.

##### Observation \#1 - As distance increases, variability increase (in general).

This is common since in general, over greater distance offsets, there is often more difference between the head and tail samples.

In some cases, such as with spatial cyclicity of the hole effect variogram model the variogram may have negative slope over somelag distance intervals

Negative slopes at lag distances greater than half the data extent are often caused by too few pairs for a reliable variogram calculation

##### Observation \#2 - Calculated with over all possible pairs separated by lag vector, $\bf{𝐡}$.

We scan through the entire data set, searching for all possible pair combinations with all other data.  We then calculate the variogram as one half the expectation of squared difference between all pairs.

More pairs results in a more reliable measure.

##### Observation \#3 - Need to plot the sill to know the degree of correlation.

**Sill** is the variance, $\sigma^2_x$

Given stationarity of the variance, $\sigma^2_x$, and variogram $\gamma(\bf{h})$:

we can define the covariance function:

\begin{equation}
C_x(\bf{h}) = \sigma^2_x - \gamma(\bf{h})
\end{equation}

The covariance measure is a measure of similarity over distance (the mirror image of the variogram as shown by the equation above).

Given a standardized distribution $\sigma^2_x = 1.0$, the covariance, $C_x(\bf{h})$, is equal to the correlogram, $\rho_x(\bf{h})$: 

\begin{equation}
\rho_x(\bf{h}) = \sigma^2_x - \gamma(\bf{h})
\end{equation}

##### Observation \#4 - The lag distance at which the variogram reaches the sill is know as the range.

At the range, knowing the data value at the tail location provides no information about a value at the head location of the lag distance vector.

##### Observation \#5 - The nugget effect, a discontinuity at the origin

Sometimes there is a discontinuity in the variogram at distances less than the minimum data spacing.  This is known as **nugget effect**.

The ratio of nugget / sill, is known as relative nugget effect (%). Modeled as a discontinuity with no correlation structure that at lags, $h \gt \epsilon$, an infinitesimal lag distance, and perfect correlation at $\bf{h} = 0$.
Caution when including nuggect effect in the variogram model as measurement error, mixing populations cause apparent nugget effect

This exercise demonstrates the semivariogram calculation with GeostatsPy. The steps include:

1. generate a 2D model with sequential Gaussian simulation
2. sample from the simulation
3. calculate and visualize experimental semivariograms

#### Objective 

In the PGE 383: Stochastic Subsurface Modeling class I want to provide hands-on experience with building subsurface modeling workflows. Python provides an excellent vehicle to accomplish this. I have coded a package called GeostatsPy with GSLIB: Geostatistical Library (Deutsch and Journel, 1998) functionality that provides basic building blocks for building subsurface modeling workflows. 

The objective is to remove the hurdles of subsurface modeling workflow construction by providing building blocks and sufficient examples. This is not a coding class per se, but we need the ability to 'script' workflows working with numerical methods.    

#### Getting Started

Here's the steps to get setup in Python with the GeostatsPy package:

1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). 
2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. 
3. In the terminal type: pip install geostatspy. 
4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. 

You will need to copy the data file to your working directory.  They are available here:

* Tabular data - sample_data.csv at https://git.io/fh4gm.

There are exampled below with these functions. You can go here to see a list of the available functions, https://git.io/fh4eX, other example workflows and source code. 

#### Load the required libraries

The following code loads the required libraries.


```python
import geostatspy.GSLIB as GSLIB                       # GSLIB utilies, visualization and wrapper
import geostatspy.geostats as geostats                 # GSLIB methods convert to Python    
```

We will also need some standard packages. These should have been installed with Anaconda 3.


```python
import os                                               # to set current working directory 
import numpy as np                                      # arrays and matrix math
import pandas as pd                                     # DataFrames
import matplotlib.pyplot as plt                         # plotting
```

If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  


#### Set the working directory

I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time).  Also, in this case make sure to place the required (see above) GSLIB executables in this directory or a location identified in the environmental variable *Path*.


```python
os.chdir("c:/PGE383")                                   # set the working directory
```

#### Loading Tabular Data

Here's the command to load our comma delimited data file in to a Pandas' DataFrame object. 


```python
df = pd.read_csv("sample_data.csv")                     # read a .csv file in as a DataFrame
#print(df.iloc[0:5,:])                                  # display first 4 samples in the table as a preview
df.head()                                               # we could also use this command for a table preview 
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>X</th>
      <th>Y</th>
      <th>Facies</th>
      <th>Porosity</th>
      <th>Perm</th>
      <th>AI</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>100.0</td>
      <td>900.0</td>
      <td>1.0</td>
      <td>0.100187</td>
      <td>1.363890</td>
      <td>5110.699751</td>
    </tr>
    <tr>
      <th>1</th>
      <td>100.0</td>
      <td>800.0</td>
      <td>0.0</td>
      <td>0.107947</td>
      <td>12.576845</td>
      <td>4671.458560</td>
    </tr>
    <tr>
      <th>2</th>
      <td>100.0</td>
      <td>700.0</td>
      <td>0.0</td>
      <td>0.085357</td>
      <td>5.984520</td>
      <td>6127.548006</td>
    </tr>
    <tr>
      <th>3</th>
      <td>100.0</td>
      <td>600.0</td>
      <td>0.0</td>
      <td>0.108460</td>
      <td>2.446678</td>
      <td>5201.637996</td>
    </tr>
    <tr>
      <th>4</th>
      <td>100.0</td>
      <td>500.0</td>
      <td>0.0</td>
      <td>0.102468</td>
      <td>1.952264</td>
      <td>3835.270322</td>
    </tr>
  </tbody>
</table>
</div>



We will work by-facies, that is separating sand and shale facies and working with them separately.  This command extracts the sand and shale 'Facies" into new DataFrames for our analysis.


```python
# Separate Sample Data by-facies - Deepcopy
df_sand = pd.DataFrame.copy(df[df['Facies'] == 1]).reset_index() # copy only 'Facies' = sand records
df_shale = pd.DataFrame.copy(df[df['Facies'] == 0]).reset_index() # copy only 'Facies' = shale records
df_sand.head()                                          # preview the sand only DataFrame
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>index</th>
      <th>X</th>
      <th>Y</th>
      <th>Facies</th>
      <th>Porosity</th>
      <th>Perm</th>
      <th>AI</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>100.0</td>
      <td>900.0</td>
      <td>1.0</td>
      <td>0.100187</td>
      <td>1.363890</td>
      <td>5110.699751</td>
    </tr>
    <tr>
      <th>1</th>
      <td>8</td>
      <td>100.0</td>
      <td>100.0</td>
      <td>1.0</td>
      <td>0.137453</td>
      <td>5.727603</td>
      <td>5823.241783</td>
    </tr>
    <tr>
      <th>2</th>
      <td>9</td>
      <td>200.0</td>
      <td>900.0</td>
      <td>1.0</td>
      <td>0.137062</td>
      <td>14.771314</td>
      <td>5621.146994</td>
    </tr>
    <tr>
      <th>3</th>
      <td>10</td>
      <td>200.0</td>
      <td>800.0</td>
      <td>1.0</td>
      <td>0.125984</td>
      <td>10.675436</td>
      <td>4292.700500</td>
    </tr>
    <tr>
      <th>4</th>
      <td>18</td>
      <td>300.0</td>
      <td>900.0</td>
      <td>1.0</td>
      <td>0.111516</td>
      <td>27.999817</td>
      <td>4183.466773</td>
    </tr>
  </tbody>
</table>
</div>



#### Summary Statistics for Tabular Data

The table includes X and Y coordinates (meters), Facies 1 and 0 (1 is sandstone and 0 interbedded sand and mudstone), Porosity (fraction), and permeability as Perm (mDarcy). 

There are a lot of efficient methods to calculate summary statistics from tabular data in DataFrames. The 'describe' command provides count, mean, minimum, maximum, and quartiles all in a nice data table. We use transpose just to flip the table so that features are on the rows and the statistics are on the columns.

Let's look at and compare the sumamry statistics for sand and shale.


```python
df_sand.describe().transpose()                          # summary table of sand only DataFrame statistics
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>count</th>
      <th>mean</th>
      <th>std</th>
      <th>min</th>
      <th>25%</th>
      <th>50%</th>
      <th>75%</th>
      <th>max</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>index</th>
      <td>162.0</td>
      <td>117.981481</td>
      <td>59.058816</td>
      <td>0.000000</td>
      <td>75.250000</td>
      <td>117.500000</td>
      <td>157.750000</td>
      <td>260.000000</td>
    </tr>
    <tr>
      <th>X</th>
      <td>162.0</td>
      <td>831.111111</td>
      <td>238.857269</td>
      <td>50.000000</td>
      <td>800.000000</td>
      <td>945.000000</td>
      <td>975.000000</td>
      <td>1005.000000</td>
    </tr>
    <tr>
      <th>Y</th>
      <td>162.0</td>
      <td>526.987654</td>
      <td>142.707797</td>
      <td>100.000000</td>
      <td>469.000000</td>
      <td>509.000000</td>
      <td>549.000000</td>
      <td>939.000000</td>
    </tr>
    <tr>
      <th>Facies</th>
      <td>162.0</td>
      <td>1.000000</td>
      <td>0.000000</td>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>1.000000</td>
    </tr>
    <tr>
      <th>Porosity</th>
      <td>162.0</td>
      <td>0.181001</td>
      <td>0.037196</td>
      <td>0.083842</td>
      <td>0.155735</td>
      <td>0.194823</td>
      <td>0.207754</td>
      <td>0.242298</td>
    </tr>
    <tr>
      <th>Perm</th>
      <td>162.0</td>
      <td>293.798999</td>
      <td>399.989890</td>
      <td>0.381032</td>
      <td>31.996547</td>
      <td>155.888123</td>
      <td>385.177730</td>
      <td>2642.999829</td>
    </tr>
    <tr>
      <th>AI</th>
      <td>162.0</td>
      <td>3441.701689</td>
      <td>969.924695</td>
      <td>1844.166880</td>
      <td>2748.296631</td>
      <td>3179.596509</td>
      <td>3998.567914</td>
      <td>6197.834381</td>
    </tr>
  </tbody>
</table>
</div>




```python
df_shale.describe().transpose()                         # summary table of shale only DataFrame statistics
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>count</th>
      <th>mean</th>
      <th>std</th>
      <th>min</th>
      <th>25%</th>
      <th>50%</th>
      <th>75%</th>
      <th>max</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>index</th>
      <td>99.0</td>
      <td>149.666667</td>
      <td>93.588330</td>
      <td>1.000000</td>
      <td>41.500000</td>
      <td>201.000000</td>
      <td>225.500000</td>
      <td>259.000000</td>
    </tr>
    <tr>
      <th>X</th>
      <td>99.0</td>
      <td>300.444444</td>
      <td>196.365820</td>
      <td>40.000000</td>
      <td>201.000000</td>
      <td>231.000000</td>
      <td>300.000000</td>
      <td>970.000000</td>
    </tr>
    <tr>
      <th>Y</th>
      <td>99.0</td>
      <td>425.111111</td>
      <td>183.665784</td>
      <td>29.000000</td>
      <td>376.000000</td>
      <td>416.000000</td>
      <td>456.000000</td>
      <td>989.000000</td>
    </tr>
    <tr>
      <th>Facies</th>
      <td>99.0</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>0.000000</td>
    </tr>
    <tr>
      <th>Porosity</th>
      <td>99.0</td>
      <td>0.100212</td>
      <td>0.014483</td>
      <td>0.058871</td>
      <td>0.091902</td>
      <td>0.101987</td>
      <td>0.109846</td>
      <td>0.141657</td>
    </tr>
    <tr>
      <th>Perm</th>
      <td>99.0</td>
      <td>3.568463</td>
      <td>6.782476</td>
      <td>0.033611</td>
      <td>0.723788</td>
      <td>1.530878</td>
      <td>3.743835</td>
      <td>52.500870</td>
    </tr>
    <tr>
      <th>AI</th>
      <td>99.0</td>
      <td>5450.493543</td>
      <td>728.871517</td>
      <td>3595.586977</td>
      <td>4897.828718</td>
      <td>5570.604405</td>
      <td>5951.816775</td>
      <td>7881.898531</td>
    </tr>
  </tbody>
</table>
</div>



The facies have strong differences in their summary statistics.

Let's transform the data grouped overall both facies (sand and shale) and separated by facies to normal score values (Gaussian distributed with a mean of 0.0 and variance of 1.0). This is required for sequential Gaussian simulation (common target for our variogram models) and the Gaussian transform assists with outliers and provides more interpretable variograms. 

Let's look at the inputs for the GeostatsPy nscore program.  Note the output include an ndarray with the transformed values (in the same order as the input data in Dataframe 'df' and column 'vcol'), and the transformation table in original values and also in normal score values. 


```python
geostats.nscore                                         # see the input parameters required by the nscore function
```




    <function geostatspy.geostats.nscore(df, vcol, wcol=0, ismooth=0, dfsmooth=0, smcol=0, smwcol=0)>



The following command will transform the Porosity and Permeabilty to standard normal. Note: we perform the transform for porsity and permeability and for each for the cases of:

1. all facies together
2. sand facies only
3. shale facies only


```python
#Transform to Gaussian by Facies
df['NPor'], tvPor, tnsPor = geostats.nscore(df, 'Porosity') # nscore transform for all facies porosity 
df_sand['NPor'], tvPorSand, tnsPorSand = geostats.nscore(df_sand, 'Porosity')  # nscore transform for sand facies porosity 
df_shale['NPor'], tvPorShale, tnsPorShale = geostats.nscore(df_shale, 'Porosity')  # nscore transform for shale facies porosity
df['NPerm'], tvPermSand, tnsPermSand = geostats.nscore(df, 'Perm')  # nscore transform for all facies permeability
df_sand['NPerm'], tvPermSand, tnsPermSand = geostats.nscore(df_sand, 'Perm')  # nscore transform for sand facies permeability
df_shale['NPerm'], tvPermShale, tnsPermShale = geostats.nscore(df_shale, 'Perm')  # nscore transform for shale facies permeability
```

Let's look at a previous of one of the DataFrames to make sure that we now have the normal score porosity and permeability.


```python
df_sand.head()                                          # preview sand DataFrame with nscore transforms
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>index</th>
      <th>X</th>
      <th>Y</th>
      <th>Facies</th>
      <th>Porosity</th>
      <th>Perm</th>
      <th>AI</th>
      <th>NPor</th>
      <th>NPerm</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>100.0</td>
      <td>900.0</td>
      <td>1.0</td>
      <td>0.100187</td>
      <td>1.363890</td>
      <td>5110.699751</td>
      <td>-2.158819</td>
      <td>-2.158819</td>
    </tr>
    <tr>
      <th>1</th>
      <td>8</td>
      <td>100.0</td>
      <td>100.0</td>
      <td>1.0</td>
      <td>0.137453</td>
      <td>5.727603</td>
      <td>5823.241783</td>
      <td>-0.817609</td>
      <td>-1.468475</td>
    </tr>
    <tr>
      <th>2</th>
      <td>9</td>
      <td>200.0</td>
      <td>900.0</td>
      <td>1.0</td>
      <td>0.137062</td>
      <td>14.771314</td>
      <td>5621.146994</td>
      <td>-0.839418</td>
      <td>-0.955141</td>
    </tr>
    <tr>
      <th>3</th>
      <td>10</td>
      <td>200.0</td>
      <td>800.0</td>
      <td>1.0</td>
      <td>0.125984</td>
      <td>10.675436</td>
      <td>4292.700500</td>
      <td>-1.031153</td>
      <td>-1.237102</td>
    </tr>
    <tr>
      <th>4</th>
      <td>18</td>
      <td>300.0</td>
      <td>900.0</td>
      <td>1.0</td>
      <td>0.111516</td>
      <td>27.999817</td>
      <td>4183.466773</td>
      <td>-1.621370</td>
      <td>-0.694044</td>
    </tr>
  </tbody>
</table>
</div>



That looks good! One way to check is to see if the relative magnitudes of the normal score transformed values match the original values.  e.g. that the normal score transform of 0.10 porosity normal score is less than the normal score transform of 0.14 porsity.  Also, the normal score transform of values close to the mean value should be close to 0.0 

Let's also check the original and transformed sand and shale porosity distributions.


```python
plt.subplot(121)                                        # plot original sand and shale porosity histograms
plt.hist(df_sand['Porosity'], facecolor='red',bins=np.linspace(0.0,0.4,50),alpha=0.2,density=True,edgecolor='black',label='Sand')
plt.hist(df_shale['Porosity'], facecolor='blue',bins=np.linspace(0.0,0.4,50),alpha=0.2,density=True,edgecolor='black',label = 'Shale')
plt.xlim([0.05,0.25]); plt.ylim([0,35.0])
plt.xlabel('Porosity (fraction)'); plt.ylabel('Frequency'); plt.title('Porosity Sand and Shale')
plt.legend(loc='upper left')
plt.grid(True)

plt.subplot(122)                                        # plot nscore transformed sand and shale histograms
plt.hist(df_sand['NPor'], facecolor='red',bins=np.linspace(-3.0,3.0,40),histtype="stepfilled",alpha=0.2,density=True,cumulative=False,edgecolor='black',label='Sand')
plt.hist(df_shale['NPor'], facecolor='blue',bins=np.linspace(-3.0,3.0,40),histtype="stepfilled",alpha=0.2,density=True,cumulative=False,edgecolor='black',label='Shale')
plt.xlim([-3.0,3.0]); plt.ylim([0,0.50])
plt.xlabel('Nscore Porosity'); plt.ylabel('Density'); plt.title('Nscore Porosity Sand and Shale')
plt.legend(loc='upper left')
plt.grid(True)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)
plt.show()
```

The normal score transform has correctly transformed the porosity over sand and shale facies to standard normal.  Let's plot the location maps of normal score transforms of porosity and permeability for all facies, sand facies and shale facies.


```python
cmap = plt.cm.plasma                    # color map
plt.subplot(231)
GSLIB.locmap_st(df,'X','Y','NPor',0,1000,0,1000,-3,3,'Nscore Porosity - All Facies','X (m)','Y (m)','Nscore Porosity',cmap)

plt.subplot(232)
GSLIB.locmap_st(df_sand,'X','Y','NPor',0,1000,0,1000,-3,3,'Nscore Porosity - Sand Facies','X (m)','Y (m)','Nscore Porosity',cmap)

plt.subplot(233)
GSLIB.locmap_st(df_shale,'X','Y','NPor',0,1000,0,1000,-3,3,'Nscore Porosity - Shale Facies','X (m)','Y (m)','Nscore Porosity',cmap)

plt.subplot(234)
GSLIB.locmap_st(df,'X','Y','NPerm',0,1000,0,1000,-3,3,'Nscore Permeability - All Facies','X (m)','Y (m)','Nscore Permeability',cmap)

plt.subplot(235)
GSLIB.locmap_st(df_sand,'X','Y','NPerm',0,1000,0,1000,-3,3,'Nscore Permeability - Sand Facies','X (m)','Y (m)','Nscore Permeability',cmap)

plt.subplot(236)
GSLIB.locmap_st(df_shale,'X','Y','NPerm',0,1000,0,1000,-3,3,'Nscore Permeability - Shale Facies','X (m)','Y (m)','Nscore Permeability',cmap)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.8, wspace=0.4, hspace=0.3)
plt.show()
```

Let's see what the parameters are for the gamv, irregular data, experimental variogram calculation program.


```python
geostats.gamv                                           # see the input parameters required by the gamv function
```




    <function geostatspy.geostats.gamv(df, xcol, ycol, vcol, tmin, tmax, xlag, xltol, nlag, azm, atol, bandwh, isill)>



We can use the location maps to help determine good variogram calculation parameters.

We are ready to calculate variogram! Let's calculate isotropic variograms for the transformed normal score porosity and permeability for sand, shale and mixed (without separating sand and shale).  Some information on the parameters that I chose:

```p
tmin = -9999.; tmax = 9999.; 
lag_dist = 100.0; lag_tol = 50.0; nlag = 7; bandh = 9999.9; azi = 0; atol = 90.0; isill = 1
```
* tmin, tmax are trimming limits - set to have no impact, no need to filter the data
* lag_dist, lag_tol are the lag distance, lag tolerance - set based on the common data spacing (100m) and tolerance as 100% of lag distance for additonal smoothing
* nlag is number of lags - set to extend just past 50 of the data extent
* bandh is the horizontal band width - set to have no effect
* azi is the azimuth -  it has not effect since we set atol, the azimuth tolerance, to 90.0
* isill is a boolean to standardize the distribution to a variance of 1 - it has no effect since the nscore transform sets the variance to 1.0

Let's try running these variograms and visualizing them.


```python
# Calculate Sample Data Isotropic Variograms
tmin = -9999.; tmax = 9999.; 
lag_dist = 100.0; lag_tol = 100.0; nlag = 7; bandh = 9999.9; azi = 0; atol = 90.0; isill = 1

lag, por_sand_gamma, por_sand_npair = geostats.gamv(df_sand,"X","Y","NPor",tmin,tmax,lag_dist,lag_tol,nlag,azi,atol,bandh,isill)
lag, por_shale_gamma, por_shale_npair = geostats.gamv(df_shale,"X","Y","NPor",tmin,tmax,lag_dist,lag_tol,nlag,azi,atol,bandh,isill)
lag, por_gamma, por_npair = geostats.gamv(df,"X","Y","NPor",tmin,tmax,lag_dist,lag_tol,nlag,azi,atol,bandh,isill)

lag, perm_sand_gamma, perm_sand_npair = geostats.gamv(df_sand,"X","Y","NPerm",tmin,tmax,lag_dist,lag_tol,nlag,azi,atol,bandh,isill)
lag, perm_shale_gamma, perm_shale_npair = geostats.gamv(df_shale,"X","Y","NPerm",tmin,tmax,lag_dist,lag_tol,nlag,azi,atol,bandh,isill)
lag, perm_gamma, perm_npair = geostats.gamv(df,"X","Y","NPerm",tmin,tmax,lag_dist,lag_tol,nlag,azi,atol,bandh,isill)

plt.subplot(121)
plt.plot(lag,por_gamma,'x',color = 'black',label = 'All')
plt.plot(lag,por_sand_gamma,'x',color = 'orange',label = 'Sand')
plt.plot(lag,por_shale_gamma,'x',color = 'brown',label = 'Shale')
plt.plot([0,2000],[1.0,1.0],color = 'black')
plt.xlabel(r'Lag Distance $\bf(h)$, (m)')
plt.ylabel(r'$\gamma \bf(h)$')
plt.title('Isotropic NSCORE Porosity Variogram')
plt.xlim([0,700])
plt.ylim([0,1.8])
plt.legend(loc='upper left')
plt.grid(True)

plt.subplot(122)
plt.plot(lag,perm_gamma,'x',color = 'black',label = 'All')
plt.plot(lag,perm_sand_gamma,'x',color = 'orange',label = 'Sand')
plt.plot(lag,perm_shale_gamma,'x',color = 'brown',label = 'Shale')
plt.plot([0,2000],[1.0,1.0],color = 'black')
plt.xlabel(r'Lag Distance $\bf(h)$, (m)')
plt.ylabel(r'$\gamma \bf(h)$')
plt.title('Isotropic NSCORE Permeaiblity Variogram')
plt.xlim([0,700])
plt.ylim([0,1.8])
plt.legend(loc='upper left')
plt.grid(True)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.3)
plt.show()
```

The experiental variograms have some interesting features:

* the range of the sand porosity and permeability is greater than the shale porosity and permeability range
* although the shale short range experimental points may be nosy due to sparce shale data

There is much more that you could try out, e.g. direction variograms.

I hope you find this code and demonstration useful. I'm always happy to discuss geostatistics, statistical modeling, uncertainty modeling and machine learning,

*Michael*

**Michael Pyrcz**, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin
On Twitter I'm the **GeostatsGuy** and on YouTube my lectures are on the channel, **GeostatsGuy Lectures**.

c1230fa1419ee35e570ff10945640d687c734bc5
ipynb
Jupyter Notebook
GeostatsPy_variogram_calculation.ipynb
AndrewAnnex/PythonNumericalDemos
11a7dec09bf08527b358fa95119811b6f73023b5
MIT
2019-03-18T10:35:08.000Z
2021-12-15T05:49:00.000Z
GeostatsPy_variogram_calculation.ipynb
AndrewAnnex/PythonNumericalDemos
11a7dec09bf08527b358fa95119811b6f73023b5
MIT
GeostatsPy_variogram_calculation.ipynb
AndrewAnnex/PythonNumericalDemos
11a7dec09bf08527b358fa95119811b6f73023b5
MIT
2019-07-09T10:25:21.000Z
2021-04-28T21:56:01.000Z
Qwen/Qwen-72B
1. YES
2. YES