Update GSTools/mcp_output/mcp_plugin/mcp_service.py

GSTools/mcp_output/mcp_plugin/mcp_service.py

@@ -1,154 +1,1206 @@
+import os
+import sys
+from typing import List, Dict, Any, Optional
+
+# Path settings to include the local source directory
+source_path = os.path.join(os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))), "source", "src")
+if source_path not in sys.path:
+    sys.path.insert(0, source_path)
+
 from fastmcp import FastMCP
+import numpy as np
+
+import gstools as gs
+from gstools import (
+    Gaussian, Exponential, Matern, Spherical, Linear,
+    Stable, Rational, Cubic, Circular, Nugget,
+    SRF, CondSRF, vario_estimate, standard_bins, Krige,
+    generate_grid, rotated_main_axes
+)
+from gstools.transform import (
+    normal_to_lognormal, normal_to_uniform, 
+    boxcox, zinnharvey, binary, discrete
+)
+from gstools.normalizer import LogNormal, BoxCox, YeoJohnson
 
 # Create the FastMCP service application
 mcp = FastMCP("gstools_service")
 
-@mcp.tool(name="list_available_tools", description="List all available GSTools functionalities")
-def list_available_tools() -> dict:
+# Mapping of model names to classes
+COVARIANCE_MODELS = {
+    "gaussian": Gaussian,
+    "exponential": Exponential,
+    "matern": Matern,
+    "spherical": Spherical,
+    "linear": Linear,
+    "stable": Stable,
+    "rational": Rational,
+    "cubic": Cubic,
+    "circular": Circular,
+    "nugget": Nugget,
+}
+
+
+# ===================== Covariance Model Tools =====================
+
+@mcp.tool(name="list_covariance_models", description="List available covariance model types")
+def list_covariance_models() -> dict:
     """
-    List all available functionalities in GSTools.
+    List all available covariance model types in GSTools.
 
     Returns:
-    - dict: A dictionary with success status and list of available tools.
+    - dict: Available model types and their descriptions.
     """
     try:
