Semivariogram parameter fitting/optimization

echopop/analysis.py

@@ -2,7 +2,10 @@
 General analysis orchestration functions that bundle related functions and procedures
 """
 
+import copy
+import os
 import warnings
+from typing import Dict, List, Optional, Tuple, Union
 
 import numpy as np
 import pandas as pd
@@ -30,6 +33,12 @@
     summarize_transect_strata,
     transect_spatial_features,
 )
+from .spatial.variogram import (
+    create_optimization_options,
+    empirical_variogram,
+    optimize_variogram,
+    validate_variogram_parameters,
+)
 from .statistics import stratified_transect_statistic
 
 
@@ -243,6 +252,106 @@ def stratified_summary(
     return stratified_results, analysis_dict
 
 
+def variogram_analysis(
+    variogram_parameters: dict,
+    optimization_parameters: dict,
+    transect_dict: dict,
+    settings_dict: dict,
+    isobath_df: pd.DataFrame,
+    fit_parameters: Union[str, List[float]],
+    initial_values: Optional[Union[List[Tuple[str, float]], Dict[str, Dict[str, float]]]],
+    lower_bounds: Union[List[Tuple[str, float]], Dict[str, Dict[str, float]]],
+    upper_bounds: Optional[Union[List[Tuple[str, float]], Dict[str, Dict[str, float]]]],
+):
+
+    # Extract specific variogram parameters
+    # ---- Number of lags
+    n_lags = variogram_parameters["n_lags"]
+    # ---- Lag resolution
+    lag_resolution = variogram_parameters["lag_resolution"]
+
+    # Compute the lag distances
+    distance_lags = np.arange(1, n_lags) * lag_resolution
+    # ---- Add to the `variogram_parameters` dictionary
+    variogram_parameters["distance_lags"] = distance_lags
+    # ---- Update the max range parameter, if necessary
+    max_range = lag_resolution * n_lags
+    # ---- Mismatched value in default parameters should be replaced
+    if "range" in variogram_parameters and variogram_parameters["range"] != max_range:
+        # ---- Update
+        variogram_parameters["range"] = max_range
+        # ---- Skip file dir in message
+        _warn_skips = (os.path.dirname(__file__),)
+        # ---- Generate warning
+        warnings.warn(
+            f"Default `range` parameter ({variogram_parameters['range']}) does not "
+            f"align with the inputs for `lag_resolution` ({lag_resolution}) "
+            f"and `n_lags` ({n_lags}). The range will be changed to {max_range}.",
+            stacklevel=2,
+            skip_file_prefixes=_warn_skips,
+        )
+
+    # Generate the optimization settings dictionary
+    optimization_settings = create_optimization_options(
+        fit_parameters,
+        variogram_parameters,
+        initial_values=initial_values,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        model=variogram_parameters["model"],
+        **optimization_parameters,
+    )
+
+    # Validate all variogram-related inputs
+    validate_variogram_parameters(variogram_parameters, fit_parameters, optimization_settings)
+
+    # Prepare the transect data
+    # ---- Create a copy of the transect dictionary
+    transect_input = copy.deepcopy(transect_dict)
+    # ---- Edit the transect data
+    transect_data = edit_transect_columns(transect_input, settings_dict)
+
+    # Standardize the transect coordinates, if necessary
+    if settings_dict["standardize_coordinates"]:
+        # ---- Transform geometry
+        transect_data, _, _ = transform_geometry(transect_data, isobath_df, settings_dict)
+        # ---- Print message if verbose
+        if settings_dict["verbose"]:
+            # ---- Print alert
+            print(
+                "Longitude and latitude coordinates (WGS84) converted to standardized "
+                "coordinates (x and y)."
+            )
+    else:
+        # ---- x
+        transect_data["x"] = "longitude"
+        # ---- y
+        transect_data["y"] = "latitude"
+
+    # Compute the empirical variogram
+    lags, gamma_h, lag_counts, lag_covariance = empirical_variogram(
+        transect_data, variogram_parameters, settings_dict
+    )
+
+    # Least-squares fitting
+    best_fit_variogram, initial_fit, optimized_fit = optimize_variogram(
+        lag_counts, lags, gamma_h, variogram_parameters, optimization_settings
+    )
+
+    # Return a dictionary of results
+    return {
+        "best_fit_parameters": best_fit_variogram,
+        "initial_fit": {
+            "parameters": dict(zip(initial_fit[0], initial_fit[1])),
+            "MAD": initial_fit[2],
+        },
+        "optimized_fit": {
+            "parameters": dict(zip(optimized_fit[0], optimized_fit[1])),
+            "MAD": optimized_fit[2],
+        },
+    }
+
+
 def krige(input_dict: dict, analysis_dict: dict, settings_dict: dict) -> tuple[pd.DataFrame, dict]:
     """
     Interpolate spatial data using ordinary kriging.

echopop/core.py

@@ -290,7 +290,7 @@
         },
         "stratified": dict(),
     },
-    "results": {"transect": dict(), "stratified": dict(), "kriging": dict()},
+    "results": {"transect": dict(), "stratified": dict(), "kriging": dict(), "variogram": dict()},
 }
 
 LAYER_NAME_MAP = {

echopop/spatial/mesh.py

@@ -346,7 +346,9 @@ def interpolate_crop_method(
 
 
 def griddify_lag_distances(
-    coordinates_1: Union[pd.DataFrame, np.ndarray], coordinates_2: Union[pd.DataFrame, np.ndarray]
+    coordinates_1: Union[pd.DataFrame, np.ndarray],
+    coordinates_2: Union[pd.DataFrame, np.ndarray],
+    angles: bool = False,
 ):
     """
     Calculate point-to-point distances between two gridded dataframes
@@ -358,6 +360,9 @@ def griddify_lag_distances(
         of values
     coordinates_2: pd.DataFrame
         Georeferenced dataframe
+    angles: bool
+        A boolean flag determining whether the point-to-point azimuth angles are also computed
+        alongside the distances.
 
     Notes
     ----------
@@ -378,8 +383,23 @@ def griddify_lag_distances(
         # ---- Differences across y-coordinates
         y_distance = np.subtract.outer(coordinates_2, coordinates_2)
 
-    # Return Euclidean distances
-    return np.sqrt(x_distance * x_distance + y_distance * y_distance)
+    # Return Euclidean distances (and angles if appropriate)
+    if angles:
+        # ---- Copy x-array
+        x_angles = x_distance.copy()
+        # ---- Replace the self-points with NaN
+        np.fill_diagonal(x_angles, np.nan)
+        # ---- Copy y-array
+        y_angles = y_distance.copy()
+        # ---- Replace the self-points with NaN
+        np.fill_diagonal(y_angles, np.nan)
+        # ---- Calculate the azimuth angle grid
+        angularity = np.arctan(y_angles / x_angles) * 180.0 / np.pi + 180 % 180
+        # ---- Return output
+        return np.sqrt(x_distance * x_distance + y_distance * y_distance), angularity
+    else:
+        # Return Euclidean distances
+        return np.sqrt(x_distance * x_distance + y_distance * y_distance)
 
 
 def stratify_mesh(input_dict: dict, kriged_mesh: pd.DataFrame, settings_dict: dict) -> pd.DataFrame: