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
from plotly.offline import download_plotlyjs, init_notebook_mode
import pyLPM
Importing dataΒΆ
Importing data as a pandas DataFrame
[4]:
wl = pd.read_csv(
'data/walker.csv',
na_values = -999)
Showing first 5 lines of the DataFrame
[5]:
wl.head()
[5]:
| Id | X | Y | V | U | T | |
|---|---|---|---|---|---|---|
| 0 | 1.0 | 11.0 | 8.0 | 0.0 | NaN | 2.0 |
| 1 | 2.0 | 8.0 | 30.0 | 0.0 | NaN | 2.0 |
| 2 | 3.0 | 9.0 | 48.0 | 224.4 | NaN | 2.0 |
| 3 | 4.0 | 8.0 | 69.0 | 434.4 | NaN | 2.0 |
| 4 | 5.0 | 9.0 | 90.0 | 412.1 | NaN | 2.0 |
Plotting locations mapsΒΆ
[6]:
pyLPM.plots.locmap(
x=wl.X,
y=wl.Y,
variable=wl['T'],
categorical=True,
title='Map for T')
[7]:
pyLPM.plots.locmap(
x=wl.X,
y=wl.Y,
variable=wl.V,
title='Map for V',
pointsize=9)
[8]:
pyLPM.plots.locmap(
x=wl.X,
y=wl.Y,
variable=wl.U,
title='Map for U',
colorscale='Bluered')
Plotting histogramsΒΆ
[9]:
pyLPM.plots.histogram(
wl.V,
title='histogram for V',
n_bins=20)
[10]:
pyLPM.plots.histogram(
wl.U,
cdf=True,
n_bins=20,
title='histogram for U')
Plotting a scatterplotΒΆ
[11]:
pyLPM.plots.scatter2d(
x=wl.V,
y=wl.U,
title='Scatterplot UxV',
xy_line=True)
[12]:
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ΒΆ
[13]:
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ΒΆ
[14]:
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 = pyLPM/gslib90/tmp/tmp.dat
columns = 1 2 0 3
tmin,tmax = -1.000000E+21 1.000000E+21
summary file = pyLPM/gslib90/tmp/tmpsum.dat
output file = pyLPM/gslib90/tmp/tmpfile.dat
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ΒΆ
[15]:
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ΒΆ
[16]:
resultsnn = pyLPM.geostatalgo.nn(
x=wl.X,
y=wl.Y,
z=None,
var=wl.V,
grid=grid)
[17]:
pyLPM.plots.pixelplot(
grid,
resultsnn,
title='NN results')
Plotting the declustered histogramsΒΆ
[18]:
wl.head()
[18]:
| 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 |
[19]:
pyLPM.plots.histogram(
wl.V,
title='Declustered histogram for V',
wt=wl["Declustering Weight"])
[20]:
pyLPM.plots.histogram(
resultsnn,
title='NN Declustered histogram for V')
VariographyΒΆ
[21]:
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,
)
[22]:
exp_variogram = pyLPM.gammapy.experimental(
dataset=wl,
X='X',
Y='Y',
head='V',
tail='V',
type_c=['Variogram', 'Variogram'],
nlags=[10,10],
lagsize=[10,10],
lineartol=[5,5],
htol=[22,22],
vtol=[22,22],
hband=[10,10],
vband=[10,10],
azimuth=[160,70],
dip=[0,0],
omni=[False,False],
show_pairs=True,
Z=None,
choice=1.0,)
[23]:
var_model = pyLPM.gammapy.modelling(
experimental_dict=exp_variogram,
rotation_max=160.,
rotation_med=0.,
rotation_min=0.,
nugget=45000,
inverted=False,
rangemax=[70.],
rangemed=[5.],
rangemin=[0.],
model=['Spherical'],
contribution=[40000.],
show_pairs=False,)
Cross validationΒΆ
[24]:
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=8,
max_oct=0,
search_radius=[120,60,1],
search_ang = [0,0,0],
discretization = [5,5,1],
krig_type='OK',
sk_mean = 0,
tmin=-1.0e21,
tmax=1.0e21,
option='cross',
debug_level=0,
debug_file='pyLPM/gslib90/tmp/debug.out',
output_file='pyLPM/gslib90/tmp/output.out',
usewine=False)
KT3D Version: 2.907
data file = pyLPM/gslib90/tmp/tmp.dat
columns = 1 2 3 0 4
0
trimming limits = -1.000000E+21 1.000000E+21
kriging option = 1
jackknife data file = pyLPM/gslib90/tmp/tmp.dat
columns = 2 3 0 4 0
debugging level = 0
debugging file = pyLPM/gslib90/tmp/debug.out
output file = pyLPM/gslib90/tmp/output.out
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 8
max per octant = 0
search radii = 120.000000 60.000000 1.000000
search anisotropy angles = 0.000000E+00 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 = 1 45000.000000
it,cc,ang[1,2,3]; 1 40000.000000 160.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 70.000000 5.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ΒΆ
[25]:
grid = pyLPM.utils.autogrid(
x=wl.X,
y=wl.Y,
z=None,
sx=10,
sy=10,
sz=1,
buffer=0)
KrigingΒΆ
[26]:
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=8,
max_oct=0,
search_radius=[120,60,1],
search_ang = [0,0,0],
discretization = [5,5,1],
krig_type='OK',
sk_mean = 0,
tmin=-1.0e21,
tmax=1.0e21,
option='grid',
debug_level=0,
debug_file='pyLPM/gslib90/tmp/debug.out',
output_file='pyLPM/gslib90/tmp/output.out',
usewine=False)
KT3D Version: 2.907
data file = pyLPM/gslib90/tmp/tmp.dat
columns = 1 2 3 0 4
0
trimming limits = -1.000000E+21 1.000000E+21
kriging option = 0
jackknife data file = pyLPM/gslib90/tmp/tmp.dat
columns = 2 3 0 4 0
debugging level = 0
debugging file = pyLPM/gslib90/tmp/debug.out
output file = pyLPM/gslib90/tmp/output.out
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 8
max per octant = 0
search radii = 120.000000 60.000000 1.000000
search anisotropy angles = 0.000000E+00 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 = 1 45000.000000
it,cc,ang[1,2,3]; 1 40000.000000 160.000000 0.000000E+00
0.000000E+00
a1 a2 a3: 70.000000 5.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.
Plotting the estimad gradesΒΆ
[27]:
pyLPM.plots.pixelplot(
grid,
estimation,
title='Kriging results')
[28]:
pyLPM.plots.pixelplot(
grid,
variance,
title='Kriging variance')
Validating the modelΒΆ
[29]:
pyLPM.plots.pixelplot(
grid,
estimation,
title='Mapa estimado',
points=[wl.X, wl.Y, wl.V])
[30]:
pyLPM.plots.qqplot(
wl.V,
estimation,
title='qq plot samples vs estimate')
[31]:
pyLPM.plots.swath_plots(
x=wl.X,
y=wl.Y,
z=None,
point_var=wl.V,
grid=grid,
grid_var=estimation,
n_bins=20)