Utilities

def vtk_export(filename, pos, fields, fieldnames=None):
    """
    Export spatial fields to VTK format.
    
    Parameters:
    - filename (str): Output VTK filename (without extension)
    - pos (array-like): Spatial positions, shape (dim, n_points)
    - fields (array-like): Field values, shape (n_fields, n_points) or (n_points,)
    - fieldnames (list): Names for field variables (default: auto-generate)
    
    Returns:
    - None: Writes VTK file to disk
    """

def vtk_export_structured(filename, pos, fields, fieldnames=None):
    """
    Export structured grid fields to VTK format.
    
    Parameters:
    - filename (str): Output VTK filename
    - pos (list): Grid coordinate arrays [x, y, z]
    - fields (array-like): Field values on structured grid
    - fieldnames (list): Field variable names
    
    Returns:
    - None: Writes structured VTK file
    """

def vtk_export_unstructured(filename, pos, fields, fieldnames=None):
    """
    Export unstructured point data to VTK format.
    
    Parameters:
    - filename (str): Output VTK filename
    - pos (array-like): Point positions
    - fields (array-like): Field values at points
    - fieldnames (list): Field variable names
    
    Returns:
    - None: Writes unstructured VTK file
    """

def to_vtk(pos, fields, fieldnames=None):
    """
    Convert spatial fields to VTK data structure.
    
    Parameters:
    - pos (array-like): Spatial positions
    - fields (array-like): Field values
    - fieldnames (list): Field variable names
    
    Returns:
    - VTK data structure for direct use with PyVista
    """

def to_vtk_structured(pos, fields, fieldnames=None):
    """Convert structured grid to VTK format."""

def to_vtk_unstructured(pos, fields, fieldnames=None):
    """Convert unstructured data to VTK format."""

def generate_grid(x, y=None, z=None, **kwargs):
    """
    Generate coordinate grid for field evaluation.
    
    Parameters:
    - x (array-like): X-coordinate values or grid specification
    - y (array-like): Y-coordinate values (default: None for 1D)
    - z (array-like): Z-coordinate values (default: None for 2D)
    - **kwargs: Additional grid parameters
    
    Returns:
    - Grid coordinate arrays for structured field generation
    """

def generate_st_grid(space_grid, time_grid, **kwargs):
    """
    Generate spatio-temporal coordinate grid.
    
    Parameters:
    - space_grid (list): Spatial coordinate arrays [x, y, z]
    - time_grid (array-like): Temporal coordinate values
    - **kwargs: Additional grid parameters
    
    Returns:
    - Spatio-temporal coordinate arrays
    """

def rotated_main_axes(dim, angles):
    """
    Calculate rotated coordinate system main axes.
    
    Parameters:
    - dim (int): Spatial dimension (2 or 3)
    - angles (float or array-like): Rotation angles
    
    Returns:
    - Rotation matrix for coordinate transformation
    """

def matrix_rotate(dim, angles):
    """Create rotation matrix."""

def matrix_isometrize(dim, anis, angles):
    """Create isometrization matrix for anisotropic scaling."""

def matrix_anisometrize(dim, anis, angles):
    """Create anisometrization matrix."""

def apply(method, field, **kwargs):
    """
    Apply transformation to field values.
    
    Parameters:
    - method (str): Transformation method name
    - field (array-like): Input field values
    - **kwargs: Method-specific parameters
    
    Returns:
    - Transformed field values
    """

def apply_function(function, field, **kwargs):
    """Apply custom function to field values."""

def binary(field, threshold=0.0, upper=1.0, lower=0.0):
    """
    Binary transformation based on threshold.
    
    Parameters:
    - field (array-like): Input field values
    - threshold (float): Threshold value
    - upper (float): Value for field > threshold
    - lower (float): Value for field <= threshold
    
    Returns:
    - Binary field values
    """

def discrete(field, values, thresholds):
    """
    Discrete transformation using multiple thresholds.
    
    Parameters:
    - field (array-like): Input field values
    - values (array-like): Output discrete values
    - thresholds (array-like): Threshold boundaries
    
    Returns:
    - Discretized field values
    """

def boxcox(field, lmbda, shift=0.0):
    """Box-Cox power transformation."""

def zinnharvey(field, conn='high'):
    """Zinn & Harvey transformation for connectivity."""

def normal_force_moments(field, mean=0.0, var=1.0):
    """Force field to have specified moments."""

