Walker Lake demoΒΆ
This block of code is a temporary hack to show plotly plots inside the online documentation
[1]:
from IPython.display import HTML
HTML('<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js"></script>')
[1]:
This block of code supress scrolling bars
[2]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
Importing necessary packages
[3]:
import pandas as pd
import pyLPM
Importing dataΒΆ
Importing data using datasets module
[4]:
wl = pyLPM.datasets.dataset('Walker_Lake')
Showing first 5 lines of the DataFrame
[5]:
wl.head()
[5]:
| X | Y | Cu | Au | |
|---|---|---|---|---|
| 0 | 11.0 | 8.0 | 0.000 | -999.0 |
| 1 | 8.0 | 30.0 | 0.000 | -999.0 |
| 2 | 9.0 | 48.0 | 2.244 | -999.0 |
| 3 | 8.0 | 69.0 | 4.344 | -999.0 |
| 4 | 9.0 | 90.0 | 4.121 | -999.0 |
Describing dataset
[6]:
wl = pd.read_csv(
filepath_or_buffer='data/walker.csv',
na_values=-999)
Plotting locations mapsΒΆ
[7]:
pyLPM.plots.locmap(
x=wl.X,
y=wl.Y,
variable=wl['T'],
categorical=True,
title='Map for T')
[8]:
pyLPM.plots.locmap(
x=wl.X,
y=wl.Y,
variable=wl.V,
title='Map for V',
pointsize=9)
[9]:
pyLPM.plots.locmap(
x=wl.X,
y=wl.Y,
variable=wl.U,
title='Map for U',
colorscale='Bluered')
Plotting histogramsΒΆ
[10]:
pyLPM.plots.histogram(
wl.V,
title='histogram for V',
n_bins=20)
[11]:
pyLPM.plots.histogram(
wl.U,
cdf=True,
n_bins=20,
title='histogram for U')
Plotting a scatterplotΒΆ
[12]:
pyLPM.plots.scatter2d(
x=wl.V,
y=wl.U,
title='Scatterplot UxV',
xy_line=True)
[13]:
pyLPM.plots.scatter2d(
x=wl.V,
y=wl.U,
variable=wl["T"],
title='Scatterplot UxV',
xy_line=False,
best_fit_line=False)
Plotting a qq plotΒΆ
[14]:
pyLPM.plots.qqplot(
x=wl.V,
y=wl.U,
title='qq plot',
x_axis='V',
y_axis='U',
pointsize=8,
figsize=(700,700),
dicretization = 200)
Running declus from GSLibΒΆ
[15]:
pyLPM.gslib.declus(
df=wl,
x='X',
y='Y',
z=None,
var='V',
x_anis=1.0,
z_anis=1.0,
n_cell=50,
min_size=5.0,
max_size=45.0,
keep_min = True,
number_offsets=100,
usewine=False)
DECLUS Version: 2.906
data file = C:\Users\Roberto Rolo\Documents\pyLPM\py
columns = 1 2 0 3
tmin,tmax = -1.000000E+21 1.000000E+21
summary file = C:\Users\Roberto Rolo\Documents\pyLPM\py
output file = C:\Users\Roberto Rolo\Documents\pyLPM\py
anisotropy = 1.000000 1.000000
minmax flag = -1
ncell min max = 50 5.000000 45.000000
offsets = 100
There are 470 data with:
mean value = 435.29860
minimum and maximum = .00000 1528.10000
size of data vol in X = 243.00000
size of data vol in Y = 283.00000
size of data vol in Z = .00000
declustered mean = 290.10570
min and max weight = .28911 2.78945
equal weighting = 1.00000
DECLUS Version: 2.906 Finished
Stop - Program terminated.
Creating a gridΒΆ
[16]:
grid = pyLPM.utils.autogrid(
x=wl.X,
y=wl.Y,
z=None,
sx=1,
sy=1,
sz=1,
buffer=0)
Running NN from Geostat algorithms moduleΒΆ
[17]:
resultsnn = pyLPM.geostatalgo.nn(
x=wl.X,
y=wl.Y,
z=None,
var=wl.V,
grid=grid)
[18]:
pyLPM.plots.pixelplot(
grid,
resultsnn,
title='NN results')
Plotting the declustered histogramsΒΆ
[19]:
wl.head()
[19]:
| Id | X | Y | V | U | T | Declustering Weight | |
|---|---|---|---|---|---|---|---|
| 0 | 1.0 | 11.0 | 8.0 | 0.0 | -999.0 | 2.0 | 2.771450 |
| 1 | 2.0 | 8.0 | 30.0 | 0.0 | -999.0 | 2.0 | 2.568608 |
| 2 | 3.0 | 9.0 | 48.0 | 224.4 | -999.0 | 2.0 | 1.948161 |
| 3 | 4.0 | 8.0 | 69.0 | 434.4 | -999.0 | 2.0 | 1.509497 |
| 4 | 5.0 | 9.0 | 90.0 | 412.1 | -999.0 | 2.0 | 1.371034 |
[20]:
pyLPM.plots.histogram(
wl.V,
title='Declustered histogram for V',
wt=wl["Declustering Weight"])
[21]:
pyLPM.plots.histogram(
resultsnn,
title='NN Declustered histogram for V')
VariographyΒΆ
Plotting a H-scatterΒΆ
[22]:
pyLPM.gammapy.hscatterplots(
dataset= wl,
X= 'X',
Y= 'Y',
head= 'V',
tail= 'V',
lagsize= 5.,
nlags=3,
azimuth= 160.,
dip= 0.,
lineartol= 2.5,
htol= 22,
vtol= 22,
hband= 15.,
vband= 15.,
Z = None,
choice = 1.0)
[23]:
pyLPM.gammapy.varmap(
dataset=wl,
X='X',
Y='Y',
head='V',
tail='V',
azimuth=0.,
dip=0.,
nlags=10,
lagdistance=10,
ndirections=12,
type_var='Variogram',
Z=None,
choice=1.0,
)
[24]:
exp_variogram = pyLPM.gammapy.experimental(
dataset=wl,
X='X',
Y='Y',
head='V',
tail='V',
type_c=['Variogram', 'Variogram'],
nlags=[15,15],
lagsize=[10,5],
lineartol=[5,2.5],
htol=[22,22],
vtol=[22,22],
hband=[10,10],
vband=[10,10],
azimuth=[166,76],
dip=[0,0],
omni=[False,False],
show_pairs=True,
Z=None,
choice=1.0,)
[25]:
var_model = pyLPM.gammapy.modelling(
experimental_dict=exp_variogram,
rotation_max=166.,
rotation_med=0.,
rotation_min=0.,
nugget=22000,
inverted=False,
rangemax=[30.,150.],
rangemed=[25.,50],
rangemin=[0.,0.],
model=['Spherical','Spherical'],
contribution=[40000.,45000.],
show_pairs=False,)
Cross validationΒΆ
[26]:
pyLPM.gslib.kt3d(
df=wl,
dh='Id',
x='X',
y='Y',
z=None,
var='V',
grid=grid,
variogram=var_model,
min_samples=4,
max_samples=18,
max_oct=3,
search_radius=[150,50,0],
search_ang = [166,0,0],
discretization = [5,5,1],
krig_type='OK',
sk_mean = 0,
tmin=-1.0e21,
tmax=1.0e21,
option='cross',
debug_level=0,
usewine=False)
KT3D Version: 2.907
data file = C:\Users\Roberto Rolo\Documents\pyLPM\py
columns = 1 2 3 0 4
0
trimming limits = -1.000000E+21 1.000000E+21
kriging option = 1
jackknife data file = C:\Users\Roberto Rolo\Documents\pyLPM\py
columns = 2 3 0 4 0
debugging level = 0
debugging file = C:\Users\Roberto Rolo\Documents\pyLPM\py
output file = C:\Users\Roberto Rolo\Documents\pyLPM\py
nx, xmn, xsiz = 244 8.000000 1.000000
ny, ymn, ysiz = 284 8.000000 1.000000
nz, zmn, zsiz = 1 0.000000E+00 1.000000
block discretization: 5 5 1
ndmin,ndmax = 4 18
max per octant = 3
search radii = 150.000000 50.000000 0.000000E+00
search anisotropy angles = 166.000000 0.000000E+00 0.000000E+00
ktype, skmean = 0 0.000000E+00
drift terms = 0 0 0 0 0
0 0 0 0
itrend = 0
external drift file = extdrift.dat
variable in external drift file = 4
nst, c0 = 2 22000.000000
it,cc,ang[1,2,3]; 1 40000.000000 166.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 30.000000 25.000000 0.000000E+00
it,cc,ang[1,2,3]; 1 45000.000000 166.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 150.000000 50.000000 0.000000E+00
Data for KT3D: Variable number 4
Number = 470
Average = 435.298600
Variance = 89738.140000
Setting up rotation matrices for variogram and search
Setting up super block search strategy
Working on the kriging
currently on estimate 47
currently on estimate 94
currently on estimate 141
currently on estimate 188
currently on estimate 235
currently on estimate 282
currently on estimate 329
currently on estimate 376
currently on estimate 423
currently on estimate 470
KT3D Version: 2.907 Finished
Stop - Program terminated.
Creating another gridΒΆ
[27]:
grid = pyLPM.utils.autogrid(
x=wl.X,
y=wl.Y,
z=None,
sx=10,
sy=10,
sz=1,
buffer=0)
KrigingΒΆ
[28]:
estimation, variance = pyLPM.gslib.kt3d(
df=wl,
dh='Id',
x='X',
y='Y',
z=None,
var='V',
grid=grid,
variogram=var_model,
min_samples=4,
max_samples=18,
max_oct=3,
search_radius=[150,50,0],
search_ang = [166,0,0],
discretization = [5,5,1],
krig_type='OK',
sk_mean = 0,
tmin=-1.0e21,
tmax=1.0e21,
option='grid',
debug_level=0,
usewine=False)
KT3D Version: 2.907
data file = C:\Users\Roberto Rolo\Documents\pyLPM\py
columns = 1 2 3 0 4
0
trimming limits = -1.000000E+21 1.000000E+21
kriging option = 0
jackknife data file = C:\Users\Roberto Rolo\Documents\pyLPM\py
columns = 2 3 0 4 0
debugging level = 0
debugging file = C:\Users\Roberto Rolo\Documents\pyLPM\py
output file = C:\Users\Roberto Rolo\Documents\pyLPM\py
nx, xmn, xsiz = 25 8.000000 10.000000
ny, ymn, ysiz = 29 8.000000 10.000000
nz, zmn, zsiz = 1 0.000000E+00 1.000000
block discretization: 5 5 1
ndmin,ndmax = 4 18
max per octant = 3
search radii = 150.000000 50.000000 0.000000E+00
search anisotropy angles = 166.000000 0.000000E+00 0.000000E+00
ktype, skmean = 0 0.000000E+00
drift terms = 0 0 0 0 0
0 0 0 0
itrend = 0
external drift file = extdrift.dat
variable in external drift file = 4
nst, c0 = 2 22000.000000
it,cc,ang[1,2,3]; 1 40000.000000 166.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 30.000000 25.000000 0.000000E+00
it,cc,ang[1,2,3]; 1 45000.000000 166.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 150.000000 50.000000 0.000000E+00
Data for KT3D: Variable number 4
Number = 470
Average = 435.298600
Variance = 89738.140000
Setting up rotation matrices for variogram and search
Setting up super block search strategy
Working on the kriging
currently on estimate 72
currently on estimate 144
currently on estimate 216
currently on estimate 288
currently on estimate 360
currently on estimate 432
currently on estimate 504
currently on estimate 576
currently on estimate 648
currently on estimate 720
KT3D Version: 2.907 Finished
Stop - Program terminated.
debug_level = 2 will plot a histogram for the negative weights
[29]:
estimation, variance = pyLPM.gslib.kt3d(
df=wl,
dh='Id',
x='X',
y='Y',
z=None,
var='V',
grid=grid,
variogram=var_model,
min_samples=4,
max_samples=18,
max_oct=3,
search_radius=[150,50,0],
search_ang = [166,0,0],
discretization = [5,5,1],
krig_type='OK',
sk_mean = 0,
tmin=-1.0e21,
tmax=1.0e21,
option='grid',
debug_level=2,
usewine=False)
KT3D Version: 2.907
data file = C:\Users\Roberto Rolo\Documents\pyLPM\py
columns = 1 2 3 0 4
0
trimming limits = -1.000000E+21 1.000000E+21
kriging option = 0
jackknife data file = C:\Users\Roberto Rolo\Documents\pyLPM\py
columns = 2 3 0 4 0
debugging level = 2
debugging file = C:\Users\Roberto Rolo\Documents\pyLPM\py
output file = C:\Users\Roberto Rolo\Documents\pyLPM\py
nx, xmn, xsiz = 25 8.000000 10.000000
ny, ymn, ysiz = 29 8.000000 10.000000
nz, zmn, zsiz = 1 0.000000E+00 1.000000
block discretization: 5 5 1
ndmin,ndmax = 4 18
max per octant = 3
search radii = 150.000000 50.000000 0.000000E+00
search anisotropy angles = 166.000000 0.000000E+00 0.000000E+00
ktype, skmean = 0 0.000000E+00
drift terms = 0 0 0 0 0
0 0 0 0
itrend = 0
external drift file = extdrift.dat
variable in external drift file = 4
nst, c0 = 2 22000.000000
it,cc,ang[1,2,3]; 1 40000.000000 166.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 30.000000 25.000000 0.000000E+00
it,cc,ang[1,2,3]; 1 45000.000000 166.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 150.000000 50.000000 0.000000E+00
Data for KT3D: Variable number 4
Number = 470
Average = 435.298600
Variance = 89738.140000
Setting up rotation matrices for variogram and search
Setting up super block search strategy
Working on the kriging
currently on estimate 72
currently on estimate 144
currently on estimate 216
currently on estimate 288
currently on estimate 360
currently on estimate 432
currently on estimate 504
currently on estimate 576
currently on estimate 648
currently on estimate 720
stimated 725 blocks
average 281.5421
variance *********
KT3D Version: 2.907 Finished
Stop - Program terminated.
Plotting the estimad gradesΒΆ
[30]:
pyLPM.plots.pixelplot(
grid,
estimation,
title='Kriging results')
[31]:
pyLPM.plots.pixelplot(
grid,
variance,
title='Kriging variance')
Validating the modelΒΆ
[32]:
pyLPM.plots.pixelplot(
grid,
estimation,
title='Mapa estimado',
points=[wl.X, wl.Y, wl.V])
[33]:
pyLPM.plots.qqplot(
wl.V,
estimation,
title='qq plot samples vs estimate')
[34]:
pyLPM.plots.swath_plots(
x=wl.X,
y=wl.Y,
z=None,
point_var=wl.V,
grid=grid,
grid_var=estimation,
n_bins=20)