<jupyter_start><jupyter_text>Variogram map test This function is under development<jupyter_code>%matplotlib inline
import pygslib  

pygslib.version.__version__
# we read the file header
# use >>> print. pygslib.gslib.read_header.__doc__ for help

#the data file is in the working directory
mydata= pygslib.gslib.read_gslib_file('../datasets/cluster.dat') 

#adding elevation and a dummy BHID
mydata['Zlocation']=0
mydata['bhid']=1


print (' \n **** 5 first rows in my datafile \n\n  ', mydata.head(n=5))<jupyter_output>**** 5 first rows in my datafile 

      Xlocation  Ylocation  Primary  Secondary  Declustering Weight  Zlocation  \
0       39.5       18.5     0.06       0.22                1.619          0   
1        5.5        1.5     0.06       0.27                1.619          0   
2       38.5        5.5     0.08       0.40                1.416          0   
3       20.5        1.5     0.09       0.39                1.821          0   
4       27.5       14.5     0.09       0.24                1.349          0   

   bhid  
0     1  
1     1  
2     1  
3     1  
4     1<jupyter_text>Defining and plotting direction over a plane<jupyter_code>#define a direction:
import numpy as np
import matplotlib.pyplot as plt

def hdir(r=1, ndir=18, refdir=0):
    theta = np.linspace(0, np.pi, ndir+1)[:-1]
    ax = plt.subplot(111, polar=True)
    
    # make clockwise like maps
    ax.set_theta_direction(-1)
    # make the plot to point north 
    ax.set_theta_zero_location("N")

    for t in theta: 
        pass
        ax.plot([0,t], [0,r], color='r', linewidth=1)
        ax.plot([0,t+np.pi], [0,r], color='r', linewidth=1)

    ax.grid(True)

    
    ax.invert_yaxis()
    ax.set_title("A line plot on a polar axis", va='bottom')
    
    
    ax.plot([0,refdir], [0,r], color='b', linewidth=3)
    ax.plot([0,refdir+np.pi], [0,r], color='b', linewidth=3)

    plt.show()
    
    
    return np.rad2deg(theta) #the last direction 
    
azm= hdir(r=1, ndir=18, refdir=0)
plt.show()
print (len (azm), azm)
ndir=18
azm= hdir(r=1, ndir=ndir, refdir=0)  #this produces a plot + directions
atol= np.ones(ndir)*180/ndir/2
dip = np.zeros(ndir)
bandwh = np.ones(ndir)*5000
bandwd = np.ones(ndir)*5000
dtol = np.ones(ndir)*5000
sills = [np.var(mydata['Primary'])]
ivtail = np.ones(ndir)
ivhead = np.ones(ndir)
ivtype = np.ones(ndir)*7


parameters = { 
'x'      :  mydata['Xlocation']   , # X coordinates, array('f') with bounds (nd), nd is number of data points
'y'      :  mydata['Ylocation'],    # Y coordinates, array('f') with bounds (nd)
'z'      :  mydata['Zlocation'],    # Z coordinates, array('f') with bounds (nd)
'bhid'   :  mydata['bhid'],         # bhid for downhole variogram, array('i') with bounds (nd)    
'vr'     :  mydata['Primary'],      # Variables, array('f') with bounds (nd,nv), nv is number of variables
'tmin'   : -1.0e21,                 # trimming limits, float
'tmax'   :  1.0e21,                 # trimming limits, float
'nlag'   :  8,                     # number of lags, int
'xlag'   :  6,                      # lag separation distance, float                
'xltol'  :  3,                      # lag tolerance, float
'azm'    : azm,                     # azimuth, array('f') with bounds (ndir)
'atol'   : atol,                    # azimuth tolerance, array('f') with bounds (ndir)
'bandwh' : bandwh,                  # bandwith h, array('f') with bounds (ndir)
'dip'    : dip,                     # dip, array('f') with bounds (ndir)
'dtol'   : dtol,              # dip tolerance, array('f') with bounds (ndir)
'bandwd' : bandwd,              # bandwith dip, array('f') with bounds (ndir)
'isill'  : 0,                       # standardize sills? (0=no, 1=yes), int
'sills'  : sills,                   # variance used to std the sills, array('f') with bounds (nv)
'ivtail' : ivtail,         # tail var., array('i') with bounds (nvarg), nvarg is number of variograms
'ivhead' : ivhead,         # head var., array('i') with bounds (nvarg)
'ivtype' : ivtype,         # variogram type, array('i') with bounds (nvarg)
'maxclp' : 50000}                   # maximum number of variogram point cloud to use, input int
#Now we are ready to calculate the variogram
pdis,pgam, phm,ptm,phv,ptv,pnump, cldi, cldj, cldg, cldh = pygslib.gslib.gamv(parameters)
#get parameters
print (type(pdis), pdis.shape)