-        tools = [
-            "random_field",
-            "cov_model",
-            "variogram",
-            "vector_field",
-            "kriging",
-            "conditioned_fields",
-            "transformations",
-            "geo_coordinates",
-            "spatio_temporal",
-            "normalizer",
-            "plurigaussian",
-            "sum_model"
-        ]
-        return {
-            "success": True,
-            "tools": tools,
-            "count": len(tools)
+        result = {
+            "gaussian": "Gaussian covariance model - smooth random fields",
+            "exponential": "Exponential covariance model - rough random fields",
+            "matern": "Matérn covariance model - adjustable smoothness via nu parameter",
+            "spherical": "Spherical covariance model - linear near origin, finite range",
+            "linear": "Linear covariance model - simple linear decrease",
+            "stable": "Stable covariance model - generalization of Gaussian",
+            "rational": "Rational quadratic covariance model",
+            "cubic": "Cubic covariance model - smooth near origin",
+            "circular": "Circular covariance model",
+            "nugget": "Nugget effect model - pure discontinuity at origin",
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="list_covariance_models", description="List all available covariance models in GSTools")
-def list_covariance_models() -> dict:
+
+@mcp.tool(name="create_covariance_model", description="Create a covariance model with specified parameters")
+def create_covariance_model(
+    model_type: str,
+    dim: int = 2,
+    var: float = 1.0,
+    len_scale: float = 10.0,
+    nugget: float = 0.0,
+    nu: Optional[float] = None,
+    alpha: Optional[float] = None
+) -> dict:
     """
-    List all available covariance models in GSTools.
+    Create a covariance model and return its properties.
+
+    Parameters:
+    - model_type (str): Type of covariance model (gaussian, exponential, matern, etc.)
+    - dim (int): Spatial dimension (default: 2)
+    - var (float): Variance/sill (default: 1.0)
+    - len_scale (float): Correlation length scale (default: 10.0)
+    - nugget (float): Nugget effect (default: 0.0)
+    - nu (float, optional): Smoothness parameter for Matérn model
+    - alpha (float, optional): Shape parameter for Stable model
 
     Returns:
-    - dict: A dictionary with success status and list of covariance models.
+    - dict: Model properties including variance, length scale, integral scale, etc.
     """
     try:
-        models = [
-            "Nugget",
-            "Gaussian",
-            "Exponential",
-            "Matern",
-            "Integral",
-            "Stable",
-            "Rational",
-            "Cubic",
-            "Linear",
-            "Circular",
-            "Spherical",
-            "HyperSpherical",
-            "SuperSpherical",
-            "JBessel"
-        ]
-        return {
-            "success": True,
-            "models": models,
-            "count": len(models)
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
+            return {
+                "success": False,
+                "result": None,
+                "error": f"Unknown model type: {model_type}. Available: {list(COVARIANCE_MODELS.keys())}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        
+        # Build kwargs
+        kwargs = {
+            "dim": dim,
+            "var": var,
+            "len_scale": len_scale,
+            "nugget": nugget,
+        }
+        
+        # Add optional parameters for specific models
+        if model_type_lower == "matern" and nu is not None:
+            kwargs["nu"] = nu
+        if model_type_lower == "stable" and alpha is not None:
+            kwargs["alpha"] = alpha
+        
+        model = ModelClass(**kwargs)
+        
+        result = {
+            "model_type": model_type,
+            "dim": model.dim,
+            "var": model.var,
+            "len_scale": float(model.len_scale),
+            "nugget": model.nugget,
+            "sill": model.sill,
+            "integral_scale": float(model.integral_scale),
+            "hankel_kw": model.hankel_kw,
         }
+        
+        # Add model-specific parameters
+        if hasattr(model, "nu"):
+            result["nu"] = model.nu
+        if hasattr(model, "alpha"):
+            result["alpha"] = model.alpha
+            
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-@mcp.tool(name="generate_random_field", description="Generate a spatial random field using GSTools")
-def generate_random_field(model_name: str, size: int, seed: int = None) -> dict:
+
+@mcp.tool(name="evaluate_covariance", description="Evaluate covariance function at given distances")
+def evaluate_covariance(
+    model_type: str,
+    distances: List[float],
+    var: float = 1.0,
+    len_scale: float = 10.0,
+    nugget: float = 0.0,
+    nu: Optional[float] = None
+) -> dict:
     """
-    Generate a spatial random field using a specified covariance model.
+    Evaluate the covariance function at given distances.
 
     Parameters:
-    - model_name: Name of the covariance model (e.g., 'Gaussian', 'Exponential')
-    - size: Size of the random field (e.g., 100 for a 100x100 field)
-    - seed: Random seed for reproducibility (optional)
+    - model_type (str): Type of covariance model
+    - distances (List[float]): List of distances to evaluate
+    - var (float): Variance (default: 1.0)
+    - len_scale (float): Length scale (default: 10.0)
+    - nugget (float): Nugget effect (default: 0.0)
+    - nu (float, optional): Smoothness for Matérn model
 
     Returns:
-    - dict: Information about the generated random field
+    - dict: Covariance and variogram values at each distance
     """
     try:
-        from gstools import SRF, Gaussian, Exponential, Matern
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
+            return {
+                "success": False,
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        kwargs = {"var": var, "len_scale": len_scale, "nugget": nugget}
+        if model_type_lower == "matern" and nu is not None:
+            kwargs["nu"] = nu
+        
+        model = ModelClass(**kwargs)
+        
+        r = np.array(distances)
+        cov_values = model.covariance(r)
+        vario_values = model.variogram(r)
+        
+        result = {
+            "distances": distances,
+            "covariance": cov_values.tolist(),
+            "variogram": vario_values.tolist(),
+            "model_type": model_type,
+            "var": var,
+            "len_scale": len_scale,
+            "nugget": nugget,
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Random Field Generation =====================
 
-        model_map = {
-            "Gaussian": Gaussian,
-            "Exponential": Exponential,
-            "Matern": Matern
+@mcp.tool(name="generate_random_field_1d", description="Generate a 1D spatial random field")
+def generate_random_field_1d(
+    model_type: str,
+    x_min: float,
+    x_max: float,
+    n_points: int,
+    var: float = 1.0,
+    len_scale: float = 10.0,
+    mean: float = 0.0,
+    seed: Optional[int] = None
+) -> dict:
+    """
+    Generate a 1D spatial random field.
+
+    Parameters:
+    - model_type (str): Type of covariance model
+    - x_min (float): Minimum x coordinate
+    - x_max (float): Maximum x coordinate
+    - n_points (int): Number of points
+    - var (float): Variance (default: 1.0)
+    - len_scale (float): Length scale (default: 10.0)
+    - mean (float): Mean of the field (default: 0.0)
+    - seed (int, optional): Random seed for reproducibility
+
+    Returns:
+    - dict: Generated field values and coordinates
+    """
+    try:
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
+            return {
+                "success": False,
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        model = ModelClass(dim=1, var=var, len_scale=len_scale)
+        
+        srf = SRF(model, mean=mean)
+        x = np.linspace(x_min, x_max, n_points)
+        
+        if seed is not None:
+            field = srf((x,), seed=seed)
+        else:
+            field = srf((x,))
+        
+        result = {
+            "x": x.tolist(),
+            "field": field.tolist(),
+            "model_type": model_type,
+            "var": var,
+            "len_scale": len_scale,
+            "mean": mean,
+            "seed": seed,
+            "n_points": n_points,
+            "field_stats": {
+                "min": float(np.min(field)),
+                "max": float(np.max(field)),
+                "mean": float(np.mean(field)),
+                "std": float(np.std(field)),
+            }
         }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="generate_random_field_2d", description="Generate a 2D spatial random field")
+def generate_random_field_2d(
+    model_type: str,
+    x_min: float,
+    x_max: float,
+    y_min: float,
+    y_max: float,
+    nx: int,
+    ny: int,
+    var: float = 1.0,
+    len_scale: float = 10.0,
+    mean: float = 0.0,
+    seed: Optional[int] = None
+) -> dict:
+    """
+    Generate a 2D spatial random field on a structured grid.
 
-        if model_name not in model_map:
+    Parameters:
+    - model_type (str): Type of covariance model
+    - x_min, x_max (float): X coordinate range
+    - y_min, y_max (float): Y coordinate range
+    - nx, ny (int): Number of grid points in x and y
+    - var (float): Variance (default: 1.0)
+    - len_scale (float): Length scale (default: 10.0)
+    - mean (float): Mean of the field (default: 0.0)
+    - seed (int, optional): Random seed
+
+    Returns:
+    - dict: Generated field as 2D array and grid info
+    """
+    try:
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
             return {
                 "success": False,
-                "error": f"Model '{model_name}' not found. Use list_covariance_models to see available options."
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        model = ModelClass(dim=2, var=var, len_scale=len_scale)
+        
+        srf = SRF(model, mean=mean)
+        x = np.linspace(x_min, x_max, nx)
+        y = np.linspace(y_min, y_max, ny)
+        
+        if seed is not None:
+            field = srf.structured((x, y), seed=seed)
+        else:
+            field = srf.structured((x, y))
+        
+        result = {
+            "x": x.tolist(),
+            "y": y.tolist(),
+            "field": field.tolist(),
+            "shape": list(field.shape),
+            "model_type": model_type,
+            "var": var,
+            "len_scale": len_scale,
+            "mean": mean,
+            "seed": seed,
+            "field_stats": {
+                "min": float(np.min(field)),
+                "max": float(np.max(field)),
+                "mean": float(np.mean(field)),
+                "std": float(np.std(field)),
             }
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Variogram Estimation =====================
+
+@mcp.tool(name="estimate_variogram_from_data", description="Estimate empirical variogram from spatial data")
+def estimate_variogram_from_data(
+    pos_x: List[float],
+    pos_y: List[float],
+    values: List[float],
+    n_bins: int = 10,
+    max_lag: Optional[float] = None,
+    estimator: str = "matheron"
+) -> dict:
+    """
+    Estimate an empirical variogram from spatial data.
 
-        model = model_map[model_name](dim=2, var=1, len_scale=10)
-        srf = SRF(model, seed=seed)
-        field = srf((size, size))
+    Parameters:
+    - pos_x (List[float]): X coordinates of data points
+    - pos_y (List[float]): Y coordinates of data points
+    - values (List[float]): Values at each point
+    - n_bins (int): Number of lag bins (default: 10)
+    - max_lag (float, optional): Maximum lag distance
+    - estimator (str): Estimator type - 'matheron' or 'cressie' (default: 'matheron')
 
-        return {
-            "success": True,
-            "model_name": model_name,
-            "field_shape": field.shape,
-            "field": field.tolist()
+    Returns:
+    - dict: Empirical variogram with bin centers and gamma values
+    """
+    try:
+        pos = np.array([pos_x, pos_y])
+        field = np.array(values)
+        
+        # Determine max lag if not provided
+        if max_lag is None:
+            x_range = np.max(pos_x) - np.min(pos_x)
+            y_range = np.max(pos_y) - np.min(pos_y)
+            max_lag = 0.5 * np.sqrt(x_range**2 + y_range**2)
+        
+        # Create bins
+        bin_edges = standard_bins(pos, max_lag=max_lag, bin_no=n_bins)
+        
+        # Estimate variogram
+        bin_center, gamma = vario_estimate(
+            pos, field, bin_edges=bin_edges, estimator=estimator
+        )
+        
+        result = {
+            "bin_center": bin_center.tolist(),
+            "gamma": gamma.tolist(),
+            "n_bins": n_bins,
+            "max_lag": max_lag,
+            "estimator": estimator,
+            "n_points": len(values),
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
-@mcp.tool(name="fit_variogram", description="Fit a variogram model to data")
-def fit_variogram(data: list, model_name: str) -> dict:
+@mcp.tool(name="fit_variogram_model", description="Fit a covariance model to empirical variogram")
+def fit_variogram_model(
+    bin_center: List[float],
+    gamma: List[float],
+    model_type: str = "gaussian",
+    nugget: bool = False
+) -> dict:
     """
-    Fit a variogram model to the given data.
+    Fit a covariance model to an empirical variogram.
 
     Parameters:
-    - data: List of data points to fit the variogram model
-    - model_name: Name of the variogram model (e.g., 'Gaussian', 'Exponential')
+    - bin_center (List[float]): Lag distances (bin centers)
+    - gamma (List[float]): Empirical variogram values
+    - model_type (str): Model type to fit (default: 'gaussian')
+    - nugget (bool): Whether to fit nugget effect (default: False)
 
     Returns:
-    - dict: Fitted variogram model parameters
+    - dict: Fitted model parameters
     """
     try:
-        from gstools.covmodel import fit_variogram
-        from gstools.covmodel.models import Gaussian, Exponential, Matern
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
+            return {
+                "success": False,
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        model = ModelClass(dim=2)
+        
+        bin_center_arr = np.array(bin_center)
+        gamma_arr = np.array(gamma)
+        
+        # Fit the model
+        model.fit_variogram(bin_center_arr, gamma_arr, nugget=nugget)
+        
+        result = {
+            "model_type": model_type,
+            "var": model.var,
+            "len_scale": float(model.len_scale),
+            "nugget": model.nugget,
+            "sill": model.sill,
+            "integral_scale": float(model.integral_scale),
+        }
+        
+        # Add model-specific parameters
+        if hasattr(model, "nu"):
+            result["nu"] = model.nu
+            
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
 
-        model_map = {
-            "Gaussian": Gaussian,
-            "Exponential": Exponential,
-            "Matern": Matern
+
+# ===================== Kriging =====================
+
+@mcp.tool(name="simple_kriging", description="Perform simple kriging interpolation")
+def simple_kriging(
+    cond_pos_x: List[float],
+    cond_pos_y: List[float],
+    cond_values: List[float],
+    target_pos_x: List[float],
+    target_pos_y: List[float],
+    model_type: str = "gaussian",
+    var: float = 1.0,
+    len_scale: float = 10.0,
+    mean: float = 0.0
+) -> dict:
+    """
+    Perform simple kriging interpolation.
+
+    Parameters:
+    - cond_pos_x, cond_pos_y (List[float]): Conditioning point coordinates
+    - cond_values (List[float]): Values at conditioning points
+    - target_pos_x, target_pos_y (List[float]): Target point coordinates
+    - model_type (str): Covariance model type (default: 'gaussian')
+    - var (float): Variance (default: 1.0)
+    - len_scale (float): Length scale (default: 10.0)
+    - mean (float): Known mean for simple kriging (default: 0.0)
+
+    Returns:
+    - dict: Kriging predictions and variances
+    """
+    try:
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
+            return {
+                "success": False,
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        model = ModelClass(dim=2, var=var, len_scale=len_scale)
+        
+        cond_pos = np.array([cond_pos_x, cond_pos_y])
+        cond_val = np.array(cond_values)
+        target_pos = np.array([target_pos_x, target_pos_y])
+        
+        krige = Krige(model, cond_pos=cond_pos, cond_val=cond_val, mean=mean)
+        predictions, variances = krige(target_pos)
+        
+        result = {
+            "predictions": predictions.tolist(),
+            "variances": variances.tolist(),
+            "target_x": target_pos_x,
+            "target_y": target_pos_y,
+            "n_cond_points": len(cond_values),
+            "n_target_points": len(target_pos_x),
+            "model_type": model_type,
+            "var": var,
+            "len_scale": len_scale,
+            "mean": mean,
         }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="ordinary_kriging", description="Perform ordinary kriging interpolation")
+def ordinary_kriging(
+    cond_pos_x: List[float],
+    cond_pos_y: List[float],
+    cond_values: List[float],
+    target_pos_x: List[float],
+    target_pos_y: List[float],
+    model_type: str = "gaussian",
+    var: float = 1.0,
+    len_scale: float = 10.0
+) -> dict:
+    """
+    Perform ordinary kriging interpolation (unknown mean).
+
+    Parameters:
+    - cond_pos_x, cond_pos_y (List[float]): Conditioning point coordinates
+    - cond_values (List[float]): Values at conditioning points
+    - target_pos_x, target_pos_y (List[float]): Target point coordinates
+    - model_type (str): Covariance model type (default: 'gaussian')
+    - var (float): Variance (default: 1.0)
+    - len_scale (float): Length scale (default: 10.0)
 
-        if model_name not in model_map:
+    Returns:
+    - dict: Kriging predictions and variances
+    """
+    try:
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
             return {
                 "success": False,
-                "error": f"Model '{model_name}' not found. Use list_covariance_models to see available options."
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        model = ModelClass(dim=2, var=var, len_scale=len_scale)
+        
+        cond_pos = np.array([cond_pos_x, cond_pos_y])
+        cond_val = np.array(cond_values)
+        target_pos = np.array([target_pos_x, target_pos_y])
+        
+        # Ordinary kriging uses unbiased=True
+        krige = Krige(model, cond_pos=cond_pos, cond_val=cond_val, unbiased=True)
+        predictions, variances = krige(target_pos)
+        
+        result = {
+            "predictions": predictions.tolist(),
+            "variances": variances.tolist(),
+            "target_x": target_pos_x,
+            "target_y": target_pos_y,
+            "n_cond_points": len(cond_values),
+            "n_target_points": len(target_pos_x),
+            "model_type": model_type,
+            "var": var,
+            "len_scale": len_scale,
+            "estimated_mean": float(np.mean(cond_values)),
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Utility Tools =====================
+
+@mcp.tool(name="generate_standard_bins", description="Generate standard binning for variogram estimation")
+def generate_standard_bins(
+    pos_x: List[float],
+    pos_y: List[float],
+    n_bins: int = 10,
+    max_lag: Optional[float] = None
+) -> dict:
+    """
+    Generate standard bin edges for variogram estimation.
+
+    Parameters:
+    - pos_x, pos_y (List[float]): Position coordinates
+    - n_bins (int): Number of bins (default: 10)
+    - max_lag (float, optional): Maximum lag distance
+
+    Returns:
+    - dict: Bin edges and bin centers
+    """
+    try:
+        pos = np.array([pos_x, pos_y])
+        
+        if max_lag is None:
+            x_range = np.max(pos_x) - np.min(pos_x)
+            y_range = np.max(pos_y) - np.min(pos_y)
+            max_lag = 0.5 * np.sqrt(x_range**2 + y_range**2)
+        
+        bin_edges = standard_bins(pos, max_lag=max_lag, bin_no=n_bins)
+        bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+        
+        result = {
+            "bin_edges": bin_edges.tolist(),
+            "bin_centers": bin_centers.tolist(),
+            "n_bins": n_bins,
+            "max_lag": max_lag,
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="get_gstools_constants", description="Get useful constants from GSTools")
+def get_gstools_constants() -> dict:
+    """
+    Get useful constants defined in GSTools for geographic computations.
+
+    Returns:
+    - dict: Earth radius and scale constants
+    """
+    try:
+        result = {
+            "EARTH_RADIUS": float(gs.EARTH_RADIUS),
+            "KM_SCALE": float(gs.KM_SCALE),
+            "DEGREE_SCALE": float(gs.DEGREE_SCALE),
+            "RADIAN_SCALE": float(gs.RADIAN_SCALE),
+            "description": {
+                "EARTH_RADIUS": "Earth radius in km (~6371)",
+                "KM_SCALE": "Scale factor for km on Earth surface",
+                "DEGREE_SCALE": "Scale factor for degrees on Earth surface",
+                "RADIAN_SCALE": "Scale factor for radians (1.0)",
             }
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
 
-        model = model_map[model_name]()
-        params = fit_variogram(data, model)
 
-        return {
-            "success": True,
-            "model_name": model_name,
-            "fitted_parameters": params
+# ===================== Conditioned Random Fields =====================
+
+@mcp.tool(name="generate_conditioned_random_field_2d", description="Generate a conditioned 2D spatial random field")
+def generate_conditioned_random_field_2d(
+    model_type: str,
+    cond_pos_x: List[float],
+    cond_pos_y: List[float],
+    cond_values: List[float],
+    x_min: float,
+    x_max: float,
+    y_min: float,
+    y_max: float,
+    nx: int,
+    ny: int,
+    var: float = 1.0,
+    len_scale: float = 10.0,
+    mean: float = 0.0,
+    seed: Optional[int] = None
+) -> dict:
+    """
+    Generate a 2D conditioned spatial random field (honors conditioning data).
+
+    Parameters:
+    - model_type (str): Type of covariance model
+    - cond_pos_x, cond_pos_y (List[float]): Conditioning point coordinates
+    - cond_values (List[float]): Values at conditioning points
+    - x_min, x_max, y_min, y_max (float): Field extent
+    - nx, ny (int): Number of grid points
+    - var (float): Variance (default: 1.0)
+    - len_scale (float): Length scale (default: 10.0)
+    - mean (float): Mean (default: 0.0)
+    - seed (int, optional): Random seed
+
+    Returns:
+    - dict: Conditioned field and grid info
+    """
+    try:
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
+            return {
+                "success": False,
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        model = ModelClass(dim=2, var=var, len_scale=len_scale)
+        
+        cond_pos = np.array([cond_pos_x, cond_pos_y])
+        cond_val = np.array(cond_values)
+        
+        # Create kriging object
+        krige = Krige(model, cond_pos=cond_pos, cond_val=cond_val, mean=mean)
+        
+        # Create conditioned SRF
+        cond_srf = CondSRF(krige)
+        
+        # Generate field
+        x = np.linspace(x_min, x_max, nx)
+        y = np.linspace(y_min, y_max, ny)
+        
+        if seed is not None:
+            field = cond_srf.structured((x, y), seed=seed)
+        else:
+            field = cond_srf.structured((x, y))
+        
+        result = {
+            "x": x.tolist(),
+            "y": y.tolist(),
+            "field": field.tolist(),
+            "shape": list(field.shape),
+            "cond_pos_x": cond_pos_x,
+            "cond_pos_y": cond_pos_y,
+            "cond_values": cond_values,
+            "n_cond_points": len(cond_values),
+            "model_type": model_type,
+            "var": var,
+            "len_scale": len_scale,
+            "mean": mean,
+            "seed": seed,
+            "field_stats": {
+                "min": float(np.min(field)),
+                "max": float(np.max(field)),
+                "mean": float(np.mean(field)),
+                "std": float(np.std(field)),
+            }
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
-# Add more tools here following the same pattern as above
-# Each tool should wrap a core functionality of GSTools
 
-if __name__ == "__main__":
-    mcp.run()
+# ===================== Field Transformations =====================
+
+@mcp.tool(name="transform_to_lognormal", description="Transform normal field to lognormal")
+def transform_to_lognormal(
+    field_values: List[float],
+    mean: float = 1.0,
+    var: float = 1.0
+) -> dict:
+    """
+    Transform a normal field to a lognormal distribution.
+
+    Parameters:
+    - field_values (List[float]): Normal field values
+    - mean (float): Target lognormal mean (default: 1.0)
+    - var (float): Target lognormal variance (default: 1.0)
+
+    Returns:
+    - dict: Transformed field values and statistics
+    """
+    try:
+        field = np.array(field_values)
+        transformed = normal_to_lognormal(field, mean=mean, var=var)
+        
+        result = {
+            "original_values": field_values,
+            "transformed_values": transformed.tolist(),
+            "target_mean": mean,
+            "target_var": var,
+            "actual_mean": float(np.mean(transformed)),
+            "actual_var": float(np.var(transformed)),
+            "actual_std": float(np.std(transformed)),
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="transform_to_uniform", description="Transform normal field to uniform distribution")
+def transform_to_uniform(
+    field_values: List[float],
+    lower: float = 0.0,
+    upper: float = 1.0
+) -> dict:
+    """
+    Transform a normal field to a uniform distribution.
+
+    Parameters:
+    - field_values (List[float]): Normal field values
+    - lower (float): Lower bound (default: 0.0)
+    - upper (float): Upper bound (default: 1.0)
+
+    Returns:
+    - dict: Transformed field values
+    """
+    try:
+        field = np.array(field_values)
+        transformed = normal_to_uniform(field, lower=lower, upper=upper)
+        
+        result = {
+            "original_values": field_values,
+            "transformed_values": transformed.tolist(),
+            "lower": lower,
+            "upper": upper,
+            "actual_min": float(np.min(transformed)),
+            "actual_max": float(np.max(transformed)),
+            "actual_mean": float(np.mean(transformed)),
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="transform_to_binary", description="Transform field to binary values")
+def transform_to_binary(
+    field_values: List[float],
+    threshold: Optional[float] = None,
+    upper_value: float = 1.0,
+    lower_value: float = 0.0
+) -> dict:
+    """
+    Transform a field to binary values (0/1 or custom values).
+
+    Parameters:
+    - field_values (List[float]): Field values
+    - threshold (float, optional): Threshold value. If None, uses mean.
+    - upper_value (float): Value for above threshold (default: 1.0)
+    - lower_value (float): Value for below threshold (default: 0.0)
+
+    Returns:
+    - dict: Binary field values and statistics
+    """
+    try:
+        field = np.array(field_values)
+        
+        if threshold is None:
+            threshold = np.mean(field)
+        
+        transformed = binary(field, threshold=threshold, 
+                           upper=upper_value, lower=lower_value)
+        
+        n_upper = np.sum(transformed == upper_value)
+        n_lower = np.sum(transformed == lower_value)
+        
+        result = {
+            "original_values": field_values,
+            "transformed_values": transformed.tolist(),
+            "threshold": float(threshold),
+            "upper_value": upper_value,
+            "lower_value": lower_value,
+            "n_upper": int(n_upper),
+            "n_lower": int(n_lower),
+            "upper_ratio": float(n_upper / len(field)),
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="transform_to_discrete", description="Transform field to discrete classes")
+def transform_to_discrete(
+    field_values: List[float],
+    thresholds: List[float],
+    values: Optional[List[float]] = None
+) -> dict:
+    """
+    Transform a field to discrete classes based on thresholds.
+
+    Parameters:
+    - field_values (List[float]): Field values
+    - thresholds (List[float]): Threshold values for binning
+    - values (List[float], optional): Class values. If None, uses [0, 1, 2, ...]
+
+    Returns:
+    - dict: Discretized field and class statistics
+    """
+    try:
+        field = np.array(field_values)
+        thresholds_arr = np.array(thresholds)
+        
+        if values is None:
+            values = list(range(len(thresholds) + 1))
+        
+        transformed = discrete(field, thresholds=thresholds_arr, values=values)
+        
+        # Count values in each class
+        class_counts = {}
+        for val in values:
+            class_counts[float(val)] = int(np.sum(transformed == val))
+        
+        result = {
+            "original_values": field_values,
+            "transformed_values": transformed.tolist(),
+            "thresholds": thresholds,
+            "class_values": values,
+            "n_classes": len(values),
+            "class_counts": class_counts,
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="apply_boxcox_transform", description="Apply Box-Cox transformation")
+def apply_boxcox_transform(
+    field_values: List[float],
+    lmbda: float = 1.0,
+    shift: float = 0.0
+) -> dict:
+    """
+    Apply Box-Cox power transformation to field.
+
+    Parameters:
+    - field_values (List[float]): Field values (must be positive)
+    - lmbda (float): Box-Cox lambda parameter (default: 1.0)
+    - shift (float): Shift to add before transformation (default: 0.0)
+
+    Returns:
+    - dict: Transformed field values
+    """
+    try:
+        field = np.array(field_values)
+        transformed = boxcox(field, lmbda=lmbda, shift=shift)
+        
+        result = {
+            "original_values": field_values,
+            "transformed_values": transformed.tolist(),
+            "lambda": lmbda,
+            "shift": shift,
+            "original_stats": {
+                "mean": float(np.mean(field)),
+                "std": float(np.std(field)),
+            },
+            "transformed_stats": {
+                "mean": float(np.mean(transformed)),
+                "std": float(np.std(transformed)),
+            }
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Normalizers =====================
+
+@mcp.tool(name="create_lognormal_normalizer", description="Create a LogNormal normalizer")
+def create_lognormal_normalizer(
+    mean: float = 1.0,
+    var: float = 1.0
+) -> dict:
+    """
+    Create a LogNormal normalizer configuration.
+
+    Parameters:
+    - mean (float): Target lognormal mean (default: 1.0)
+    - var (float): Target lognormal variance (default: 1.0)
+
+    Returns:
+    - dict: Normalizer parameters
+    """
+    try:
+        normalizer = LogNormal(mean=mean, var=var)
+        
+        result = {
+            "normalizer_type": "LogNormal",
+            "mean": mean,
+            "var": var,
+            "description": "Transforms normal fields to lognormal distribution"
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Grid Generation Tools =====================
+
+@mcp.tool(name="generate_grid_2d", description="Generate a 2D structured grid")
+def generate_grid_2d(
+    x_min: float,
+    x_max: float,
+    y_min: float,
+    y_max: float,
+    nx: int,
+    ny: int
+) -> dict:
+    """
+    Generate a 2D structured grid with GSTools generate_grid.
+
+    Parameters:
+    - x_min, x_max (float): X coordinate range
+    - y_min, y_max (float): Y coordinate range
+    - nx, ny (int): Number of grid points
+
+    Returns:
+    - dict: Grid coordinates and meshgrid arrays
+    """
+    try:
+        x = np.linspace(x_min, x_max, nx)
+        y = np.linspace(y_min, y_max, ny)
+        
+        grid_x, grid_y = generate_grid([x, y])
+        
+        result = {
+            "x": x.tolist(),
+            "y": y.tolist(),
+            "grid_x": grid_x.tolist(),
+            "grid_y": grid_y.tolist(),
+            "nx": nx,
+            "ny": ny,
+            "n_points": nx * ny,
+            "x_extent": [x_min, x_max],
+            "y_extent": [y_min, y_max],
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="calculate_rotated_axes", description="Calculate rotated main axes")
+def calculate_rotated_axes(
+    dim: int,
+    angles: List[float]
+) -> dict:
+    """
+    Calculate rotation matrix for anisotropic models.
+
+    Parameters:
+    - dim (int): Spatial dimension (2 or 3)
+    - angles (List[float]): Rotation angles in radians
+
+    Returns:
+    - dict: Rotation matrix and transformed axes
+    """
+    try:
+        angles_arr = np.array(angles)
+        rotation_matrix = rotated_main_axes(dim, angles_arr)
+        
+        result = {
+            "dim": dim,
+            "angles_rad": angles,
+            "angles_deg": (np.array(angles) * 180 / np.pi).tolist(),
+            "rotation_matrix": rotation_matrix.tolist(),
+            "matrix_shape": list(rotation_matrix.shape),
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Advanced Covariance Features =====================
+
+@mcp.tool(name="fit_variogram_with_multiple_models", description="Fit and compare multiple variogram models")
+def fit_variogram_with_multiple_models(
+    bin_center: List[float],
+    gamma: List[float],
+    model_types: Optional[List[str]] = None,
+    nugget: bool = False
+) -> dict:
+    """
+    Fit multiple covariance models and compare their fit quality.
+
+    Parameters:
+    - bin_center (List[float]): Lag distances
+    - gamma (List[float]): Empirical variogram values
+    - model_types (List[str], optional): Models to fit. If None, fits all available.
+    - nugget (bool): Whether to fit nugget effect (default: False)
+
+    Returns:
+    - dict: Fitted parameters for all models and comparison metrics
+    """
+    try:
+        if model_types is None:
+            model_types = ["gaussian", "exponential", "matern", "spherical"]
+        
+        bin_center_arr = np.array(bin_center)
+        gamma_arr = np.array(gamma)
+        
+        results = {}
+        for model_type in model_types:
+            model_type_lower = model_type.lower()
+            if model_type_lower not in COVARIANCE_MODELS:
+                continue
+            
+            try:
+                ModelClass = COVARIANCE_MODELS[model_type_lower]
+                model = ModelClass(dim=2)
+                model.fit_variogram(bin_center_arr, gamma_arr, nugget=nugget)
+                
+                # Calculate fitted values and RMSE
+                fitted_gamma = model.variogram(bin_center_arr)
+                rmse = float(np.sqrt(np.mean((gamma_arr - fitted_gamma)**2)))
+                
+                results[model_type] = {
+                    "var": model.var,
+                    "len_scale": float(model.len_scale),
+                    "nugget": model.nugget,
+                    "rmse": rmse,
+                    "fitted_values": fitted_gamma.tolist(),
+                }
+                
+                if hasattr(model, "nu"):
+                    results[model_type]["nu"] = model.nu
+                    
+            except Exception as e:
+                results[model_type] = {"error": str(e)}
+        
+        # Find best fit
+        best_model = min(
+            [k for k in results if "rmse" in results[k]], 
+            key=lambda k: results[k]["rmse"]
+        ) if any("rmse" in v for v in results.values()) else None
+        
+        result = {
+            "models": results,
+            "best_model": best_model,
+            "best_rmse": results[best_model]["rmse"] if best_model else None,
+            "n_models_tested": len(model_types),
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="evaluate_anisotropic_covariance", description="Evaluate anisotropic covariance model")
+def evaluate_anisotropic_covariance(
+    model_type: str,
+    distances_x: List[float],
+    distances_y: List[float],
+    var: float = 1.0,
+    len_scale_x: float = 10.0,
+    len_scale_y: float = 5.0,
+    angle_deg: float = 0.0
+) -> dict:
+    """
+    Evaluate anisotropic covariance with different length scales in each direction.
+
+    Parameters:
+    - model_type (str): Type of covariance model
+    - distances_x, distances_y (List[float]): Distances in x and y directions
+    - var (float): Variance (default: 1.0)
+    - len_scale_x, len_scale_y (float): Length scales in x and y (default: 10.0, 5.0)
+    - angle_deg (float): Rotation angle in degrees (default: 0.0)
+
+    Returns:
+    - dict: Anisotropic covariance values
+    """
+    try:
+        model_type_lower = model_type.lower()
+        if model_type_lower not in COVARIANCE_MODELS:
+            return {
+                "success": False,
+                "result": None,
+                "error": f"Unknown model type: {model_type}"
+            }
+        
+        ModelClass = COVARIANCE_MODELS[model_type_lower]
+        
+        # Create anisotropic model
+        anis = len_scale_y / len_scale_x  # anisotropy ratio
+        angle_rad = angle_deg * np.pi / 180
+        
+        model = ModelClass(
+            dim=2, 
+            var=var, 
+            len_scale=len_scale_x,
+            anis=anis,
+            angles=angle_rad
+        )
+        
+        # Create position arrays
+        dx = np.array(distances_x)
+        dy = np.array(distances_y)
+        
+        # Calculate distances
+        r = np.sqrt(dx**2 + dy**2)
+        
+        # Evaluate covariance
+        cov_values = model.covariance(np.array([dx, dy]))
+        
+        result = {
+            "model_type": model_type,
+            "var": var,
+            "len_scale_x": len_scale_x,
+            "len_scale_y": len_scale_y,
+            "anisotropy_ratio": anis,
+            "rotation_angle_deg": angle_deg,
+            "distances_x": distances_x,
+            "distances_y": distances_y,
+            "euclidean_distances": r.tolist(),
+            "covariance": cov_values.tolist(),
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+def create_app() -> FastMCP:
+    """
+    Create and return the FastMCP application instance.
+
+    Returns:
+    - The FastMCP instance for the service.
+    """
+    return mcp