GeostatsPy - Geostatistical Analysis
Quick Reference
import geostatspy.GSLIB as GSLIB
import geostatspy.geostats as geostats
import pandas as pd
df = pd.read_csv('data.csv')
df['npor'], tvpor, tnspor = geostats.nscore(df, 'porosity') # Transform
lag, gamma, npairs = geostats.gamv(df, 'X', 'Y', 'npor', # Variogram
tmin=-9999, tmax=9999, xlag=50, xltol=25, nlag=15,
azm=0, atol=22.5, bandwh=9999, bandwd=9999)
vario = GSLIB.make_variogram(nug=0.0, nst=1, it1=1, cc1=1.0, # Model
azi1=0, hmaj1=300, hmin1=300)
est, var = geostats.kb2d(df, 'X', 'Y', 'npor', ..., vario=vario) # Krige
Key Functions
| Category |
Functions |
| Visualization |
locmap, pixelplt, hist |
| Variogram |
gamv, vmodel |
| Kriging |
kb2d, kb3d |
| Simulation |
sgsim, sisim |
| Transforms |
nscore, backtr |
| Declustering |
declus |
Common Operations
1. Normal Score Transform
df['npor'], tvpor, tnspor = geostats.nscore(df, 'porosity')
original = geostats.backtr(nscore_data, tvpor, tnspor, zmin=0, zmax=0.3)
2. Experimental Variogram
lag, gamma, npairs = geostats.gamv(
df, 'X', 'Y', 'npor',
tmin=-9999, tmax=9999, # Trimming limits
xlag=50, xltol=25, # Lag distance, tolerance
nlag=15, azm=0, atol=22.5, bandwh=9999, bandwd=9999)
3. Variogram Model
# Types: 1=spherical, 2=exponential, 3=gaussian
vario = GSLIB.make_variogram(
nug=0.0, nst=1, # Nugget, number of structures
it1=1, cc1=1.0, # Type, sill contribution
azi1=0, hmaj1=300, hmin1=300) # Azimuth, major/minor range
4. Kriging (kb2d)
est, var = geostats.kb2d(
df, 'X', 'Y', 'npor', tmin=-9999, tmax=9999,
nx=50, xmn=25, xsiz=50, # Grid X: ncells, origin, size
ny=50, ymn=25, ysiz=50, # Grid Y
nxdis=1, nydis=1, ndmin=1, ndmax=10,
radius=500, ktype=0, skmean=0.0, vario=vario) # ktype: 0=simple, 1=ordinary
5. Sequential Gaussian Simulation
sim = geostats.sgsim(
df, 'X', 'Y', 'npor', wcol=-1, scol=-1,
tmin=-9999, tmax=9999, itrans=0,
ismooth=0, dession=0, dmession=0,
zmin=-4, zmax=4, ltail=1, ltpar=0, utail=1, utpar=0,
nsim=1, nx=50, xmn=25, xsiz=50, ny=50, ymn=25, ysiz=50,
nz=1, zmn=0, zsiz=1, seed=73073,
ndmin=1, ndmax=10, nodmax=10, radius=500, radius1=500,
sang1=0, sang2=0, sang3=0, mxctx=10, mxcty=10, mxctz=1,
ktype=0, vario=vario)
6. Declustering
wts, cell_size, ncut = geostats.declus(
df, 'X', 'Y', 'porosity', iminmax=1, noff=10, ncell=20, cmin=10, cmax=500)
declustered_mean = np.average(df['porosity'], weights=wts)
Variogram Models
| Code |
Model |
Use Case |
| 1 |
Spherical |
Most common, finite range |
| 2 |
Exponential |
Reaches sill asymptotically |
| 3 |
Gaussian |
Very smooth, parabolic near origin |
| 4 |
Power |
Unbounded, fractal-like |
Key Parameters
| Parameter |
Description |
nug |
Nugget effect (measurement error + micro-scale variation) |
sill |
Total variance (nugget + structure contributions) |
range |
Distance where correlation becomes negligible |
azimuth |
Direction of maximum continuity (degrees from N) |
ktype |
0=simple kriging (known mean), 1=ordinary kriging |
When to Use vs Alternatives
| Use Case |
Tool |
Why |
| GSLIB-style workflows |
GeostatsPy |
Direct port of GSLIB programs to Python |
| SGSIM / SISIM simulation |
GeostatsPy |
Full GSLIB simulation engine |
| Declustering spatial data |
GeostatsPy |
Built-in declus function |
| Modern variogram API |
scikit-gstat |
Cleaner API, sklearn integration |
| Kriging only (no simulation) |
pykrige |
Focused API, universal kriging support |
| Random field generation |
gstools |
Flexible covariance models, field generation |
| Large-scale 3D geomodelling |
SGeMS / Petrel |
GUI-based, industrial workflows |
| Indicator simulation |
GeostatsPy (sisim) |
Categorical property simulation |
Choose GeostatsPy when: You need GSLIB-compatible workflows in Python, especially
for sequential simulation (SGSIM/SISIM), declustering, or if you are familiar with
GSLIB parameter conventions. Best for reservoir characterization workflows.
Choose scikit-gstat when: You prefer a modern scikit-learn-style API for variogram
analysis and kriging, with better integration into Python data science workflows.
Choose pykrige when: You only need kriging interpolation (no simulation) and want
universal kriging with external drift or regression kriging capabilities.
Common Workflows
Variogram Analysis and Kriging Interpolation
Common Issues
| Issue |
Solution |
| Variogram doesn't reach sill |
Increase nlag or xlag to capture full range |
| Kriging produces negative values |
Back-transform after kriging, not before |
| SGSIM artifacts |
Check grid definition (xmn, xsiz) matches data extent |
| Too few variogram pairs |
Increase atol (angular tolerance) or xltol (lag tolerance) |
| Hole effect in variogram |
May indicate periodicity; try nested structures |
Tips
- Always transform to normal scores - Most methods assume Gaussian
- Start isotropic - Add anisotropy only if justified by directional variograms
- Check variogram pairs - Ensure enough pairs at each lag
- Multiple realizations - Use 50-100+ for uncertainty quantification
- Back-transform last - Apply to final results only
References
- GSLIB Programs - GSLIB program reference
- Simulation Methods - SGSIM, SISIM, and other methods
Scripts
- scripts/kriging_example.py - Complete kriging workflow
majiayu000