def normal_to_lognormal(field, mean=1.0, var=1.0):
    """Transform normal field to log-normal distribution."""

def normal_to_uniform(field, a=0.0, b=1.0):
    """Transform normal field to uniform distribution."""

def normal_to_arcsin(field, a=0.0, b=1.0):
    """Transform normal field to arcsine distribution."""

def normal_to_uquad(field, a=0.0, b=1.0):
    """Transform normal field to U-quadratic distribution."""

def confidence_scaling(dim, nugget, var, len_scale, anis, alpha):
    """Calculate confidence scaling for kriging variance."""

def inc_gamma(s, x):
    """
    Incomplete gamma function.
    
    Parameters:
    - s (float): Shape parameter
    - x (array-like): Input values
    
    Returns:
    - Incomplete gamma function values
    """

def inc_gamma_low(s, x):
    """Lower incomplete gamma function."""

def exp_int(s, x):
    """Exponential integral function."""

def inc_beta(a, b, x):
    """Incomplete beta function."""

def tplstable_cor(r, len_scale, hurst, alpha):
    """TPL stable correlation function."""

def tpl_exp_spec_dens(k, len_scale, hurst, len_low):
    """TPL exponential spectral density."""

def tpl_gau_spec_dens(k, len_scale, hurst, len_low):
    """TPL Gaussian spectral density."""

class Normalizer:
    """
    Base class for data normalization.
    
    Methods:
    - normalize(data): Transform data to normalized form
    - denormalize(data): Transform back to original scale  
    - fit(data): Fit normalization parameters to data
    """
    def normalize(self, data): ...
    def denormalize(self, data): ...
    def fit(self, data): ...

class LogNormal(Normalizer):
    """Log-normal distribution normalizer."""

class BoxCox(Normalizer):
    """Box-Cox transformation normalizer."""

class BoxCoxShift(Normalizer):
    """Box-Cox transformation with shift parameter."""

class YeoJohnson(Normalizer):
    """Yeo-Johnson transformation normalizer."""

class Modulus(Normalizer):
    """Modulus transformation normalizer."""

class Manly(Normalizer):
    """Manly transformation normalizer."""

def apply_mean_norm_trend(field, mean, normalizer, trend, pos):
    """Apply mean, normalization, and trend to field."""

def remove_trend_norm_mean(field, mean, normalizer, trend, pos):
    """Remove trend, normalization, and mean from field."""

EARTH_RADIUS: float = 6371.0
"""Earth radius in kilometers for geographic coordinate calculations."""

KM_SCALE: float = 6371.0
"""Kilometer scaling factor for geographic projections."""

DEGREE_SCALE: float = 57.29577951308232
"""Conversion factor from radians to degrees (180/π)."""

RADIAN_SCALE: float = 1.0
"""Radian scaling factor (identity)."""

class RNG:
    """
    Random number generator with seed control.
    
    Parameters:
    - seed (int): Random seed for reproducibility
    - **kwargs: Additional RNG parameters
    """
    def __init__(self, seed=None, **kwargs): ...

class MasterRNG:
    """Master random number generator for coordinated seeding."""

def dist_gen(dist, **kwargs):
    """Distribution generator factory function."""

NUM_THREADS: int = None
"""Number of threads for parallel computation (None = auto-detect)."""

USE_GSTOOLS_CORE: bool = True  
"""Whether to use GSTools core library for acceleration."""

import gstools as gs
import numpy as np

# Generate field data
model = gs.Gaussian(dim=2, var=1.0, len_scale=10.0)
srf = gs.SRF(model, seed=12345)

# Structured grid export
x = np.arange(0, 50, 1.0)
y = np.arange(0, 30, 1.0)
field = srf.structured([x, y])

gs.vtk_export_structured('field_structured', [x, y], field, ['temperature'])

# Unstructured point export
pos = np.random.rand(2, 1000) * 50
field = srf.unstructured(pos)

gs.vtk_export_unstructured('field_points', pos, field, ['temperature'])

# Simple 1D grid
x = gs.generate_grid(np.linspace(0, 100, 101))

# 2D structured grid
x = np.arange(0, 50, 1.0)
y = np.arange(0, 30, 1.0)
grid_2d = gs.generate_grid(x, y)

# 3D grid
z = np.arange(0, 20, 1.0)
grid_3d = gs.generate_grid(x, y, z)

# Spatio-temporal grid
time = np.linspace(0, 10, 11)
st_grid = gs.generate_st_grid([x, y], time)