# for example
nvrg = pdis.shape[0]
ndir = pdis.shape[1]
nlag = pdis.shape[2]-2

%matplotlib inline
import matplotlib.pyplot as plt

# knowing that this data is stored in a 3D matrix (nvarg, ndir, nlag+2)

#plot variogram 1 , direction 1,2,3

v=1

ndir

for d in range(ndir):
    dip=parameters['dip'][d]
    azm=parameters['azm'][d]
    plt.plot (pdis[v, d, 1:], pgam[v, d, 1:], '-o', label=str(dip) + '-->' + str(azm))

plt.legend()
plt.grid(True)

plt.show()<jupyter_output><empty_output><jupyter_text>Variogram map  > **Warning:** remove undefined lags (replace zero with nan?)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

#get an array of direction + complementary direction = direction + 180
ndir2= len(parameters['azm'])*2
azimuths = np.zeros (ndir2 + 1) #this is to repeat the first direction
zeniths = pdis[0, 0, 1:]
values = np.zeros ((ndir2 +1  ,len(pdis[0, 0, 1:])))

for i in range (ndir2/2):
    azimuths[i] = parameters['azm'][i]
    azimuths[i+ndir2/2] = parameters['azm'][i]+180
    for ii in range(len(pdis[0, 0, 1:])):
        values[i,ii] = pgam[0, i, ii+1]
        values[i+ndir2/2,ii] = pgam[0, i, ii+1]        

azimuths[ndir2] = azimuths[ndir2/2]+180
values[ndir2,:] = values[0,:]
        
#prepare grid        
r, theta = np.meshgrid(zeniths, azimuths)

#-- Plot... ------------------------------------------------
ax = plt.subplot(111, polar=True)
# make clockwise like maps
ax.set_theta_direction(-1)
# make the plot to point north 
ax.set_theta_zero_location("N")
cont=ax.contourf(np.deg2rad(theta), r, values)
ax.contour(np.deg2rad(theta), r, values, colors='k')
ax.plot(np.deg2rad(theta.flatten()), r.flatten(), color= '0.6', linewidth=0.2)
ax.set_rmax(50)
plt.colorbar(cont)
plt.show() 

#-- Plot... with zoom ------------------------------------------------
ax = plt.subplot(111, polar=True)
# make clockwise like maps
ax.set_theta_direction(-1)
# make the plot to point north 
ax.set_theta_zero_location("N")
cont=ax.contourf(np.deg2rad(theta), r, values)
ax.contour(np.deg2rad(theta), r, values, colors='k')
ax.plot(np.deg2rad(theta.flatten()), r.flatten(), color= '0.6', linewidth=0.2)
ax.set_rmax(10)
plt.colorbar(cont)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16018306636155608, 'fraction_numerical': 0.03768115942028986, 'mean_word_length': 2.528525296017223, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '551672', 'license': 'MIT', 'n_tokens_mistral': 2546, 'n_tokens_neox': 2413, 'n_words': 627, 'path': 'pygslib/Ipython_templates/broken/gam_raw_variongram_map_horizontal.ipynb', 'repo_name': 'cmrajan/pygslib', 'size': '6555'}