from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
from ipywidgets import GridspecLayout
import ipywidgets as widgets
from IPython.display import display
from IPython.display import clear_output
from plotly.subplots import make_subplots
import plotly.express as px
import plotly.offline as pyo
import plotly.io as pio
import plotly.express as px
import plotly.graph_objects as go
import matplotlib.pyplot as plt
from matplotlib import colors,ticker,cm
import matplotlib
from sklearn.neighbors import KNeighborsRegressor
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import fitpack
from scipy.interpolate import griddata
from scipy import interpolate
import numpy as np
import pandas as pd
import time
import math
import numba
import warnings
import math
import os
from pyLPM import plots
# INITIALIZE OFFLINE NOTEBOOK MODE
pyo.init_notebook_mode()
##############################################################################################################################################
# global variables
global return_exp_var # return of experimental variograms
return_exp_var = []
global return_model_var #return of model variograms
return_model_var = {}
##############################################################################################################################################
# Calculate the distance coordinates
@numba.jit(fastmath=False)
def _hdist(distancex, distancey, distancez):
"""euclidian distance between samples
Args:
distancex (np.array): distances in x coordinates
distancey (np.array): distances in y coordinates
distancez (np.array): distances in z coordinantes
Returns:
np.array: list of euclidian distances between samples
"""
dist =np.zeros(distancex.shape[0])
for i in range(distancex.shape[0]):
dist[i] = np.sqrt(distancex[i]**2 + distancey[i]**2 + distancez[i]**2) + 0.0000000001
return dist
##############################################################################################################################################
# Calculate the distances in planar coordinates
@numba.jit(fastmath=False)
def _xydist(distancex, distancey):
"""Planar distances between samples
Args:
distancex (np.array): distances in x coordinates
distancey (np.array): distances in y coordinates
Returns:
np.array: List of planar distances
"""
dist =np.zeros(distancex.shape[0])
for i in range(distancex.shape[0]):
dist[i] = np.sqrt(distancex[i]**2 + distancey[i]**2 ) + 0.0000000001
return dist
##############################################################################################################################################
# Calculate the x coordinates distances
@numba.jit(fastmath=False)
def _xdist(pairs):
"""Summary
Args:
pairs (lst): Sample points
Returns:
np.array: x coordinate distances
"""
dist =np.zeros(pairs.shape[0])
for i in range(pairs.shape[0]):
dist[i] = (pairs[i][0][0] - pairs[i][1][0])
return dist
##############################################################################################################################################
# Calculate the y coordinates distances
@numba.jit(fastmath=False)
def _ydist(pairs):
"""Summary
Args:
pairs (lst): Sample points
Returns:
np.array: y coordinate distances
"""
dist =np.zeros(pairs.shape[0])
for i in range(pairs.shape[0]):
dist[i] = (pairs[i][0][1] - pairs[i][1][1])
return dist
##############################################################################################################################################
# Calculate the z coordinates distances
@numba.jit(fastmath=False)
def _zdist(pairs):
"""Summary
Args:
pairs (lst): Sample points
Returns:
np.array: z coordinate distances
"""
dist =np.zeros(pairs.shape[0])
for i in range(pairs.shape[0]):
dist[i] = (pairs[i][0][2] - pairs[i][1][2])
return dist
##############################################################################################################################################
# Calculate pairs combination where vertical and horizontal distances are less than the maximum distance
@numba.jit(fastmath=False)
def _combin(points,n, max_dist):
"""Summary
Args:
points (lst): sample points
n (integer): number of samples
max_dist (float): maximum permissible distance
Returns:
lst: List of tuples containing the permissible pairs
"""
dist =[]
p = 0
for i in range(0,n):
for j in range((i+1),n):
if (points[i][0] - points[j][0]) < max_dist:
if (points[i][1] - points[j][1]) < max_dist:
if (points[i][2] - points[j][2]) < max_dist:
dist.append((points[i] , points[j]))
p += 1
return np.array(dist)
##############################################################################################################################################
# Rotate data according azimuth and dip directions
@numba.jit(fastmath=False)
def _rotate_data(xh, yh, zh, azimute, dip):
"""Summary
Args:
xh (lst): x coordinate distances
yh (lst): y coordinate distances
zh (lst): z coordinate distances
azimute (lst): azimuth directions
dip (lst): dip directions
Returns:
lst: Rotate coordiantes Xrot, Yrot and Zrot
"""
xhrot = []
yhrot = []
zhrot = []
for i in range(0, len(xh)):
# ROTACIONE PRIMEIRAMENTE NO AZIMUTE
xrot = math.cos(math.radians(azimute))*xh[i] - math.sin(math.radians(azimute))*yh[i]
yrot = math.sin(math.radians(azimute))*xh[i] + math.cos(math.radians(azimute))*yh[i]
yrot2 = math.cos(math.radians(dip))*yrot - math.sin(math.radians(dip))*zh[i]
zrot = math.sin(math.radians(dip))*yrot + math.cos(math.radians(dip))*zh[i]
xhrot.append(xrot)
yhrot.append(yrot2)
zhrot.append(zrot)
return xhrot, yhrot, zhrot
##############################################################################################################################################
# Calculate all distances in dataset
def _distances(dataset, nlags, x_label, y_label, z_label,
lagdistance, head_property, tail_property, choice):
'''distances
Args:
dataset (pd.DataFrame): pandas dataset containing all spatial data
nlags (int): Number of lags
x_label (string): Label for X direction
y_label (string): Label for Y direction
z_label (string): Label for Z direction
lagdistance (float): Lag size
head_property (string): Label for head variable
tail_property (string): Label for the tail variable
choice (float): Percentual of datasest sampling 0.0-1.0
Returns:
pd.DataFrame: Distances of permissible pairs
'''
# Variables definition
max_dist = (nlags + 1)*lagdistance # Define the maximum distance search
X = dataset[x_label].values
Y = dataset[y_label].values
Z = dataset[z_label].values
HEAD = dataset[head_property].values
TAIL = dataset[tail_property].values
# If random option is selected, select a number of data according some proportion
if (choice != 1.0):
points = np.array(list(zip(X,Y,Z,HEAD,TAIL)))
points = np.array([points[np.random.randint(0,len(points))] for i in range(int(points.shape[0]*choice))])
else:
points = np.array(list(zip(X,Y,Z,HEAD,TAIL)))
# Define the permissible combination of samples according the maximum distance
pairs = _combin(points,points.shape[0], max_dist)
# Define distance variables and dataset
distancex = _xdist(pairs)
distancey = _ydist(pairs)
distancez = _zdist(pairs)
distanceh = _hdist(distancex, distancey, distancez)
distancexy = _xydist(distancex, distancey)
head_1 = np.array([pair[0][3] for pair in pairs])
head_2 = np.array([pair[0][4] for pair in pairs])
tail_1 = np.array([pair[1][3] for pair in pairs])
tail_2 = np.array([pair[1][4] for pair in pairs])
# Return the distances dataframe
distance_dataframe = np.array([distancex,
distancey,
distancez,
distancexy,
distanceh,
head_1,
head_2,
tail_1,
tail_2,]).T
return distance_dataframe[distanceh[:,] < max_dist ] # Only values less than the maximum distance
##############################################################################################################################################
def _distances_varmap(dataset, X,Y,Z, nlags, lagdistance, head_property, tail_property,
choice):
'''distances
Args:
dataset (pd.DataFrame): DataFrame containing all spatial data
X (string): Label for the X coordinates
Y (string): Label for the Y coordinates
Z (string): Label for the Z coordinates
nlags (int): Number of lags to calculate distannces
lagdistance (float): Size of the lag
head_property (string): Label for the head property
tail_property (string): Label for the tail property
choice (float): Perncentual of data sample
Return:
np.array: All permissible dataset distances
'''
# Define coordinates
X = np.array(X)
Y = np.array(Y)
Z = np.array(Z)
# Define the maximum distance
max_dist = (nlags + 1)*lagdistance
# Define tail and head properties
HEAD = dataset[head_property].values
TAIL = dataset[tail_property].values
# If random option is selected, select a number of data according some proportion
if (choice != 1.0):
np.random.seed(0)
np.random.seed(138276)
points = np.array(list(zip(X,Y,Z,HEAD,TAIL)))
points = np.array([points[np.random.randint(0,len(points))] for i in range(int(points.shape[0]*choice))])
else:
points = np.array(list(zip(X,Y,Z,HEAD,TAIL)))
# Define the permissible combination of samples according the maximum distance
pairs = _combin(points,points.shape[0], max_dist)
# Define distance variables and dataset
distancex = _xdist(pairs)
distancey = _ydist(pairs)
distancez = _zdist(pairs)
distanceh = _hdist(distancex, distancey, distancez)
distancexy = _xydist(distancex, distancey)
head_1 = np.array([pair[0][3] for pair in pairs])
head_2 = np.array([pair[0][4] for pair in pairs])
tail_1 = np.array([pair[1][3] for pair in pairs])
tail_2 = np.array([pair[1][4] for pair in pairs])
# Return the distances dataframe
distance_dataframe = np.array([distancex,
distancey,
distancez,
distancexy,
distanceh,
head_1,
head_2,
tail_1,
tail_2,]).T
return distance_dataframe[distanceh[:,] < max_dist ] # Only values less than the maximum distance
##############################################################################################################################################
def __permissible_pairs (lag_multiply, lagdistance, lineartolerance, check_azimuth,
check_dip, check_bandh, check_bandv, htol, vtol, hband, vband, omni,
dist):
'''permissible_pairs
Args:
lag_multiply (double): Mutliplemaxi of lag distance
lagdistance (TYPE): Description
lineartolerance (TYPE): Description
check_azimuth (TYPE): Description
check_dip (TYPE): Description
check_bandh (TYPE): Description
check_bandv (TYPE): Description
htol (TYPE): Description
vtol (TYPE): Description
hband (TYPE): Description
vband (TYPE): Description
omni (TYPE): Description
dist (TYPE): Description
Returns:
distances (numpy array): Returns the permissible sample pairs for omnidirecional functions
'''
# Define the minimum range and the maximum rannge
minimum_range = lag_multiply*lagdistance - lineartolerance
maximum_range = lag_multiply*lagdistance + lineartolerance
# Filter samples acoording distances tolerances
if omni == False:
filter_dist = dist[(dist[:,4] >= minimum_range) &
(dist[:,4] <= maximum_range) &
(check_azimuth >= htol) &
(check_dip >= vtol) &
(check_bandh < hband)&
(check_bandv < vband)]
else:
filter_dist = dist[(dist[:,4] >= minimum_range) &
(dist[:,4] <= maximum_range)]
return filter_dist
##############################################################################################################################################
def __calculate_experimental(dist, lag_multiply, lagdistance, lineartolerance, check_azimuth,
check_dip, check_bandh, check_bandv, htol, vtol, hband, vband, omni,type_var):
'''calculate_experimental
Args:
dist (pd.DataFrame): Pandas Dataframe containing all spatial data
lag_multiply (int): Number of lag mutliplication
lagdistance (float): Lag distance size
lineartolerance (float): Linear tolerance length
check_azimuth (float): check if azimuth direction is permissible
check_dip (float): check if dip direction is permissible
check_bandh (float): check if horizontal band is permissible
check_bandv (float): check if vertical band is permissible
htol (float): Horizontal angular tolerance
vtol (float): Vertical angular tolerance
hband (float): Horizontal band width
vband (float): Vertical band width
omni (bool): Boolean for omnidirecional variogram (True= use, False =don`t use)
type_var (string): Experimental Continuity function
type_var:
'Variogram'
'Covariogram'
'Correlogram'
'PairWise'
'Relative_Variogram'
Returns:
value (double): Experimental continuity function value
number_of_pairs (int): number of pairs used to calculate the experimental function
average_distace (double): average distance of experimental continuity function value
Raises:
Exception: Experimental continuity function not in permissible functions
'''
# Calculate the permissible pairs for experimental continuity functions
points = __permissible_pairs(lag_multiply, lagdistance, lineartolerance, check_azimuth,
check_dip, check_bandh, check_bandv, htol, vtol, hband, vband, omni,
dist)
# Show error if experimental variogram is not valid
if type_var not in['Variogram','Covariogram','Correlogram','PairWise','Relative_Variogram']:
raise Exception("Experimental continuity function not in admissible functions")
# Calculate experimental variogram according the specific experimental function
if points.size != 0:
number_of_pairs = float(points.shape[0])
average_distance = points[:,4].mean()
value_exp = 0
if type_var == 'Variogram':
value_exp = ((points[:,5] - points[:,7])*(points[:,6] - points[:,8]))/(2*number_of_pairs)
value_exp = value_exp.sum()
if type_var == 'Covariogram':
value_exp = ((points[:,5] - points[:,5].mean())*(points[:,8]-points[:,8].mean()))/number_of_pairs
value_exp = value_exp.sum()
if type_var == 'Correlogram':
value_exp = ((points[:,5] - points[:,5].mean())*(points[:,8]-points[:,8].mean()))/(number_of_pairs*points[:,5].var()*points[:,8].var())
value_exp = value_exp.sum()
if type_var == 'PairWise':
value_exp = 2*((points[:,5] - points[:,7])/(points[:,6] + points[:,8]))**2/number_of_pairs
value_exp = value_exp.sum()
if type_var == 'Relative_Variogram':
average_tail = (points[:,7] + points[:,8])/2
average_head = (points[:,5] + points[:,6])/2
value_exp = 4*((points[:,5] - points[:,7])*(points[:,6] - points[:,8]))/(number_of_pairs*(average_head + average_tail)**2)
value_exp = value_exp.sum()
return [value_exp, number_of_pairs, average_distance]
return [np.nan , np.nan, np.nan]
##############################################################################################################################################
def _calculate_experimental_function(dataset, x_label, y_label, z_label,
type_var, lagdistance, lineartolerance, htol, vtol, hband, vband, azimuth, dip, nlags,
head_property, tail_property, choice, omni):
'''calculate_experimental_function
Args:
dataset (pd.DataFrame): Pandas dataframe containing all values
x_label (string): Label for x coordinates
y_label (string): Label for y coordinates
z_label (string): Label for z coordinates
type_var (string): Experimental continuity function
type_var:
'Variogram'
'Covariogram'
'Correlogram'
'PairWise'
'Relative_Variogram'
lagdistance (float): Size of lag size
lineartolerance (float): Linear tolerance
htol (float): Horizontal angular tolerance
vtol (float): Vertical angular tolerance
hband (float): Horizontal band width
vband (float): Vertical band width
azimuth (float): Azimuth direction in degrees
dip (float): Float direction in degrees
nlags (int): Number of lags to calculate experimental functions
head_property (string): String with name of Head Property
tail_property (string): String with anme of Tail Property
choice (float): Percentual of random sampling 0-1
omni (bool): Omnidirecional option True =Use omnidirecional , False = Don't use omnidirecional
Returns:
df (pandas.DataFrame): Pandas Dataframe containing the experimental continuity functions of all lags
'''
# Calculate all permissible distances
dist = _distances(dataset, nlags, x_label, y_label, z_label, lagdistance, head_property, tail_property, choice)
# Define the angles and tolerances
cos_Azimuth = np.cos(np.radians(90-azimuth))
sin_Azimuth = np.sin(np.radians(90-azimuth))
cos_Dip = np.cos(np.radians(90-dip))
sin_Dip = np.sin(np.radians(90-dip))
htol = np.abs(np.cos(np.radians(htol)))
vtol= np.abs(np.cos(np.radians(vtol)))
check_azimuth = np.abs((dist[:,0]*cos_Azimuth + dist[:,1]*sin_Azimuth)/dist[:,3])
check_dip = np.abs((dist[:,3]*sin_Dip + dist[:,2]*cos_Dip)/dist[:,4])
check_bandh = np.abs(cos_Azimuth*dist[:,1]- sin_Azimuth*dist[:,0])
check_bandv = np.abs(sin_Dip*dist[:,2] - cos_Dip*dist[:,3])
# Create experimental variogram dataset
number_of_int = range(1, (nlags +1))
value_exp = np.array([__calculate_experimental(dist, i, lagdistance, lineartolerance, check_azimuth,
check_dip, check_bandh, check_bandv, htol, vtol, hband, vband, omni,type_var) for i in number_of_int])
df = pd.DataFrame(value_exp,
columns = ['Spatial continuity', 'Number of pairs', 'Average distance'])
df = df.dropna()
return df
##############################################################################################################################################
def _calculate_experimental_function_varmap(type_var, lineartolerance, htol, vtol, hband, vband,
azimuth, dip,dataset, nlags,
X_rot, Y_rot,Z_rot,
lagdistance, head_property, tail_property, choice ,
omni = False, plot_graph=False, show_pairs=False):
'''calculate_experimental_function
Args:
type_var (string): Type of experimental variogram function
type_var:
'Variogram'
'Covariogram'
'Correlogram'
'PairWise'
'Relative_Variogram'
lineartolerance (TYPE): Description
htol (TYPE): Description
vtol (TYPE): Description
hband (TYPE): Description
vband (TYPE): Description
azimuth (TYPE): Description
dip (TYPE): Description
dataset (TYPE): Description
nlags (TYPE): Description
X_rot (TYPE): Description
Y_rot (TYPE): Description
Z_rot (TYPE): Description
lagdistance (TYPE): Description
head_property (TYPE): Description
tail_property (TYPE): Description
choice (TYPE): Description
omni (bool, optional): Description
plot_graph (bool, optional): Description
show_pairs (bool, optional): Description
Returns:
df (pandas.DataFrame): Pandas Dataframe containing the experimental continuity functions of all lags
'''
# Calculate all permissible distances
dist = _distances_varmap(dataset, X_rot, Y_rot,Z_rot, nlags, lagdistance, head_property, tail_property,choice)
# Define the angles and tolerances
cos_Azimuth = np.cos(np.radians(90-azimuth))
sin_Azimuth = np.sin(np.radians(90-azimuth))
cos_Dip = np.cos(np.radians(90-dip))
sin_Dip = np.sin(np.radians(90-dip))
htol = np.abs(np.cos(np.radians(htol)))
vtol= np.abs(np.cos(np.radians(vtol)))
check_azimuth = np.abs((dist[:,0]*cos_Azimuth + dist[:,1]*sin_Azimuth)/dist[:,3])
check_dip = np.abs((dist[:,3]*sin_Dip + dist[:,2]*cos_Dip)/dist[:,4])
check_bandh = np.abs(cos_Azimuth*dist[:,1]- sin_Azimuth*dist[:,0])
check_bandv = np.abs(sin_Dip*dist[:,2] - cos_Dip*dist[:,3])
# Create experimental variogram dataset
number_of_int = range(1, (nlags +1))
value_exp = np.array([__calculate_experimental(dist, i, lagdistance, lineartolerance, check_azimuth,
check_dip, check_bandh, check_bandv, htol, vtol, hband, vband, omni,type_var) for i in number_of_int])
df = pd.DataFrame(value_exp,
columns = ['Spatial continuity', 'Number of pairs', 'Average distance'])
df = df.dropna()
if plot_graph == True:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(df['Average distance'].values, df['Spatial continuity'].values)
ax.set_xlabel('Lag distance (h)')
ax.set_ylabel(type_c)
ax.set_title('Experimental continuity function : {} , azimuth {} and dip {} '.format(str(type_c), str(azimuth), str(dip)))
if show_pairs == True:
x, y = df['Average distance'].values, df['Spatial continuity'].values
for i, j in enumerate(df['Number of pairs'].values):
ax.annotate(str(j), xy =(x[i], y[i]), xytext =(x[i], (y[i]+0.05*y[i])))
ax.set_ylim((min(y),1.10*max(y)))
plt.grid()
plt.show()
return df
##############################################################################################################################################
def _modelling(experimental_dataframe,azimuths, dips, rotation_reference,
model_func, ranges, contribution, nugget, inverted= False, plot_graph = True ):
'''plot_experimental_function)_omni
Args:
experimental_dataframe (pd.DataFrame): pandas Data Frame containing the experimental continuity functions
azimuths (lst): List of azimuth directions
dips (lst): List of dip directions
rotation_reference (lst): List containing the azimuth and dips for all experimental variograms [azm, dip]
model_func (str): String containing the experimental function model
model_func:
'Exponential'
'Gaussian'
'Spherical'
ranges (lst): List containing the ranges for each structure for experimental continuity functionns
contribution (lst): List containing the variogram model contribution for each structure
nugget (float): The nugget effect
inverted (bool, optional): Invert variogram modelling to a covariogram modelling
plot_graph (bool, optional): Option to plot graph
Returns:
plot (matplotlib.pyplot): Plot of omnidirecional experimental continuity function and the spatial continuity model for one direction
Raises:
ValueError: Number of principal directions must be less or equal 3 or Variogram structures not have the same size
'''
# Give errors for modelling inconsistences
if len(ranges[0]) != 3:
raise ValueError("Number of principal directions must be less or equal 3")
if len(ranges) == len(contribution) and len(ranges) == len(model_func):
pass
else:
raise ValueError("Variogram structures must be have the same size")
# Calculate cartesian coordiantes according polar approach
y = math.cos(math.radians(dips))*math.cos(math.radians(azimuths))
x = math.cos(math.radians(dips))*math.sin(math.radians(azimuths))
z = math.sin(math.radians(dips))
# Define reference plane direction s
angle_azimuth = math.radians(rotation_reference[0])
angle_dip = math.radians(rotation_reference[1])
angle_rake = math.radians(rotation_reference[2])
# Define rotation matrix
rotation1 = np.array([[math.cos(angle_azimuth), -math.sin(angle_azimuth), 0],
[math.sin(angle_azimuth), math.cos(angle_azimuth), 0],
[0,0,1]])
rotation2 = np.array([[1, 0, 0],
[0, math.cos(angle_dip), math.sin(angle_dip)],
[0,-math.sin(angle_dip),math.cos(angle_dip)]])
rotation3 = np.array([[math.cos(angle_rake), 0, -math.sin(angle_rake)],
[0, 1, 0],
[math.sin(angle_rake),0,math.cos(angle_rake)]])
rot1 = np.dot(rotation1.T, np.array([x,y,z]))
rot2 = np.dot(rotation2.T, rot1)
rot3= np.dot(rotation3.T,rot2)
# Rotate samples
rotated_range =[]
for i in ranges:
rangex = float(i[1])
rangey = float(i[0])
rangez = float(i[2])
rotated = (np.multiply(rot3, np.array([rangex, rangey, rangez]).T))
rotated_range.append(math.sqrt(rotated[0]**2+rotated[1]**2+rotated[2]**2))
# Calculate model function
distancemax = experimental_dataframe['Average distance'].max()
distances = np.linspace(0, distancemax, 200)
model = []
if inverted == False:
for i in distances:
soma = 0
for j, o, l in zip(contribution, model_func, rotated_range):
if o == 'Exponential':
soma += j*(1-math.exp(-3*i/float(l)))
if o == 'Gaussian':
soma += j*(1-math.exp(-3*(i/float(l))**2))
if o == 'Spherical':
if i <= l:
soma += j*(1.5*i/float(l)-0.5*(i/float(l))**3)
else:
soma += j
soma += nugget
model.append(soma)
else:
for i in distances:
soma = 0
for j, o, l in zip(contribution, model_func, rotated_range):
if o == 'Exponential':
soma += (j+nugget)*(math.exp(-3*i/float(l)))
if o == 'Gaussian':
soma += (j+nugget)*(math.exp(-3*(i/float(l)**2)))
if o == 'Spherical':
if i <= l:
soma += (j+nugget)*(1 - (1.5*i/float(l)-0.5*(i/float(l))**3))
else:
soma += 0
model.append(soma)
df = pd.DataFrame(np.array([distances, model]).T, columns= ['distances', 'model'])
return df
##############################################################################################################################################
def _varmap(azimuth, dip, lineartolerance, htol, vtol, hband, vband, type_var, dataset, x_label, y_label, z_label, nlags,
lagdistance, head_property, tail_property, ndirections, wait_please, choice = 1.00, interactive = False):
"""Summary
Args:
azimuth (float): Azimuth of reference plane
dip (dip): Dip of the reference plane
lineartolerance (TYPE): Linear tolerance
htol (float): Horizontal tolerance
vtol (float): Vertical tolerance
hband (float): Horizontal bandwidth
vband (float): Vertical bandwidth
type_var (float): Experimental continuity function
type_var:
'Variogram'
'Covariogram'
'Correlogram'
'PairWise'
'Relative_Variogram'
dataset (pd.DataFrame): pandas DataFrame containing the data
x_label (str): Label for x coordinates
y_label (str): Label for y coordinates
z_label (str): Label for z coordinates
nlags (int): Number of lags to calculate experimental variogram
lagdistance (float): Lag size
head_property (str): String for the head property
tail_property (str): String for the tail property
ndirections (int): Number of directions for calculate experimental function
wait_please (obj): Object to update progress bar
choice (float, optional): Values to random select experimental variograms
interactive (bool, optional): Enable interactive varmaps
"""
# Define X, Y, Z coordinates
X = dataset[x_label].values
Y = dataset[y_label].values
Z = dataset[z_label].values
# Set progress bar to zero
if interactive == True:
wait_please.value = 0
# rotate dataset acording the reference plane
Xrot, Yrot, Zrot = _rotate_data(X, Y, Z, azimuth, dip)
# Calculate experimental functions
lag_adm = []
azimute_adm = []
continuidade = []
for i in np.arange(0,360,ndirections):
if interactive == True:
wait_please.value += 1
df_exp = _calculate_experimental_function_varmap(type_var, lineartolerance, htol, vtol, hband, vband,
i, 0,dataset, nlags, Xrot, Yrot,Zrot,lagdistance, head_property, tail_property, choice)
var_val = df_exp['Spatial continuity'].values
continuidade = np.append(continuidade,var_val)
values_1 = df_exp['Average distance'].values
lag_adm = np.append(lag_adm,values_1)
values_2 = np.array([np.radians(i) for j in range(df_exp.shape[0])])
azimute_adm = np.append(azimute_adm, values_2)
# Define the interplation parameters
gdiscrete = 40
ncontour = 40
# define the x and y coordinates in cartesian way
x = np.array(lag_adm)*np.sin(azimute_adm)
y = np.array(lag_adm)*np.cos(azimute_adm)
# Define the maximum grid size
maximo = max([max(x), max(y)])
# Filter values less than a circle distance
ray = maximo
truzao = []
for j, i in enumerate(x):
if math.sqrt(x[j]**2 + y[j]**2) > ray:
truzao.append(False)
else:
truzao.append(True)
x= x[truzao]
y =y[truzao]
continuidade = continuidade[truzao]
x= x.ravel()
y= y.ravel()
# Create grid and interpolation
Xi = np.linspace(-maximo,maximo,gdiscrete)
Yi = np.linspace(-maximo,maximo,gdiscrete)
f = plt.figure(figsize=(15,10))
left, bottom, width, height= [0,0.1, 0.7, 0.7]
ax = plt.axes([left, bottom, width, height])
pax = plt.axes([left, bottom, width, height],
projection='polar',
facecolor='none')
pax.set_theta_zero_location("N")
pax.set_theta_direction(-1)
pax.set_xticks(np.pi/180. * np.linspace(0, 360, 36, endpoint=False))
ax.set_aspect(1)
ax.axis('Off')
# grid the data.
Vi = griddata((x, y), np.array(continuidade), (Xi[None,:], Yi[:,None]), method='linear')
cf = ax.contourf(Xi,Yi,Vi, ncontour, cmap=plt.cm.jet)
# Define gradient
gradient = np.linspace(1,0, 256)
gradient = np.vstack((gradient, gradient))
# Create matplotlib graph
cax = plt.axes([0.72,0.1, 0.05, 0.7])
cax.xaxis.set_major_locator(plt.NullLocator())
cax.yaxis.tick_right()
cax.imshow(gradient.T, aspect='auto', cmap=plt.cm.jet)
cax.set_yticks(np.linspace(0,256,12))
cax.set_yticklabels(list(map(str, cf.get_array()))[::-1])
plt.show()
##############################################################################################################################################
# EM FASE DE TESTE
def _covariogram_map_3d(property_value, dataset, x_label, y_label, z_label, neighbors, division = 20, alpha= 0.7,
cutx =[-np.inf, np.inf],cuty =[-np.inf,np.inf],cutz =[-np.inf,np.inf], size =20 ):
'''covariogram_map_3d
Args:
property_value (TYPE): Description
dataset (TYPE): Description
x_label (TYPE): Description
y_label (TYPE): Description
z_label (TYPE): Description
neighbors (int): Number of neighbors using in KNearest neighbors
division (int, optional): Description
alpha (float, optional): Description
cutx (list, optional): list containing the minimum cutsize and the maximum cutsize for x coordinates
cuty (list, optional): list containing the minimum cutsize and the maximum cutsize for y coordinates
cutz (list, optional): list containing the minimum cutsize and the maximum cutsize for z coordinates
size (int, optional): Description
'''
X = dataset[x_label].values
Y = dataset[y_label].values
Z = dataset[z_label].values
R = dataset[property_value].values
max_x, min_x = max(X), min(X)
max_y, min_y = max(Y), min(Y)
max_z, min_z = max(Z), min(Z)
cordinatesx = np.linspace(min_x,max_x,division)
cordinatesy = np.linspace(min_y,max_y, division)
cordinatesz = np.linspace(min_z,max_z,division)
cordinates = np.array([np.array([i,j,k]).T for i in cordinatesx for j in cordinatesy for k in cordinatesz])
estimates = np.zeros(len(cordinates))
nb = KNeighborsRegressor(n_neighbors=neighbors).fit(np.array([X,Y,Z]).T, R)
for i, j in zip(range(len(estimates)), cordinates):
estimates[i] = nb.predict(j.reshape(1, -1))[0]
estimates = estimates.reshape((division,division,division))
fft_u = np.conjugate(np.fft.fftn(estimates,axes=(0,1,2), norm ='ortho'))
fft_u_roll = np.fft.fftn(estimates,axes=(0,1,2), norm ='ortho')
product = np.multiply(fft_u_roll,fft_u)
Covariance = np.fft.fftshift(np.fft.ifftn(product,axes=(0,1,2), norm ='ortho'))
Covariance = np.real(Covariance)
Covariance = Covariance.reshape(-1,1)[:,0]/(division*division*division)
filter_cord = (cordinates[:,0]>cutx[0]) & (cordinates[:,0]<cutx[1]) & (cordinates[:,1]>cuty[0]) & (cordinates[:,1]<cuty[1]) & (cordinates[:,2]>cutz[0]) & (cordinates[:,2]<cutz[1])
cx = cordinates[:,0][filter_cord]
cy = cordinates[:,1][filter_cord]
cz = cordinates[:,2][filter_cord]
Covariance = Covariance[filter_cord]
df_varmap = pd.DataFrame(np.array([cx, cy, cz, Covariance]).T, columns=['cx','cy','cz','Covariance'])
fig = px.scatter_3d(df_varmap, x='cx', y='cy', z='cz',
color='Covariance')
fig.show()
##############################################################################################################################################
def _modelling_to_interact(**kargs):
"""Summary
Args:
**kargs: Args to modelling variogram interactively
Returns:
TYPE: Variogram model
"""
# Extract model parameters
rotation_reference = [kargs.get('rotation_max') ,kargs.get('rotation_med'),kargs.get('rotation_min') ]
ranges =[]
for i in range(kargs.get('nstructures')):
mx_rg = 'rangemax_{}'.format(str(i))
md_rg = 'rangemed_{}'.format(str(i))
mi_rg = 'rangemin_{}'.format(str(i))
ranges.append([kargs.get(mx_rg), kargs.get(md_rg), kargs.get(mi_rg)])
model_func = []
contribution = []
for i in range(kargs.get('nstructures')):
mld = 'model_{}'.format(str(i))
contr = 'contribution_{}'.format(str(i))
model_func.append(kargs.get(mld))
contribution.append(kargs.get(contr))
nugget = kargs.get('nugget')
inverted = kargs.get('inverted')
azimuths = kargs.get('azimuths')
dips = kargs.get('dips')
# Calculate models for each experimental directions
dfs = []
for j, i in enumerate(kargs.get('experimental_dataframe')):
dfs.append(_modelling(i, azimuths[j], dips[j], rotation_reference,model_func ,ranges,contribution,nugget, inverted))
# Plot graph
size_row = 1 if len(dfs) < 4 else int(math.ceil(len(dfs)/4))
size_cols = 4 if len(dfs) >= 4 else int(len(dfs))
titles = ["Azimuth {} - Dip {}".format(azimuths[j], dips[j]) for j in range(len(dfs))]
fig = make_subplots(rows=size_row, cols=size_cols, subplot_titles=titles)
count_row = 1
count_cols = 1
for j, i in enumerate(dfs):
fig.add_trace(go.Scatter(x=kargs.get('experimental_dataframe')[j]['Average distance'], y=kargs.get('experimental_dataframe')[j]['Spatial continuity'],
mode='markers',
name='Experimental',
marker= dict(color =kargs.get('experimental_dataframe')[j]['Number of pairs'] ),
text =kargs.get('experimental_dataframe')[j]['Number of pairs'].values if kargs.get('show_pairs')== True else None) , row=count_row, col=count_cols)
fig.add_trace(go.Scatter(x=i['distances'], y=i['model'],
mode='lines',
name='Model'), row=count_row, col=count_cols )
fig.update_xaxes(title_text="Distance", row=count_row, col=count_cols, automargin = True)
fig.update_yaxes(title_text="Variogram", row=count_row, col=count_cols, automargin=True)
fig.update_layout(autosize=True)
count_cols += 1
if count_cols > 4:
count_cols = 1
count_row += 1
fig.show()
# Store variogram model
global return_model_var
return_model_var = {'number_of_structures' : kargs.get('nstructures'),
'rotation_reference': rotation_reference,
'models': model_func,
'nugget': nugget,
'contribution': contribution,
'ranges' : ranges}
return return_model_var
##############################################################################################################################################
[docs]def interactive_modelling(experimental_dataframe: callable, number_of_structures: int , show_pairs: bool = False):
"""Opens the interactive modeling controlls. Store the results in a global variable `gammapy.return_model_var`.
Args:
experimental_dataframe (DataFrame): Experimental variogram DataFrame. Assessed by `gammapy.return_exp_var`
number_of_structures (int): number of structures
show_pairs (bool, optional): show number of pairs flag. Defaults to False.
Returns:
dict: Variogram model
"""
warnings.filterwarnings('ignore')
azimuths = experimental_dataframe['Directions'][0]
dips = experimental_dataframe['Directions'][1]
rotation_max = widgets.FloatSlider(description='Rotação Azim', min= 0, max= 360, step=1,continuous_update=False)
rotation_med = widgets.FloatSlider(description='Rotação Dip',min= 0, max= 360, step=1,continuous_update=False)
rotation_min = widgets.FloatSlider(description='Rotação Rake',min= 0, max= 360, step=1,continuous_update=False)
nugget = widgets.FloatSlider(description='Pepita',min= 0, max= 2*max(experimental_dataframe['Values'][0]['Spatial continuity']), step=10*max(experimental_dataframe['Values'][0]['Spatial continuity'])/1000,continuous_update=False)
inverted = widgets.Checkbox(description='Inverter a modelagem')
u1 = widgets.HBox([rotation_max, rotation_med, rotation_min])
grid1 = GridspecLayout(number_of_structures, 1)
values1 = []
for i in range(number_of_structures):
values1.append([widgets.FloatSlider(description='Alcance máximo', min= 0, max= 2*max(experimental_dataframe['Values'][0]['Average distance']), step=max(experimental_dataframe['Values'][0]['Average distance'])/1000.0,continuous_update=False),
widgets.FloatSlider(description='Alcance médio', min= 0, max= 2*max(experimental_dataframe['Values'][0]['Average distance']), step=max(experimental_dataframe['Values'][0]['Average distance'])/1000.0,continuous_update=False),
widgets.FloatSlider(description='Alcance mÃnimo',min= 0, max= 2*max(experimental_dataframe['Values'][0]['Average distance']), step=max(experimental_dataframe['Values'][0]['Average distance'])/1000.0,continuous_update=False)])
grid1[i,0] = widgets.HBox(values1[i])
u4 = widgets.HBox([nugget, inverted])
values2 = []
grid2 = GridspecLayout((number_of_structures+1), 1)
for i in range(number_of_structures):
values2.append([widgets.Dropdown(options= ['Spherical','Exponential','Gaussian']),
widgets.FloatSlider(description='Contribution',min= 0, max= 2*max(experimental_dataframe['Values'][0]['Spatial continuity']), step=max(experimental_dataframe['Values'][0]['Spatial continuity'])/1000,continuous_update=False)])
grid2[i,0] = widgets.HBox(values2[i])
grid2[number_of_structures,0] = u4
children = [u1,grid1,grid2]
accordion = widgets.Accordion(children=children)
accordion.set_title(0, 'Rotações dos eixos principais')
accordion.set_title(1, 'Alcance dos eixos principais ')
accordion.set_title(2, 'Função, contribuição e efeito pepita')
outputs = {'experimental_dataframe': fixed(experimental_dataframe['Values']),
'azimuths' : fixed(azimuths),
'dips' : fixed(dips),
'nstructures': fixed(number_of_structures),
'rotation_max' : rotation_max,
'rotation_med' :rotation_med,
'rotation_min' : rotation_min,
'nugget': nugget,
'inverted' : inverted}
for i in range(number_of_structures):
outputs['rangemax_{}'.format(str(i))] = values1[i][0]
outputs['rangemed_{}'.format(str(i))] = values1[i][1]
outputs['rangemin_{}'.format(str(i))] = values1[i][2]
outputs['model_{}'.format(str(i))] = values2[i][0]
outputs['contribution_{}'.format(str(i))] = values2[i][1]
outputs['show_pairs'] = fixed(show_pairs)
out = widgets.interactive_output(_modelling_to_interact, outputs)
display(accordion, out)
##############################################################################################################################################
##############################################################################################################################################
[docs]def modelling(experimental_dict: callable, rotation_max: float, rotation_med: float, rotation_min: float,
nugget: int ,inverted: bool, rangemax: list, rangemed: list, rangemin: list, model: list, contribution:list,
show_pairs: bool = False):
"""Variogram modeling function without interactive controlls.
Args:
experimental_dict (dict): experimental variogram
rotation_max (float): azimuth
rotation_med (float): dip
rotation_min (float): rake
nugget (float): nugget contribution
inverted (bool): inverted model flag
rangemax (lst): list of ranges in maximum continuity direction for each structure
rangemed (lst): list of ranges in intermediate continuity direction for each structure
rangemin (lst): list of ranges in minumim continuity direction for each structure
model (lst of str): lst of models for each structure. `Spherical`, `Exponential` or `Gaussian`.
contribution (lst of floats): list of contributions for each structure
show_pairs (bool, optional): Show pairs flag. Defaults to False.
Returns:
dict: variogram model DataFrame
"""
warnings.filterwarnings('ignore')
azimuths = experimental_dict['Directions'][0]
dips = experimental_dict['Directions'][1]
number_of_structures = len(rangemax)
outputs = {'experimental_dataframe': experimental_dict['Values'],
'azimuths' : azimuths,
'dips' : dips,
'nstructures': number_of_structures,
'rotation_max' : rotation_max,
'rotation_med' :rotation_med,
'rotation_min' : rotation_min,
'nugget': nugget,
'inverted' : inverted}
for i in range(number_of_structures):
outputs['rangemax_{}'.format(str(i))] = rangemax[i]
outputs['rangemed_{}'.format(str(i))] = rangemed[i]
outputs['rangemin_{}'.format(str(i))] = rangemin[i]
outputs['model_{}'.format(str(i))] = model[i]
outputs['contribution_{}'.format(str(i))] = contribution[i]
outputs['show_pairs'] = show_pairs
return _modelling_to_interact(**outputs)
##############################################################################################################################################
[docs]def interactive_varmap(dataset:callable, X:str, Y:str, head:str, tail:str, Z:str =None, choice:float =1.0):
"""Opens interactive variogram map controls.
Args:
dataset (DataFrame): Data points DataFrame
X (str): x column coordinates name
Y (str): y column cordinates name
head (str): head property name
tail (str): tail property name
Z (str, optional): z coordinates column name. Defaults to None.
choice (float, optional): pool a random number of data to calculate the variogram. Defaults to 1.0.
"""
# set same seed
# Define the Ipython Widgets
warnings.filterwarnings('ignore')
azimuth = widgets.BoundedFloatText(value=0,min = 0, max = 360, description='Azimuth:',disabled=False)
dip = widgets.BoundedFloatText(value=0,min= 0 , max=360, description='Dip:',disabled=False)
nlags = widgets.BoundedIntText(value=5,min= 1, max= 1000000, description='Nlags:',disabled=False)
lagdistance = widgets.BoundedFloatText(value=7.5, min= 0.0000001, max=1000000, description='Lag size:',disabled=False)
ndirections = widgets.Dropdown(options=[9, 18, 36],value=18,description='Ndirections:',disabled=False)
type_var = widgets.Dropdown(options=['Variogram', 'Covariogram', 'Correlogram', 'PairWise', 'Relative_Variogram'],
value='Variogram',description='Number:',disabled=False)
execute = widgets.Button(description='Execute',icon='check')
output = widgets.Output()
# Defining progress bar
max_count = 1000
wait_please = widgets.IntProgress(min=0, max=ndirections.value) # instantiate the bar
# Define tolerance properties according Ipython values
lineartolerance = lagdistance.value/2.0
hband = lagdistance.value/2.0
vband= lagdistance.value/2.0
htol = 360/(ndirections.value)
vtol = 360/(ndirections.value)
# Creating Hbox layout
u1 = widgets.HBox([azimuth, dip])
u2 = widgets.HBox([nlags, lagdistance])
u3 = widgets.HBox([type_var, ndirections])
u4 = widgets.HBox([execute, wait_please])
# Creating grid layout
grid = GridspecLayout(4, 1)
grid[0,0] = u1
grid[1,0] = u2
grid[2,0] = u3
# Creating accordion layout
children = [grid]
accordion = widgets.Accordion(children=children)
accordion.set_title(0, 'Varmap parameters')
display(accordion, u4) # display values
# Calling varmap according 2D or 3D dataset
if Z == None:
dataset['Z'] = np.zeros(dataset[X].values.shape[0])
def on_varmap(change):
"""Summary
Args:
change (TYPE): Description
"""
_varmap(azimuth.value,0, lineartolerance, htol, vtol, hband, vband, type_var.value, dataset, X, Y, 'Z', nlags.value,
lagdistance.value, head, tail, ndirections.value, wait_please, choice, True)
else:
def on_varmap(change):
"""Summary
Args:
change (TYPE): Description
"""
_varmap(azimuth.value, dip.value, lineartolerance, htol, vtol, hband, vband, type_var.value, dataset, X, Y, Z, nlags.value,
lagdistance.value, head, tail, ndirections.value, wait_please, choice, True)
execute.on_click(on_varmap)
##############################################################################################################################################
[docs]def varmap(dataset: callable, X: str, Y: str, head: str, tail: str, azimuth: float,
dip: float, nlags: int, lagdistance: float, ndirections: int, type_var: str, Z:str =None, choice =1.0):
"""Plots a variogram map.
Args:
dataset (DataFrame): Point set DataFrame
X (str): x coordinates column name
Y (str): y coordinates column name
head (str): head variable name
tail (str): tail variable name
azimuth (float): azimuth
dip (float): dip
nlags (int): number of lags
lagdistance (float): lag distance
ndirections (int): number of directions
type_var (str): covariance funcrion type. `Variogram`, `Correlogram` ...
Z (str, optional): z coordinates variable name. Defaults to None.
choice (float, optional): pool a random number of data to calculate the variogram. Defaults to 1.0.
"""
# set same seed
# Define tolerance properties according Ipython values
warnings.filterwarnings('ignore')
lineartolerance = lagdistance/2.0
hband = lagdistance/2.0
vband= lagdistance/2.0
htol = 360/(ndirections)
vtol = 360/(ndirections)
wait_please = 0
# Calling varmap according 2D or 3D dataset
if Z == None:
dataset['Z'] = np.zeros(dataset[X].values.shape[0])
_varmap(azimuth,0, lineartolerance, htol, vtol, hband, vband, type_var, dataset, X, Y, 'Z', nlags,
lagdistance, head, tail, ndirections, wait_please, choice)
else:
_varmap(azimuth, dip, lineartolerance, htol, vtol, hband, vband, type_var, dataset, X, Y, Z, nlags,
lagdistance, head, tail, ndirections, wait_please, choice)
##############################################################################################################################################
[docs]def interactive_experimental(dataset:callable, X:str, Y:str, head:str, tail:str, ndirections:int, show_pairs =False, Z =None, choice =1.0):
"""Calculates experimental variogram. Store the results in a global variable `gammapy.return_exp_var`.
Args:
dataset (DataFrame): data points DataFrame
X (str): x coordinates column name
Y (str): y coordinates column name
head (str): head variable
tail (str): tail variable
ndirections (int): number of directions
show_pairs (bool, optional): show number of pairs flag. Defaults to False.
Z (str, optional): z coordinates column name. Defaults to None.
choice (float, optional): pool a random number of data to calculate the variogram. Defaults to 1.0.
Returns:
dict: experimental variogram DataFrame
"""
warnings.filterwarnings('ignore') # Ignore numpy erros
# If dataset 2D, fill Z values with zeros
if Z == None:
dataset['Z'] = np.zeros(dataset[X].values.shape[0])
Z = 'Z'
# Creating Ipython Widgets
widgets_l = []
vboxes = []
for i in range(ndirections):
widgets_l.append([widgets.Dropdown(options=['Variogram', 'Covariogram', 'Correlogram','PairWise', 'RelativeVariogram' ],value='Variogram',description='Type:',disabled=False),
widgets.BoundedFloatText(value=0,min = 0, max = 360, description='Azimuth:',disabled=False),
widgets.BoundedFloatText(value=0,min = 0, max = 360, description='Dip:',disabled=False),
widgets.BoundedIntText(value=10,min = 1, max = 10000, description='nlags:',disabled=False),
widgets.BoundedFloatText(value=10,min = 0.000000001, max = 1000000000, description='lagsize:',disabled=False),
widgets.BoundedFloatText(value=10,min = 0.000000001, max = 1000000000, description='lineartol:',disabled=False),
widgets.BoundedFloatText(value=10,min = 0.000000001, max = 90, description='htol:', disabled=False),
widgets.BoundedFloatText(value=10,min = 0.000000001, max = 90, description='vtol:', disabled=False),
widgets.BoundedFloatText(value=10,min = 0.000000001, max = 90, description='hband:', disabled=False),
widgets.BoundedFloatText(value=10,min = 0.000000001, max = 1000000000, description='vband:', disabled=False),
widgets.Dropdown(options=[True, False ],value=False,description='Omni:',disabled=False)])
vboxes.append(widgets.VBox(widgets_l[i], layout=widgets.Layout(width='370px')))
execute = widgets.Button(description='Execute',icon='check')
output = widgets.Output()
# Defining progress bar
max_count = 1000
wait_please = widgets.IntProgress(min=0, max=ndirections) # instantiate the bar
# More boxes and layouts
grid = widgets.HBox(vboxes, layout=widgets.Layout(width='5000px'))
box2 = widgets.HBox([execute, wait_please], display='flex',align_items='stretch',width='50%')
box3 = widgets.VBox([grid, box2], display='flex',align_items='stretch',width='50%')
children = [box3]
accordion = widgets.Accordion(children=children)
accordion.set_title(0, 'Experimental parameters')
display(accordion)
# Define variograms with change
returning = []
def on_variogram(change):
# fill parameters
type_c, nlags, lagsize, lineartol,htol, vtol, hband, vband, azimuth, dip, omni = [], [], [], [], [], [] , [] ,[] ,[] ,[], []
for i in range(ndirections):
type_c.append(widgets_l[i][0].value)
azimuth.append(widgets_l[i][1].value)
dip.append(widgets_l[i][2].value)
nlags.append(widgets_l[i][3].value)
lagsize.append(widgets_l[i][4].value)
lineartol.append(widgets_l[i][5].value)
htol.append(widgets_l[i][6].value)
vtol.append(widgets_l[i][7].value)
hband.append(widgets_l[i][8].value)
vband.append(widgets_l[i][9].value)
omni.append(widgets_l[i][10].value)
# Calculate experimental variograms
dfs = []
for i in range(ndirections):
wait_please.value += 1
dfs.append(_calculate_experimental_function(dataset, X, Y, Z,
type_c[i], lagsize[i], lineartol[i], htol[i], vtol[i], hband[i], vband[i], azimuth[i], dip[i],
nlags[i], head, tail, choice, omni[i]))
# Store experimental variograms
global return_exp_var
returning = {'Directions': [azimuth, dip],
'Values' : dfs}
return_exp_var = returning
# Plot variograms
size_row = 1 if len(dfs) < 4 else int(math.ceil(len(dfs)/4))
size_cols = 4 if len(dfs) >= 4 else int(len(dfs))
titles = ["Azimuth {} - Dip {}".format(azimuth[j], dip[j]) if omni[j]==False else 'Omni' for j in range(len(dfs))]
fig = make_subplots(rows=size_row, cols=size_cols, subplot_titles=titles)
count_row = 1
count_cols = 1
for j, i in enumerate(dfs):
fig.add_trace(go.Scatter(x=dfs[j]['Average distance'], y=dfs[j]['Spatial continuity'],
mode='markers',
name='Experimental' ,
marker= dict(color =dfs[j]['Number of pairs']),
text =dfs[j]['Number of pairs'].values if show_pairs== True else None,
textposition='bottom center') , row=count_row, col=count_cols)
fig.update_xaxes(title_text="Distance", row=count_row, col=count_cols, automargin = True)
fig.update_yaxes(title_text="Variogram", row=count_row, col=count_cols, automargin=True)
fig.update_layout(autosize=True)
count_cols += 1
if count_cols > 4:
count_cols = 1
count_row += 1
fig.show()
execute.on_click(on_variogram)
##############################################################################################################################################
##############################################################################################################################################
[docs]def experimental(dataset:callable, X:str, Y:str, head:str, tail:str,type_c:str,
nlags:list , lagsize:list, lineartol:list ,htol: list, vtol: list, hband: list, vband:list, azimuth:list, dip:list, omni:bool,
show_pairs:bool =False, Z:str =None, choice:float =1.0):
"""Calculates experimental variogram. Store the results in a global variable `gammapy.return_exp_var`.
Args:
dataset (DataFrame): point set DataFrame
X (str): x coordinates column name
Y (str): y coordinates column name
head (str): head variable name
tail (str): tail variable name
type_c (lst of str): covariance funcrion type list. `Variogram`, `Correlogram` ...
nlags (lst): list of lag number for each direction
lagsize (lst): list of lag size for each direction
lineartol (lst): list of linear tolerance for each direction
htol (lst): list of horizontal tolerance for each direction
vtol (lst): liist of vertical tolerance for each direction
hband (lst): list of horizontal band for each direction
vband (lst): list of vertical band for each direction
azimuth (lst): azimuth
dip (lst): dip
omni (bool): omnidirectionl flag
show_pairs (bool, optional): show pairs flag. Defaults to False.
Z (str, optional): z coordinates column name. Defaults to None.
choice (float, optional): pool a random number of data to calculate the variogram. Defaults to 1.0.
Returns:
dict: experimental variogram DataFrame
"""
# Ignore numpy erros
warnings.filterwarnings('ignore')
# define the number of directions in experimental variogram
ndirections = len(nlags)
# Fill dataset with 0 z coordinates if 2D case
if Z == None:
dataset['Z'] = np.zeros(dataset[X].values.shape[0])
Z = 'Z'
# Return the dictionary containing all information of experimental variograms
dfs = []
returning = {'Directions': [azimuth, dip],
'Values' : dfs}
# Calculate experimental variograms for each directions
for i in range(ndirections):
dfs.append(_calculate_experimental_function(dataset, X, Y, Z,
type_c[i], lagsize[i], lineartol[i], htol[i], vtol[i], hband[i], vband[i], azimuth[i], dip[i],
nlags[i], head, tail, choice, omni[i]))
# Plot experimental variograms
size_row = 1 if len(dfs) < 4 else int(math.ceil(len(dfs)/4))
size_cols = 4 if len(dfs) >= 4 else int(len(dfs))
titles = ["Azimuth {} - Dip {}".format(azimuth[j], dip[j]) if omni[j]==False else 'Omni' for j in range(len(dfs))]
fig = make_subplots(rows=size_row, cols=size_cols, subplot_titles=titles)
count_row = 1
count_cols = 1
for j, i in enumerate(dfs):
fig.add_trace(go.Scatter(x=dfs[j]['Average distance'], y=dfs[j]['Spatial continuity'],
mode='markers',
name='Experimental' ,
marker= dict(color =dfs[j]['Number of pairs']),
text =dfs[j]['Number of pairs'].values if show_pairs== True else None,
textposition='bottom center') , row=count_row, col=count_cols)
fig.update_xaxes(title_text="Distance", row=count_row, col=count_cols, automargin = True)
fig.update_yaxes(title_text="Variogram", row=count_row, col=count_cols, automargin=True)
fig.update_layout(autosize=True)
count_cols += 1
if count_cols > 4:
count_cols = 1
count_row += 1
fig.show()
return returning
##############################################################################################################################################
[docs]def hscatterplots (dataset, X:str, Y:str, head:str, tail:str, lagsize: float, nlags: int,
azimuth: float, dip: float, lineartol: float,htol: float, vtol: float, hband: float, vband: float , Z:str= None, choice:float=1.0, figsize=(700,700)):
"""Plots H-scatterplot for n lags
Args:
dataset (pd.DataFrame): Pandas DataFrame
X (str): x coordinates column name
Y (str): y coordinates column name
head (str): head variable name
tail (str): tail variable name
lagsize (float): lag size
nlags (int): number of lags
azimuth (float): azimuth
dip (float): dip
lineartol (float): linear tolerance
htol (float): angular horizontal tolerance in degrees
vtol (float): angular vertical tolerance in degrees
hband (float): horinztal bandwidth
vband (float): vertical bandwidth
Z (str, optional): z coordinates column name. Defaults to None.
choice (float, optional): pool a random number of data to calculate the variogram. Defaults to 1.0. """
# If dataset 2D, fill Z values with zeros
warnings.filterwarnings('ignore')
if Z == None:
dataset['Z'] = np.zeros(dataset[X].values.shape[0])
Z = 'Z'
# Calculate all permissible distances
dist = _distances(dataset, nlags, X, Y, Z, lagsize, head, tail, choice)
# Define the angles and tolerances
cos_Azimuth = np.cos(np.radians(90-azimuth))
sin_Azimuth = np.sin(np.radians(90-azimuth))
cos_Dip = np.cos(np.radians(90-dip))
sin_Dip = np.sin(np.radians(90-dip))
htol = np.abs(np.cos(np.radians(htol)))
vtol= np.abs(np.cos(np.radians(vtol)))
check_azimuth = np.abs((dist[:,0]*cos_Azimuth + dist[:,1]*sin_Azimuth)/dist[:,3])
check_dip = np.abs((dist[:,3]*sin_Dip + dist[:,2]*cos_Dip)/dist[:,4])
check_bandh = np.abs(cos_Azimuth*dist[:,1]- sin_Azimuth*dist[:,0])
check_bandv = np.abs(sin_Dip*dist[:,2] - cos_Dip*dist[:,3])
store = []
lagmultiply = []
for i in range(1,(nlags+1)):
lagmultiply.append(i)
points = __permissible_pairs (i, lagsize, lineartol, check_azimuth,
check_dip, check_bandh, check_bandv, htol, vtol, hband, vband, False,
dist)
average_tail = (points[:,7] + points[:,8])/2
average_head = (points[:,5] + points[:,6])/2
store.append([average_head, average_tail])
plots.plot_hscat(store, lagsize, lagmultiply, figsize)