# 2D rotation matrix
angles = np.pi/4  # 45 degrees
rot_matrix = gs.rotated_main_axes(2, angles)

# Transform coordinates
coords = np.array([[1, 0], [0, 1]])  # Unit vectors
rotated_coords = rot_matrix @ coords.T

# Anisotropic transformation
anis = [1.0, 0.5]  # Anisotropy ratios
angles = np.pi/6   # 30 degrees
anis_matrix = gs.matrix_isometrize(2, anis, angles)

# Import transform functions
from gstools.transform import normal_to_lognormal, binary, boxcox
# Or: import gstools.transform as gst

# Generate normal field
field = srf.structured([x, y])

# Log-normal transformation
lognormal_field = normal_to_lognormal(field, mean=2.0, var=1.0)

# Binary transformation
binary_field = binary(field, threshold=0.5)

# Box-Cox transformation
boxcox_field = boxcox(field, lmbda=0.5)

# Custom transformation using normalizer
from gstools.normalizer import LogNormal
normalizer = LogNormal()
normalizer.fit(field)
normalized = normalizer.normalize(field)
denormalized = normalizer.denormalize(normalized)

# Generate categorical field
categories = ['sand', 'clay', 'rock']
thresholds = [-0.5, 0.5]

from gstools.transform import discrete
categorical = discrete(
    field, 
    values=categories, 
    thresholds=thresholds
)

# Count category frequencies
unique, counts = np.unique(categorical, return_counts=True)
for cat, count in zip(unique, counts):
    print(f"{cat}: {count} ({100*count/len(field):.1f}%)")

# Incomplete gamma function for Matérn covariance
nu = 1.5
distances = np.linspace(0.1, 10, 100)
gamma_values = gs.tools.inc_gamma(nu, distances)

# Exponential integral
exp_int_values = gs.tools.exp_int(1, distances)

# TPL correlation function
tpl_corr = gs.tools.tplstable_cor(
    distances, len_scale=5.0, hurst=0.5, alpha=1.5
)

# Set number of threads globally
gs.config.NUM_THREADS = 4

# Use GSTools core for acceleration
gs.config.USE_GSTOOLS_CORE = True

# Check current configuration
print(f"Threads: {gs.config.NUM_THREADS}")
print(f"Use core: {gs.config.USE_GSTOOLS_CORE}")

# Work with lat/lon coordinates
lat = np.linspace(40, 45, 50)
lon = np.linspace(-75, -70, 60)

# Use Earth radius for distance calculations
earth_radius_km = gs.EARTH_RADIUS
degree_to_km = earth_radius_km * np.pi / 180

# Convert degrees to kilometers
lat_km = lat * degree_to_km
lon_km = lon * degree_to_km

# Geographic field generation
geo_model = gs.Exponential(
    dim=2, 
    var=1.0, 
    len_scale=100.0,  # 100 km correlation length
    latlon=True
)
geo_srf = gs.SRF(geo_model)
geo_field = geo_srf.structured([lon, lat])

# Import transform functions
from gstools.transform import normal_force_moments, zinnharvey, apply_function

# Force specific moments
field_moments = normal_force_moments(
    field, mean=5.0, var=2.0
)

# Zinn-Harvey connectivity transformation
connected_field = zinnharvey(field, conn='high')

# Apply custom function
def custom_transform(x):
    return np.exp(x) / (1 + np.exp(x))  # Sigmoid

transformed = apply_function(custom_transform, field)

# Complete geostatistical workflow using utilities
import gstools as gs
import numpy as np

# 1. Generate synthetic data
model = gs.Matern(dim=2, var=2.0, len_scale=15.0, nu=1.5)
srf = gs.SRF(model, seed=42)

# 2. Create evaluation grid
x = np.arange(0, 100, 2.0)
y = np.arange(0, 80, 2.0)
field = srf.structured([x, y])

# 3. Apply transformation
from gstools.transform import normal_to_lognormal
log_field = normal_to_lognormal(field, mean=3.0, var=1.0)

# 4. Export for visualization
gs.vtk_export_structured(
    'workflow_result', 
    [x, y], 
    [field, log_field], 
    ['original', 'lognormal']
)

# 5. Generate summary statistics
print(f"Original field - Mean: {np.mean(field):.2f}, Std: {np.std(field):.2f}")
print(f"Log-normal field - Mean: {np.mean(log_field):.2f}, Std: {np.std(log_field):.2f}")