<jupyter_start><jupyter_text>Daytum +2 Course: Data Analytics, Geostatistics and Machine Learning Deep Dive Trend Modeling Demonstration and Exercise  GoalCalculate data-driven trend for a spatial dataset. DescriptionHere's a simple, documented workflow, demonstration of trend calculation and diagnostics for subsurface modeling workflows. This should help you get started with building subsurface models that integrate trends.   Trend ModelingTrend modeling is the modeling of local features, based on data and interpretation, that are deemed certain (known).  The trend is substracted from the data, leaving a residual that is modeled stochastically with uncertainty (treated as unknown).* geostatistical spatial estimation methods will make an assumption concerning stationarity    * in the presence of significant nonstationarity we can not rely on spatial estimates based on data + spatial continuity model* if we observe a trend, we should model the trend.    * then model the residuals stochastically    Trend modeling significantly improves our models by accounting for the nonstationary map-able component of the spatial property.  Consider this spatial model in 1D with and without trend.INSERT FIGURE FROM 04_Spatial_Continuity.ppts SLIDE 50.At a distance beyond the range of spatial continuity away from the data, the best estimate approaches the global mean or the trend model.  WIthout the trend model we would systemally underestimate in the high regions and overestimate in the low regions.Steps: 1. model trend consistent with data and intepretation at all locations within the area of itnerest, integrate all available information and expertise.\begin{equation}m(\bf{u}_\beta), \, \forall \, \beta \in \, AOI\end{equation}2. substract trend from data at the $n$ data locations to formulate a residual at the data locations.\begin{equation}y(\bf{u}_{\alpha}) = z(\bf{u}_{\alpha}) - m(\bf{u}_{\alpha}), \, \forall \, \alpha = 1, \ldots, n\end{equation}3. characterize the statistical behavoir of the residual $y(\bf{u}_{\alpha})$ integrating any information sources and interpretations.  For example the global cumulative distribution function and a measure of spatial continuity shown here.\begin{equation}F_y(y) \quad \gamma_y(\bf{h})\end{equation}4. model the residual at all locations with $L$ multiple realizations.\begin{equation}Y^\ell(\bf{u}_\beta),  \, \forall \, \beta \, \in \, AOI; \, \ell = 1, \ldots, L\end{equation}5. add the trend back in to the stochastic residual realizations to calculate the multiple realizations, $L$, of the property of interest based on the composite model of known deterministic trend, $m(\bf{u}_\alpha)$ and unknown stochastic residual, $y(\bf{u}_\alpha)$ \begin{equation}Z^\ell(\bf{u}_\beta) = Y^\ell(\bf{u}_\beta) + m(\bf{u}_\beta),  \, \forall \, \beta \in \, AOI; \, \ell = 1, \ldots, L\end{equation}6. check the model, including quantification of the proportion of variance treated as known (trend) and unknown (residual).\begin{equation}\sigma^2_{Z} = \sigma^2_{Y} + \sigma^2_{m} + 2 \cdot C_{Y,m}\end{equation}given $C_{Y,m} \to 0$:\begin{equation}\sigma^2_{Z} = \sigma^2_{Y} + \sigma^2_{m}\end{equation}I can now describe the proportion of variance allocated to known and unknown components as follows:\begin{equation}Prop_{Known} = \frac{\sigma^2_{m}}{\sigma^2_{Y} + \sigma^2_{m}} \quad Prop_{Unknown} = \frac{\sigma^2_{Y}}{\sigma^2_{Y} + \sigma^2_{m}}\end{equation}I provide some practical, data-driven methods for trend model, but I should indicate that:1. trend modeling is very important in reservoir modeling as it has large impact on local model accuracy and on the undertainty model2. trend modeling is used in almost every subsurface model, unless the data is dense enough to impose local trends3. trend modeling includes a high degree of expert judgement combined with the integration of various information sourcesWe limit ourselves to simple data-driven methods, but acknowledge much more is needed. In fact, trend modeling requires a high degree of knowledge concerning local geoscience and engineering data and knowledge.   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 StartedHere'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_biased.csv at https://git.io/fh0CWThere 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.<jupyter_code>import geostatspy.GSLIB as GSLIB          # GSLIB utilies, visualization and wrapper
import geostatspy.geostats as geostats    # GSLIB methods convert to Python<jupyter_output><empty_output><jupyter_text>We will also need some standard packages. These should have been installed with Anaconda 3.<jupyter_code>import numpy as np                        # ndarrys for gridded data
import pandas as pd                       # DataFrames for tabular data
import os                                 # set working directory, run executables
import matplotlib.pyplot as plt           # for plotting
from scipy import stats                   # summary statistics
import math                               # trig etc.
import scipy.signal as signal             # kernel for moving window calculation<jupyter_output><empty_output><jupyter_text>Set the working directoryI 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).<jupyter_code>os.chdir("c:/PGE383")             # set the working directory<jupyter_output><empty_output><jupyter_text>Loading Tabular DataHere's the command to load our comma delimited data file in to a Pandas' DataFrame object.<jupyter_code>df = pd.read_csv('sample_data_biased.csv')     # load our data table (wrong name!)<jupyter_output><empty_output><jupyter_text>It worked, we loaded our file into our DataFrame called 'df'. But how do you really know that it worked? Visualizing the DataFrame would be useful and we already leard about these methods in this demo (https://git.io/fNgRW). We can preview the DataFrame by printing a slice or by utilizing the 'head' DataFrame member function (with a nice and clean format, see below). With the slice we could look at any subset of the data table and with the head command, add parameter 'n=13' to see the first 13 rows of the dataset.<jupyter_code>print(df.iloc[0:5,:])                   # display first 4 samples in the table as a preview
df.head(n=13)                           # we could also use this command for a table preview<jupyter_output>X    Y  Facies  Porosity       Perm
0  100  900       1  0.115359   5.736104
1  100  800       1  0.136425  17.211462
2  100  600       1  0.135810  43.724752
3  100  500       0  0.094414   1.609942
4  100  100       0  0.113049  10.886001<jupyter_text>Summary Statistics for Tabular DataThe table includes X and Y coordinates (meters), Facies 1 and 2 (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.<jupyter_code>df.describe().transpose()<jupyter_output><empty_output><jupyter_text>Visualizing Tabular Data with Location Maps It is natural to set the x and y coordinate and feature ranges manually. e.g. do you want your color bar to go from 0.05887 to 0.24230 exactly? Also, let's pick a color map for display. I heard that plasma is known to be friendly to the color blind as the color and intensity vary together (hope I got that right, it was an interesting Twitter conversation started by <NAME> from Agile if I recall correctly). We will assume a study area of 0 to 1,000m in x and y and omit any data outside this area.<jupyter_code>xmin = 0.0; xmax = 1000.0               # range of x values
ymin = 0.0; ymax = 1000.0               # range of y values
pormin = 0.05; pormax = 0.25;           # range of porosity values
nx = 100; ny = 100; csize = 10.0
cmap = plt.cm.plasma                    # color map<jupyter_output><empty_output><jupyter_text>Let's try out locmap. This is a reimplementation of GSLIB's locmap program that uses matplotlib. I hope you find it simpler than matplotlib, if you want to get more advanced and build custom plots lock at the source. If you improve it, send me the new code. Any help is appreciated. To see the parameters, just type the command name:<jupyter_code>GSLIB.locmap<jupyter_output><empty_output><jupyter_text>Now we can populate the plotting parameters and visualize the porosity data.<jupyter_code>plt.subplot(111)
GSLIB.locmap_st(df,'X','Y','Porosity',xmin,xmax,ymin,ymax,pormin,pormax,'Well Data - Porosity','X(m)','Y(m)','Porosity (fraction)',cmap)
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Let's get some declustering weights. For more information see the demonstration on declustering.<jupyter_code>wts, cell_sizes, dmeans = geostats.declus(df,'X','Y','Porosity',iminmax = 1, noff= 10, ncell=100,cmin=10,cmax=2000)
df['Wts'] = wts                            # add weights to the sample data DataFrame
df.head()                                  # preview to check the sample data DataFrame<jupyter_output>There are 289 data with:
   mean of      0.13474387540138408 
   min and max  0.058547873 and 0.228790002
   standard dev 0.03767982164385207<jupyter_text>Trend by Convolution / Local Window AverageLet's first attempt a convolution-based trend model, this is a moving window average of the local data.We have a convenience function that takes data with X and Y locations in a DataFrame and makes a sparse 2D array.  All cells without a data value are assigned to NumPy's NaN (null values, missing value). Let's see the inputs for this command.<jupyter_code>GSLIB.DataFrame2ndarray<jupyter_output><empty_output><jupyter_text>Let's make an sparse array with the appropriate parameters. The reason we are doing this is that convolution programs in general work with ndarrays and not from DataFrames.<jupyter_code>por_grid = GSLIB.DataFrame2ndarray(df,'X','Y','Porosity',xmin, xmax, ymin, ymax, csize)<jupyter_output><empty_output><jupyter_text>We have a ndarray (por_grid) with the data assigned to grid cells. Now we need a kernel. The kernel represents the weights within the moving window. If we use constant 1.0 in the moving window will get discontinuities in our trend model. A Gaussian kernel (weights highest in the middle of the window and decreasing to 0,0 at the edge) is useful to get a smooth trend. We can use the SciPy package's signal functions. Of course we will have to import that package and then we can make our kernel. Here's an example below. There shouldnt be any surprises.<jupyter_code>gkern1d = signal.gaussian(53,15).reshape(53, 1)
gkern2d = np.outer(gkern1d, gkern1d)
print('We have made a kernel of size, number of grid cells (ny, nx) ' + str(gkern2d.shape))

plt.subplot(111)
GSLIB.pixelplt_st(gkern2d,xmin=-265,xmax=265,ymin=-265,ymax=265,step=10,vmin=0,vmax=1,title='Kernel',xlabel='X(m)',ylabel='Y(m)',vlabel='weight',cmap=cmap)
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.6, top=0.8, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output>We have made a kernel of size, number of grid cells (ny, nx) (53, 53)<jupyter_text>Now we need to convolve our sparse data assigned to a ndarray with our Gaussian kernel. There are many functions available for convolution. But we have a problem as we want to apply our Gaussian kernel to a sparse ndarray full of missing values. It turns out this is a common issue for our friends in Astronomy and so their Astropy package has a convolution method that will work well. I figured out the following (so you don't have to!).<jupyter_code>import astropy.convolution.convolve as convolve
porosity_trend = convolve(por_grid,gkern2d,boundary='extend',nan_treatment='interpolate',normalize_kernel=True)<jupyter_output><empty_output><jupyter_text>No errors? It worked? Let's look at the results. We can plot and compare the original porosity data and the resulting trend to check for consistency.<jupyter_code>plt.subplot(131)
GSLIB.locmap_st(df,'X','Y','Porosity',xmin,xmax,ymin,ymax,pormin,pormax,'Well Data - Porosity','X(m)','Y(m)','Porosity (fraction)',cmap)

plt.subplot(132)
GSLIB.pixelplt_st(porosity_trend,xmin,xmax,ymin,ymax,csize,pormin,pormax,'Porosity Trend','X(m)','Y(m)','Porosity (fraction)',cmap)

plt.subplot(133)
GSLIB.locpix_st(porosity_trend,xmin,xmax,ymin,ymax,csize,pormin,pormax,df,'X','Y','Porosity','Porosity Data and Trend','X(m)','Y(m)','Porosity (fraciton)',cmap)

plt.subplots_adjust(left=0.0, bottom=0.0, right=3.0, top=1.2, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Other Methods for Trend CalculationThere are a variety of other methods for trend calculation. I will just mention them here.1. hand-drawn, expert interpretation - many 3D modeling packages allow for experts to draw trends and allow for fast interpolation to build an exhaustive trend model.2. kriging - kriging provides best linear unbiased estimates between data given a spatial continuity model (more on this when we cover spatial estimation).  One note of caution is that kriging is exact; therefore it will over fit unless it is use with averaged data values (e.g. over the vertical) or with a block kriging option (kriging at a volume support larger than the data).3. regression - fit a function as a function of X, Y coordinates. This could be extended to more complicated prediction models from machine learning.  Adding Trend to DataFrameLet's add the porosity trend to our DataFrame. We have a sample program in GeostatsPy that takes a 2D ndarray and extracts the values at the data locations and adds them as a new column.  Then we can do a little math to calculate and add the porosity residual also and visualize this all together as a final check.<jupyter_code>df = GSLIB.sample(porosity_trend,xmin,ymin,csize,"Por_Trend",df,'X','Y')
df['Por_Res'] = df['Porosity'] - df['Por_Trend'] # calculate the residual and add to DataFrame<jupyter_output><empty_output><jupyter_text>Let's check out the DataFrame and confirm that we have everything now. We will need trend and residual in our DataFrame to support all subsequent modeling steps.<jupyter_code>df.head()<jupyter_output><empty_output><jupyter_text>Trend DiagnotisticsLet's go back to the convolution trend and check it (to demonstrate the method of trend checking). Note, I haven't tried to perfect the result. I'm just demonstrating the method. In addition to the previous visualization, let's look at the distributions and summary statistics of the original declustered porosity data and the trend.<jupyter_code>plt.subplot(131)
GSLIB.hist_st(df['Porosity'],pormin,pormax,False,False,20,df['Wts'],'Porosity (fraction)','Declustered Porosity')

plt.subplot(132)
GSLIB.hist_st(porosity_trend.flatten(),pormin,pormax,False,False,20,None,'Porosity Trend (fraction)','Porosity Trend')

plt.subplot(133)
GSLIB.hist_st(df['Por_Res'],-.1,.1,False,False,20,df['Wts'],'Porosity Residual (fraction)','Declustered Porosity Residual')

plt.subplots_adjust(left=0.0, bottom=0.0, right=3.0, top=1.5, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>We can also look at the summary statistics. Here's a function that calculates the weighted standard deviation (and the average). We can use this with the data and declustering weights and figure out the allocation of variance between the trend and the residual.  We will calculated the proportion of variance in the trend and the residual.  This includes the following steps:* calculate the representative weighted standard deviation of the original data, convert to variance, $\sigma$* calculate the standard deviation of the trend model, convert to variance, $\sigma_t$* calculate the variance of the residual at the data locations, $\sigma_r$ * back calculate the covariance form the relationship:\begin{equation}\sigma = \sigma_t + \sigma_r + 2 \times C_{tr}\end{equation}<jupyter_code># Weighted average and standard deviation
def weighted_avg_and_std(values, weights):                  # from <NAME>, stack overflow
    average = np.average(values, weights=weights)
    variance = np.average((values-average)**2, weights=weights)
    return (average, math.sqrt(variance))

wavg_por,wstd_por = weighted_avg_and_std(df['Porosity'],df['Wts']) 
var_por_trend = np.var(porosity_trend.flatten())
wavg_por_res,wstd_por_res = weighted_avg_and_std(df['Por_Res'],df['Wts']) 

prop_trend = var_por_trend/(wstd_por**2)
prop_resid = (wstd_por_res**2)/(wstd_por**2)
prop_cov = 1 - prop_trend - prop_resid

variance_prop = [prop_trend, prop_resid, prop_cov]
labels = 'Trend', 'Residual', '2 x Covariance'
colors = 'orange', 'red', 'green'
explode = (0.1, 0.1,0.1)

fig1, ax1 = plt.subplots()
ax1.pie(variance_prop, explode=explode, colors = colors, labels=labels, autopct='%1.1f%%',
        shadow=True, startangle=90)
ax1.axis('equal')  # Equal aspect ratio ensures that pie is drawn as a circle.

plt.show()

print('Declustered Porosity Data: Average ' + str(round(wavg_por,4)) + ', Var ' + str(round(wstd_por**2,5)))
print('Porosity Trend: Average            ' + str(round(wavg_por_trend,4)) + ', Var ' + str(round(wstd_por_trend**2,5)))
print('Porosity Residual: Average         ' + str(round(wavg_por_res,4)) + ', Var ' + str(round(wstd_por_res**2,5)))
print('Proportion Trend / Known:          ' + str(round(wstd_por_trend**2/(wstd_por**2),3)))
print('Proportion Residual / Unknown:     ' + str(round((wstd_por**2 - wstd_por_trend**2)/(wstd_por**2),3)))<jupyter_output><empty_output><jupyter_text>When we analyze the trend contributions we observe that:* 48% of the variance is captured by the trend model* 41% of the variance is captured in the residual* 11% of the vairance remains in the coveriance termThat looks good. A quick check, confirm that the Porosity column is equal to the Por_Trend + the Por_Res columns.  As a final check let's visualize the original porosity data, porosity trends at the data locations and the porosity residuals.<jupyter_code>plt.subplot(131)
GSLIB.locmap_st(df,'X','Y','Porosity',xmin,xmax,ymin,ymax,pormin,pormax,'Well Data - Porosity','X(m)','Y(m)','Porosity (fraction)',cmap)

plt.subplot(132)
GSLIB.locmap_st(df,'X','Y','Por_Trend',xmin,xmax,ymin,ymax,pormin,pormax,'Well Data - Porosity Trend','X(m)','Y(m)','Porosity (fraction)',cmap)

plt.subplot(133)
GSLIB.locmap_st(df,'X','Y','Por_Res',xmin,xmax,ymin,ymax,-0.01,0.01,'Well Data - Porosity Residual','X(m)','Y(m)','Porosity (fraction)',cmap)

plt.subplots_adjust(left=0.0, bottom=0.0, right=3.0, top=1.2, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Does it look correct? There is a strong degree of consistency between the porosity data and trend and the porosity residual no longer has a trend, it has been detrended.We should also check for conditional bias.  Let's look at the porosity residuals vs. the porosity trend.<jupyter_code>plt.subplot(111)
plt.scatter(df['Por_Trend'],df['Por_Res'], s=30, alpha = 0.2, edgecolors = "black", facecolors = 'red')
plt.xlabel('Porosity Trend (fraction)')
plt.ylabel('Porosity Residual (residual fraction)')
plt.title('Porosity Residual vs. Porosity Trend')
plt.plot([0,2000],[0,0],color = 'black')
plt.ylim(-0.1,0.1)
plt.xlim(0.075,0.2)
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.2, top=1.2, wspace=0.2, hspace=0.2)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09332503351848305, 'fraction_numerical': 0.025186745834131392, 'mean_word_length': 4.946753986332574, 'pattern_counts': {'":': 0, '<': 84, '<?xml version=': 0, '>': 84, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '427318', 'license': 'Apache-2.0', 'n_tokens_mistral': 6679, 'n_tokens_neox': 6151, 'n_words': 2405, 'path': '2Day/Daytum_Trends_v2.ipynb', 'repo_name': 'GeostatsGuy/curriculum', 'size': '20895'}