Geostatistics

Geostatistics#

Mesh Cropping and Coordinate Transformations#

echopop.geostatistics.hull_crop(transects: DataFrame, mesh: DataFrame, num_nearest_transects: int = 3, mesh_buffer_distance: float = 2.5, projection: str = 'epsg:4326', coordinate_names: tuple[str, str] = ('longitude', 'latitude')) → DataFrame#

Crop the kriging mesh using convex hull polygons generated from survey transects.

This function creates a survey boundary by generating convex hulls around each transect and its nearest neighbors, then filters the mesh to include only cells within the buffered survey area. This approach provides a more flexible alternative to region-based cropping.

Parameters:
transectspandas.DataFrame

Georeferenced survey transect data used for defining the spatial extent for the kriging mesh grid. Must contain columns: 'longitude', 'latitude', 'transect_num'.

meshpandas.DataFrame

Complete kriging mesh DataFrame that is subsequently cropped. Must contain columns: 'longitude', 'latitude'.

num_nearest_transectsint, default=3

The number of nearest-neighbor transects used for defining the local extent around each transect. These convex hulls are then combined to generate the full survey extent hull. Higher values create more inclusive boundaries.

mesh_buffer_distancefloat, default=2.5

Buffer distance in nautical miles applied to the survey polygon before filtering mesh cells. This ensures adequate coverage around the survey boundary.

projectionstr, default=’epsg:4326’

EPSG projection code for the input coordinate system. Default is WGS84.

coordinate_namesTuple[str, str], default=(“longitude”, “latitude”)

Names of the coordinate columns when using DataFrames. Expected format: (x_col, y_col).

Returns:
pandas.DataFrame

Cropped mesh DataFrame containing only cells within the buffered survey extent. The 'geometry' column is removed from the output.

Notes

The function performs the following steps:

  1. Converts the mesh DataFrame to a geopandas.GeoDataFrame with point geometries

  2. Transforms coordinates from WGS84 to UTM for accurate distance calculations

  3. Generates survey extent polygon using echopop.geostatistics.transect_extent()

  4. Applies buffer distance (converted from nautical miles to meters)

  5. Filters mesh cells to those within the buffered polygon

  6. Returns the filtered mesh without geometry column

The UTM transformation ensures accurate distance calculations for the convex hull generation and buffering operations. The buffer distance helps ensure adequate mesh coverage around the survey boundary.

This method is particularly useful for irregularly shaped survey areas where region-based cropping may be too restrictive or complex.

Examples

>>> cropped_mesh = hull_crop(
...     transects, mesh,
...     num_nearest_transects=5,
...     mesh_buffer_distance=3.0
... )
>>> print(f"Original mesh size: {len(mesh)}")
>>> print(f"Cropped mesh size: {len(cropped_mesh)}")
echopop.geostatistics.transect_coordinate_centroid(spatial_grouped: GeoSeries)#

Calculate the centroid of a given spatial group.

This function computes the geometric centroid of a collection of spatial points, which is useful for determining the center point of transect lines or other spatial groupings.

Parameters:
spatial_grouped: geopandas.GeoSeries

A geopandas.GeoSeries comprising coordinates (i.e. points). Each element should be a shapely.Point geometry representing spatial locations.

Returns:
Point

A shapely.Point object representing the centroid of all input coordinates.

Notes

The function uses geopandas.GeoSeries.union_all() to combine all geometries before calculating the centroid, which ensures proper handling of the spatial reference system.

Examples

>>> import geopandas as gpd
>>> from shapely.geometry import Point
>>> points = gpd.GeoSeries([Point(0, 0), Point(1, 1), Point(2, 0)])
>>> centroid = transect_coordinate_centroid(points)
>>> print(f"Centroid: ({centroid.x:.1f}, {centroid.y:.1f})")
Centroid: (1.0, 0.3)
echopop.geostatistics.transform_coordinates(data: DataFrame, x_offset: float = 0.0, y_offset: float = 0.0, coordinate_names: tuple[str, str] = ('longitude', 'latitude'), reference: DataFrame | None = None, delta_x: float | None = None, delta_y: float | None = None) → tuple[DataFrame, float | None, float | None]#

Transform the x- and y-coordinates of a georeferenced dataset.

Parameters:
datapandas.DataFrame

DataFrame with coordinates

x_offsetfloat, default=0.

Offset to apply to the x-coordinates that corresponds to coordinate_names[0]

y_offsetfloat, default=0.

Offset to apply to the y-coordinates that corresponds to coordinate_names[0]

coordinate_namesTuple[str, str], default=(“longitude”, “latitude”)

Names of the coordinate columns when using DataFrames. Expected format: (x_col, y_col).

referencepandas.DataFrame, optional

Reference DataFrame with x and y coordinates for interpolation that is used as an additional offset to the x-axis.

delta_xfloat, optional

Total x-axis distance used for standardizing coordinates. Will use the full range of the original x-axis if not provided.

delta_yfloat, optional

Total y-axis distance used for standardizing coordinates. Will use the full range of the original y-axis if not provided.

Returns:
pandas.DataFrame

DataFrame with the new transformed coordinates ‘x’ and ‘y’.

float or None

Distance of the pre-transformed x-axis coordinates that can be used to transform other georeferenced datasets (assuming shared projections).

float or None

Distance of the pre-transformed y-axis coordinates that can be used to transform other georeferenced datasets (assuming shared projections).

echopop.geostatistics.transect_extent(transects: DataFrame, projection: str, num_nearest_transects: int, **kwargs)#

Compute the spatial extent of survey transects using convex hull generation.

This function creates a polygon representing the spatial extent of survey transects by generating convex hulls around each transect and its nearest neighbors, then ‘unioning’ all hulls to create the overall survey boundary.

Parameters:
transectspandas.DataFrame

Dataframe containing survey transect data with columns:

  • 'longitude': Longitude coordinates

  • 'latitude': Latitude coordinates

  • 'transect_num': Transect identifier numbers

projectionstr

EPSG projection code string (e.g., 'epsg:4326' for WGS84)

num_nearest_transectsint

Number of nearest neighbor transects to include when generating the convex hull around each transect

Returns:
shapely.Point or shapely.LineString or shapely.Polygon

A shapely.Point, shapely.LineString, or shapely.Polygon object that matches the geometry-type of the input transect data. This object represents the union of all transect convex hulls, defining the overall spatial extent of the survey area.

Notes

The function performs the following steps:

  1. Converts the DataFrame to a geopandas.GeoDataFrame with point geometries

  2. Transforms coordinates from WGS84 to UTM for accurate distance calculations

  3. Calculates centroids for each transect

  4. For each transect, finds the nearest neighbor transects

  5. Generates convex hulls around each transect and its neighbors

  6. Returns the union of all convex hulls

The resulting polygon can be used for spatial filtering, mesh cropping, or defining survey boundaries for analysis.

Examples

>>> import pandas as pd
>>> transect_data = pd.DataFrame({
...     'longitude': [-125.0, -125.1, -125.2],
...     'latitude': [48.0, 48.1, 48.2],
...     'transect_num': [1, 2, 3]
... })
>>> extent = transect_extent(transect_data, 'epsg:4326', 2)
>>> print(f"Extent type: {type(extent)}")

Extent type: shapely.Polygon

echopop.geostatistics.utm_string_generator(longitude: float, latitude: float)#

Generate UTM EPSG projection string from longitude/latitude coordinates.

This function converts WGS84 coordinates to the appropriate UTM zone EPSG code string, automatically determining the correct UTM zone based on longitude and the hemisphere (north/south) based on latitude.

Parameters:
longitudefloat

Longitude coordinate in decimal degrees (WGS84)

latitudefloat

Latitude coordinate in decimal degrees (WGS84)

Returns:
str

EPSG code string for the appropriate UTM zone (e.g., ‘32610’ for UTM zone 10N)

Notes

UTM zones are numbered from 1 to 60, with each zone spanning 6 degrees of longitude. The zone number is calculated as: floor((longitude + 180) / 6) + 1

EPSG codes follow the pattern:

  • 326XX for northern hemisphere (where XX is the zero-padded zone number)

  • 327XX for southern hemisphere (where XX is the zero-padded zone number)

Examples

>>> utm_code = utm_string_generator(-125.0, 48.0)
>>> print(f"UTM EPSG code: {utm_code}")
UTM EPSG code: 32610

>>> utm_code = utm_string_generator(-125.0, -48.0)
>>> print(f"UTM EPSG code: {utm_code}")
UTM EPSG code: 32710
echopop.geostatistics.wgs84_to_utm(geodataframe: GeoDataFrame)#

Transform a geopandas.GeoDataFrame from WGS84 to the appropriate UTM coordinate system.

This function automatically determines the correct UTM zone based on the median longitude and latitude of the geopandas.GeoDataFrame and transforms the coordinate reference system (CRS) in place. The transformation improves accuracy for distance and area calculations.

Parameters:
geodataframegeopandas.GeoDataFrame

GeoDataFrame containing spatial data with WGS84 coordinates. Must contain columns with 'lat' and 'long' in their names (case-insensitive).

Returns:
None

The function modifies the geopandas.GeoDataFrame in place by changing its CRS.

Notes

The function performs the following steps: 1. Automatically detects longitude and latitude columns by name 2. Calculates the median coordinates to determine the appropriate UTM zone 3. Generates the UTM EPSG code using utm_string_generator() 4. Transforms the GeoDataFrame to the new CRS in place

The transformation uses the median coordinates to ensure the UTM zone is appropriate for the entire dataset, which is particularly important for datasets spanning multiple UTM zones.

Examples

>>> import geopandas as gpd
>>> from shapely.geometry import Point
>>> df = gpd.GeoDataFrame({
...     'longitude': [-125.0, -125.1],
...     'latitude': [48.0, 48.1],
...     'geometry': [Point(-125.0, 48.0), Point(-125.1, 48.1)]
... }, crs='epsg:4326')
>>> print(f"Original CRS: {df.crs}")
Original CRS: EPSG:4326
>>> wgs84_to_utm(df)
>>> print(f"Transformed CRS: {df.crs}")
Transformed CRS: EPSG:32610
echopop.geostatistics.projection.reproject_dataset(data: DataFrame, crs_out: str, coordinate_names: tuple[str, str] = ('longitude', 'latitude'), projection: str = 'epsg:4326') → DataFrame#

Transform coordinates using a new projection via GeoPandas.

This function converts coordinates from one coordinate reference system (CRS) to another using GeoPandas for accurate cartographic projections. It creates a temporary geopandas.GeoDataFrame to perform the transformation and returns the results as a pandas.DataFrame with new 'x' and 'y' columns containing the projected coordinates.

Parameters:
data_dfpandas.DataFrame

DataFrame containing coordinate data to be reprojected.

crs_outstr

Target Coordinate Reference System (CRS) string (e.g., 'epsg:32610' for UTM Zone 10N).

coordinate_namesTuple[str, str], default=(“longitude”, “latitude”)

Names of the coordinate columns in the input DataFrame. Expected format: (x_col, y_col).

projectionstr, default=’epsg:4326’

Input Coordinate Reference System (CRS) string representing the original projection of the coordinate data (default is WGS84 geographic coordinates).

Returns:
pandas.DataFrame

DataFrame with the original data plus new 'x' and 'y' columns containing the reprojected coordinates in the target CRS.

Notes

This function leverages GeoPandas for accurate coordinate transformations, which is more reliable than simple mathematical conversions for cartographic projections. The geometry column created during processing is automatically dropped from the output pandas.DataFrame.

Examples

>>> # Project from WGS84 to UTM Zone 10N
>>> df_projected = reproject_dataset(
...     data=survey_data,
...     crs_out='epsg:32610',
...     coordinate_names=('longitude', 'latitude')
... )
>>> print(df_projected[['x', 'y']].head())

>>> # Project from UTM back to WGS84
>>> df_geo = reproject_dataset(
...     data=utm_data,
...     crs_out='epsg:4326',
...     coordinate_names=('x', 'y'),
...     projection='epsg:32610'
... )

(Semi)Variogram and Covariance Models#

class echopop.geostatistics.Variogram(lag_resolution: float, n_lags: int, coordinate_names: tuple[str, str] = ('x', 'y'))#

Class for calculating and fitting variograms to quantify spatial variance.

A variogram (or semivariogram) is a fundamental tool in geostatistics that describes the spatial correlation structure of a variable as a function of distance. The empirical semivariogram \(\\hat{\\gamma}(h)\) is defined as:

\[\begin{split}\\hat{\\gamma}(h) = \\frac{1}{2N(h)} \\sum\\limits_{i<j:\\; h_{ij}\\approx h} \\bigl(z(\\mathbf{x}_j)-z(\\mathbf{x}_i)\\bigr)^2\end{split}\]

here \(z(\\mathbf{x}_i)\) and \(z(\\mathbf{x}_j)\) denote the observed (realized) values of the variable at sample locations \(\\mathbf{x}_i\) and \(\\mathbf{x}_j\), respectively, and \(N(h)\) is the number of pairs separated by distance h.

The class implements a standardized semivariogram approach where the semivariance is normalized by the product of head and tail standard deviations:

\[\begin{split}\\hat{\\gamma}_{\\text{std}}(h) = \\frac{1}{2N(h)} \\sum_{i<j:\\; h_{ij}\\approx h} \\frac{[z(\\mathbf{x}_j) - z(\\mathbf{x}_i)]^2}{\\sigma(\\mathbf{x}_j) \\cdot \\sigma(\\mathbf{x}_i)}\end{split}\]

where \(\\sigma(\\mathbf{x}_i)\) and \(\\sigma(\\mathbf{x}_j)\) are local standard deviations at \(\\mathbf{x}_i\) and \(\\mathbf{x}_j\). This standardization helps account for local variance heterogeneity and provides more robust variogram estimates for irregularly distributed data [1].

Parameters:
lag_resolutionfloat

The distance interval represented by each lag. Must be positive and determines the spatial resolution of the variogram analysis.

n_lagsint

The number of lags used for computing the semivariogram. Must be positive. Typical values range from 10-30 depending on data density and spatial extent.

coordinate_namesTuple[str, str], default=(‘x’, ‘y’)

Column names for the spatial coordinates in input DataFrames. Format: (horizontal_axis, vertical_axis).

Attributes:
gammanumpy.ndarray[float] or None

The computed semivariance values for each lag distance.

lagsnumpy.ndarray[float] or None

The lag distances used in the variogram analysis.

lag_countsnumpy.ndarray[int] or None

The number of data pairs contributing to each lag estimate.

lag_covariancefloat or None

The mean lag covariance between head and tail points.

variogram_params_optimizeddict or None

Optimized parameters from theoretical variogram model fitting.

variogram_params_initialdict or None

Initial parameters used for theoretical variogram model fitting.

variogram_fit_initialfloat or None

Mean absolute deviation of initial model fit.

variogram_fit_optimizedfloat or None

Mean absolute deviation of optimized model fit.

Methods

calculate_empirical_variogram(data, variable, azimuth_filter=True,

azimuth_angle_threshold=180.0, force_lag_zero=True)

Compute the empirical variogram from transect data.

fit_variogram_model(model, model_parameters, optimizer_kwargs={})

Fit a theoretical variogram model to the empirical variogram using weighted least squares.

Notes

The class implements several key applications for fisheries acoustic survey data:

  1. Standardized Semivariogram: Uses local variance normalization to handle heteroscedastic data common in biological surveys.

  2. Azimuth Filtering: Supports directional variogram analysis for anisotropic spatial patterns often observed in marine ecosystems.

  3. Robust Lag Estimation: Uses weighted binning approaches that account for irregular sampling patterns typical in acoustic transect surveys.

The theoretical variogram models available include single-family models exponential, Gaussian, spherical, Bessel functions) and composite models that can capture complex spatial structures with both short-range correlation and periodic patterns [2], [3].

References

[1]

Cressie, N.A.C. (1993). Statistics for Spatial Data. John Wiley & Sons.

[2]

Chilès, J.P., and Delfiner, P. (2012). Geostatistics: Modeling Spatial Uncertainty. 2nd ed. John Wiley & Sons.

[3]

Rivoirard, J., Simmonds, J., Foote, K.G., Fernandes, P., and Bez, N. (2000). Geostatistics for Estimating Fish Abundance. Blackwell Science.

Examples

Basic usage for computing an empirical variogram:

>>> import pandas as pd
>>> import numpy as np
>>> from echopop.nwfsc_feat.variogram import Variogram
>>>
>>> # Create sample spatial data
>>> data = pd.DataFrame({
...     'x': np.random.uniform(0, 100, 50),
...     'y': np.random.uniform(0, 100, 50),
...     'biomass_density': np.random.exponential(2, 50)
... })
>>>
>>> # Initialize variogram
>>> vario = Variogram(lag_resolution=5.0, n_lags=15)
>>>
>>> # Compute empirical variogram
>>> vario.calculate_empirical_variogram(data, 'biomass_density')
>>>
>>> # Fit theoretical model
>>> from lmfit import Parameters
>>> params = Parameters()
>>> params.add('sill', value=2.0, min=0.1)
>>> params.add('nugget', value=0.1, min=0.0)
>>> params.add('correlation_range', value=20.0, min=1.0)
>>>
>>> best_fit = vario.fit_variogram_model('exponential', params)
calculate_empirical_variogram(data: DataFrame, variable: str, azimuth_filter: bool = True, azimuth_angle_threshold: float = 180.0, force_lag_zero: bool = True) → None#

Compute the empirical variogram from transect data.

Calculates the standardized semivariogram using the method described in Rivoirard et al. (2000) [1], which is particularly suitable for fisheries acoustic data with irregular sampling patterns and heterogeneous variance.

Parameters:
datapandas.DataFrame

DataFrame containing georeferenced coordinates and the target variable. Must include columns specified in coordinate_names and variable.

variablestr

Column name of the variable for variogram computation (e.g., 'biomass_density').

azimuth_filterbool, default=True

If True, computes azimuth angles between point pairs to enable directional filtering. Useful for detecting spatial anisotropy.

azimuth_angle_thresholdfloat, default=180.0

Maximum azimuth angle deviation (in degrees) for including point pairs. Values less than 180° enable directional variogram analysis.

force_lag_zerobool, default=True

If True, forces the nugget effect to zero by prepending lag 0 with \(\\gamma(0) = 0\). This assumes no measurement error or micro-scale variation.

Notes

The method uses a triangular mask to ensure each point pair is counted only once, avoiding redundant calculations in the distance matrix. For acoustic survey data, typical azimuth_angle_threshold values of 45-90° can help identify along-track vs. across-track spatial patterns.

The lag counts provide important information about the reliability of each lag estimate - lags with very few pairs should be interpreted cautiously.

References

[1]

Rivoirard, J., Simmonds, J., Foote, K.G., Fernandes, P., and Bez, N. (2000). Geostatistics for Estimating Fish Abundance. Blackwell Science.

fit_variogram_model(model: str | list[str], model_parameters: Parameters, optimizer_kwargs: dict[str, Any] | None = None) → dict[str, Any]#

Fit a theoretical variogram model to the empirical variogram.

Uses weighted least squares optimization. Fits parametric models of the form:

\[\begin{split}\\gamma(h) = C_0 + C_1 \\cdot \\mathscr{f}(h; \\theta)\end{split}\]

where \(C_0\) is the nugget (the value at zero lag), \(C_1\) is the partial sill (\(C - C_0\)), and \(\\mathscr{f}(h; \\theta)\) is a correlation function that depends on lag \(h\) and model-specific parameters \(\\theta\) (e.g., range, sill, power)

The optimization minimizes the weighted residual sum of squares:

\[\begin{split}\\min_{\\theta} \\sum_{b=1}^{n} \\hat{w}_b \\left[ \\gamma_\\text{empirical}(h_b) - \\gamma_\\text{model}(h_b; \\theta) \\right]^2\end{split}\]

where \(n\) is the number of lag bins and \(b\) indexes those lag bins (\(b=1, \\dots, n\)). The weights \(\\hat{w}_b\) are associated with each lag bins and are normalized where:

\[\begin{split}\\hat{w}_b = \\frac{N(h_b)}{\\sum\\limits_{b=1}^{n} N(h_b)}\end{split}\]
Parameters:
modelstr or List[str]

Theoretical variogram model name(s). Single string for basic models (e.g., 'exponential', 'gaussian', 'spherical'). List of two strings for composite models (e.g., ['bessel', 'exponential']).

model_parameterslmfit.Parameters

Initial parameter values and constraints for optimization. Required parameters depend on the chosen model, but the formatting should follow lmfit.parameter.Parameters specifications. See echopop.geostatistics.compute_variogram() for details on available models and their respective parameter requirements.

optimizer_kwargsdict, default={}

Additional arguments passed to lmfit.minimizer.minimize(). Common options include 'max_nfev' (maximum function evaluations) and solver-specific parameters.

Returns:
dict

Dictionary of optimized parameter values with parameter names as keys.

Notes

The method uses Trust Region Reflective algorithm for bounded optimization, which handles parameter constraints robustly [1].

References

[1]

Branch, M.A., Coleman, T.F., and Li, Y. (1999). A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems. SIAM Journal on Scientific Computing, 21(1), 1-23.

echopop.geostatistics.compute_variogram(distance_lags: ndarray[float], variogram_parameters: dict[str, float] | None = None, model: str | list[str] | None = None, **kwargs)#

Compute the theoretical semivariogram.

Parameters:
distance_lagsnumpy.ndarray[float]

An array of lag distances.

variogram_parametersOptional[Dict[str, float]]

An optional dictionary that can contain values for variogram model parameters (see the below table associated with the argument model). Alternatively, these parameters can be entered directly and are contained within kwargs. Possible parameters include:

  • sill (Sill): The asymptotic value as lags approach infinity.

  • nugget (Nugget): The semivariogram y-intercept that corresponds to variability at lag distances shorter than the lag resolution.

  • correlation_range (Correlation length scale/range): The ascending rate for the semivariogram.

  • hole_effect_range (Hole effect range): The (normalized) length scale/range that ‘holes’ are observed, which represent ‘null’ (or very small) points compared to their neighboring lags.

  • decay_power (Decay term exponent): An exponential term that is used in certain generalized exponential (or related) semivariogram models that modulates the ascending rate for a semivariogram.

  • enhance_semivariance (Semivariance enhancement): A boolean term that determines whether the correlation decay in certain cosine-related variogram models are enhanced (or not) at further lag distances.

modelOptional[Union[str, list]]

A string or list of model names. A single name represents a single family model. Two inputs represent the desired composite model (e.g. the composite J-Bessel and exponential model). Available variogram models and their respective arguments include (alongside distance_lags):

variogram() model

Input

Parameters

cubic()

'cubic'

  • sill

  • nugget

  • correlation_range

exponential()

'exponential'

  • sill

  • nugget

  • correlation_range

gaussian()

'gaussian'

  • sill

  • nugget

  • correlation_range

jbessel()

'jbessel'

  • sill

  • nugget

  • hole_effect_range

kbessel()

'kbessel'

  • sill

  • nugget

  • hole_effect_range

linear()

'linear'

  • nugget

  • sill

matern()

'matern'

  • sill

  • nugget

  • correlation_range

  • smoothness_parameter

nugget()

'nugget'

  • nugget

  • sill

pentaspherical()

'pentaspherical'

  • sill

  • nugget

  • correlation_range

power()

'power'

  • nugget

  • sill

  • power_exponent

quadratic()

'quadratic'

  • sill

  • nugget

  • correlation_range

  • shape_parameter

sinc()

'sinc'

  • sill

  • nugget

  • hole_effect_range

spherical()

'spherical'

  • sill

  • nugget

  • correlation_range

bessel_exponential()

['bessel', 'exponential']

  • sill

  • nugget

  • correlation_range

  • decay_power

  • hole_effect_range

bessel_gaussian()

['bessel', 'gaussian']

  • sill

  • nugget

  • correlation_range

  • decay_power

  • hole_effect_range

cosine_exponential()

['cosine', 'exponential']

  • sill

  • nugget

  • correlation_range

  • hole_effect_range

  • enhance_semivariance

cosine_gaussian()

['cosine', 'gaussian']

  • sill

  • nugget

  • correlation_range

  • hole_effect_range

exponential_linear()

['exponential', 'linear']

  • sill

  • nugget

  • correlation_range

  • hole_effect_range

  • decay_power

gaussian_linear()

['gaussian', 'linear']

  • sill

  • nugget

  • correlation_range

  • hole_effect_range
Returns:
variogramnp.ndarray

An array containing the (normalized) semivariance for each lag bin.

echopop.geostatistics.fit_variogram(lags: ndarray[float], lag_counts: ndarray[int], gamma: ndarray[float], model_parameters: Parameters, model: str | list[str] | None = None, optimizer_kwargs: dict[str, Any] | None = None) → tuple[dict[str, Any], float, float]#

Fit theoretical variogram models to empirical semivariogram data using weighted least squares.

This function performs non-linear optimization to find the best-fitting parameters for theoretical variogram models. The optimization uses weighted least squares where weights are proportional to the lag counts, giving more influence to lags with more data pairs.

Parameters:
lagsnumpy.ndarray[float]

Array of lag distances from empirical variogram computation.

lag_countsnumpy.ndarray[int]

Number of data point pairs contributing to each lag estimate. Used as optimization weights.

gammanumpy.ndarray[float]

Empirical semivariogram values (standardized semivariance) at each lag.

modelstr or List[str]

Theoretical variogram model specification. Single string for basic models (e.g., 'exponential', 'gaussian', 'spherical', 'jbessel', 'linear'). List of two strings for composite models (e.g., ['bessel', 'exponential'], ['bessel', 'gaussian'], ['cosine', 'exponential']).

model_parameterslmfit..parameters.Parameters

Parameter object containing initial values, bounds, and constraints for optimization. Required parameters depend on the selected model.

optimizer_kwargsdict, default={}

Additional keyword arguments passed to lmfit.minimizer.minimize(). Common options include 'max_nfev' for maximum function evaluations and solver-specific parameters.

Returns:
Tuple[dict, dict, float]

  • Optimized parameter values as dictionary

  • Initial parameter values as dictionary

  • Mean absolute deviation of optimized fit

Notes

The optimization minimizes the weighted residual sum of squares:

\[\begin{split}\\min_{\\theta} \\sum_{b=1}^{n} \\hat{w}_b \\left[ \\gamma_\\text{empirical}(h_b) - \\gamma_\\text{model}(h_b; \\theta) \\right]^2\end{split}\]

where \(n\) is the number of lag bins and \(b\) indexes those lag bins (\(b=1, \\dots, n\)). The weights \(\\hat{w}_b\) are associated with each lag bins and are normalized where:

\[\begin{split}\\hat{w}_b = \\frac{N(h_b)}{\\sum\\limits_{b=1}^{n} N(h_b)}\end{split}\]

The function uses Trust Region Reflective algorithm (default in lmfit.minimizer.minimize()) which handles parameter bounds robustly and is suitable for the non-linear nature of variogram models.

References

[1]

Cressie, N. (1993). Statistics for Spatial Data. Wiley.

[2]

Newville, M., et al. (2014). LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python. Zenodo.

Kriging and Spatial Interpolation#

class echopop.geostatistics.Kriging(mesh: DataFrame, kriging_params: dict[str, Any], variogram_params: dict[str, Any], coordinate_names: tuple[str, str] = ('x', 'y'))#

Class for kriging spatially distributed data.

Class for performing ordinary kriging to predict population values and other metrics at un-sampled locations defined by a mesh grid.

Kriging is a geostatistical interpolation technique that provides the Best Linear Unbiased Predictor (BLUP) for spatial data. The method uses the spatial correlation structure, characterized by a semivariogram model, to optimally weight nearby observations when predicting values at unsampled locations.

For ordinary kriging, the prediction at location x₀ is given by:

\[\begin{split}\\mathbf{z}^*(\\mathbf{u}) = \\sum_{b=1}^n \\lambda_b(\\mathbf{u}) z(\\mathbf{u}_b)\end{split}\]

where \(n\) is the number of observed locations, \(z(\\mathbf{u}_b)\) are the known values of a spatially varying variable at locations \(\\mathbf{u}_b\), and \(\\lambda_b(\\mathbf{u})\) are the kriging weights assigned to all known locations \(\\mathbf{u}\). This kriging approach is subject to the unbiasedness constraint where:

\[\begin{split}\\sum\\limits_{b=1}^n \\lambda_b = 1\end{split}\]

The kriging weights \(\\lambda_b\) are obtained by solving the kriging system:

\[\begin{split}\\mathbf{\\Gamma} \\hat{\\mathbf{\\lambda}} = \\hat{\\mathbf{\\gamma}}_\\mathbf{u}\end{split}\]

where:

\[\begin{split}\\underbrace{ \\begin{bmatrix} \\gamma_{1,1} & \\cdots & \\gamma_{1,n} & 1 \\\\ \\vdots & \\ddots & \\vdots & \\vdots \\\\ \\gamma_{n,1} & \\cdots & \\gamma_{n,n} & 1 \\\\ 1 & \\cdots & 1 & 0 \\end{bmatrix} }_{\\mathbf{\\Gamma}} \\underbrace{ \\begin{bmatrix} \\lambda_1 \\\\ \\vdots \\\\ \\lambda_n \\\\ \\mu \\end{bmatrix} }_{\\hat{\\mathbf{\\lambda}}} = \\underbrace{ \\begin{bmatrix} \\gamma(\\mathbf{u}_1 - \\mathbf{u}) \\\\ \\vdots \\\\ \\gamma(\\mathbf{u}_n - \\mathbf{u}) \\\\ 1 \\end{bmatrix} }_{\\mathbf{\\hat{\\gamma}}_\\mathbf{u}}\\\\\end{split}\]

Here, \(\\gamma_{i,j}\) is the semivariance between locations \(i\) and \(j\), \(\\gamma_{i,\\mathbf{u}} = \\gamma(\\mathbf{u}_i - \\mathbf{u})\) is the semivariance between data point \(i\) and target location \(\\mathbf{u}\), and \(\\mu\) is the Lagrange multiplier enforcing the unbiasedness constraint [1], [2].

Parameters:
meshpandas.DataFrame

DataFrame containing the mesh grid used for interpolating the values from data. This DataFrame must contain coordinate columns as specified in coordinate_names.

Optional columns include:

  • area (float): Cell areas in square nautical miles for projection calculations

  • fraction (float): Fraction of cell area if using default cell areas

coordinate_namesTuple[str, str], default=(“x”, “y”)

Names of the coordinate columns that are shared between input data and mesh. Format: (horizontal_coordinate, vertical_coordinate).

kriging_paramsDict[str, Any]

Dictionary containing kriging parameters required for interpolation:

  • aspect_ratio (float): Ratio of minor to major axis correlation ranges for anisotropy handling. Values near 1 indicate isotropy, smaller values indicate directional elongation.

  • k_min (int): Minimum number of nearest neighbors for kriging (typically 3-8).

  • k_max (int): Maximum number of nearest neighbors for kriging (typically 8-20).

  • search_radius (float): Maximum distance for neighbor search in coordinate units.

variogram_paramsDict[str, Any]

Dictionary containing variogram model parameters:

  • model (str or List[str]): Variogram model (e.g., 'exponential', 'gaussian', 'spherical', or composite models like ['bessel', 'exponential']).

  • nugget (float): Nugget effect representing micro-scale variability.

  • sill (float): Total variance (nugget + partial sill).

  • correlation_range (float): Correlation length scale.

  • Additional parameters depending on model choice (e.g., hole_effect_range, decay_power).

Attributes:
meshpandas.DataFrame

Original mesh grid for interpolation.

coordinate_namesTuple[str, str]

Coordinate column names.

kriging_paramsDict[str, Any]

Kriging parameters dictionary.

variogram_paramsDict[str, Any]

Variogram model parameters dictionary.

mesh_croppedpandas.DataFrame or None

Cropped mesh grid to prevent extrapolation beyond survey boundaries.

survey_cvfloat or None

Overall survey coefficient of variation after kriging.

Methods

crop_mesh([crop_function, coordinate_names])

Crop the mesh grid to prevent extrapolation beyond survey boundaries.

krige(transects, variable[, extrapolate, ...])

Perform ordinary kriging interpolation and project results onto the mesh grid.

register_search_strategy(name, strategy)

Register a custom search strategy function.

list_search_strategies()

List all available search strategies.

Notes

Key Features:

1. Adaptive Search Strategy: Uses k-nearest neighbor search with fallback extrapolation for sparse data regions.

2. Numerical Stability: Employs truncated SVD for matrix inversion to handle ill-conditioned covariance matrices.

3. Anisotropy Support: Handles directional correlation through aspect ratio parameterization.

4. Mesh Cropping: Prevents extrapolation beyond survey boundaries using convex hull or custom boundary functions.

Typical Parameter Ranges:

  • k_min: 3-8 (minimum for stable estimates)

  • k_max: 8-20 (balance between locality and stability)

  • search_radius: 2-5× correlation range

  • aspect_ratio: 0.1-1.0 (lower values = stronger anisotropy)

Performance Considerations:

  • Computational complexity: \(\\mathcal{O}(n \\times m \\times k^3)\) where \(n\) = mesh points, \(m\) = data points, \(k\) = neighbors

  • Memory usage scales with mesh size and neighbor count

  • Large search radii increase computational cost but improve spatial continuity

References

[1]

Cressie, N.A.C. (1993). Statistics for Spatial Data. John Wiley & Sons.

[2]

Chilès, J.P., and Delfiner, P. (2012). Geostatistics: Modeling Spatial Uncertainty. 2nd ed. John Wiley & Sons.

Examples

Basic kriging workflow for fisheries acoustic data:

>>> import pandas as pd
>>> import numpy as np
>>> from echopop.nwfsc_feat.kriging import Kriging
>>>
>>> # Create sample transect data
>>> transects = pd.DataFrame({
...     'longitude': np.random.uniform(-125, -120, 100),
...     'latitude': np.random.uniform(40, 45, 100),
...     'biomass_density': np.random.exponential(50, 100)
... })
>>>
>>> # Create kriging mesh
>>> lon_grid, lat_grid = np.meshgrid(
...     np.linspace(-125, -120, 20),
...     np.linspace(40, 45, 20)
... )
>>> mesh = pd.DataFrame({
...     'longitude': lon_grid.flatten(),
...     'latitude': lat_grid.flatten(),
...     'area': np.full(400, 1.0)  # 1 nmi² cells
... })
>>>
>>> # Define kriging parameters
>>> kriging_params = {
...     'k_min': 4,
...     'k_max': 12,
...     'search_radius': 10.0,
...     'aspect_ratio': 0.8
... }
>>>
>>> # Define variogram parameters (from fitted model)
>>> variogram_params = {
...     'model': 'exponential',
...     'nugget': 0.1,
...     'sill': 2.5,
...     'correlation_range': 8.0
... }
>>>
>>> # Initialize kriging object
>>> krig = Kriging(
...     mesh=mesh,
...     coordinate_names=('longitude', 'latitude'),
...     kriging_params=kriging_params,
...     variogram_params=variogram_params
... )
>>>
>>> # Perform kriging interpolation
>>> results = krig.krige(
...     transects=transects,
...     variable='biomass_density',
...     extrapolate=False  # Requires mesh cropping first
... )
crop_mesh(crop_function: ~collections.abc.Callable = <function hull_crop>, coordinate_names: tuple[str, str] = ('longitude', 'latitude'), **kwargs) → None#

Crop the mesh grid to prevent extrapolation beyond survey boundaries.

Mesh cropping is essential for avoiding unreliable kriging predictions in regions with sparse or no data coverage. The default method uses convex hull boundaries around survey transects, but custom boundary functions can be specified for more complex survey geometries.

Parameters:
crop_functionCallable, default=hull_crop

Function that defines the survey boundary for mesh subsetting. The default hull_crop creates convex hull polygons around transect data with optional buffering. Custom functions must accept mesh as a keyword argument. See echopop.geostatistics.hull_crop() for more details.

coordinate_namesTuple[str, str], default=(“longitude”, “latitude”)

Column names for spatial coordinates in the cropping function. Uses geographic coordinates by default for compatibility with hull_crop.

**kwargs

Additional arguments passed to the cropping function. For hull_crop:

  • 'transects' (pandas.DataFrame): Survey transect data for boundary definition

  • 'num_nearest_transects'(int): Number of neighbors for local hull creation

  • 'mesh_buffer_distance' (float): Buffer distance in nautical miles

  • 'projection' (str): EPSG code for coordinate system (e.g., 'epsg:4326')

Returns:
None

Updates the mesh_cropped attribute with the subsetted mesh.

Notes

The cropping process typically involves:

  1. Boundary Definition: Creating polygons around survey transects

  2. Buffer Application: Expanding boundaries to include nearby areas

  3. Mesh Intersection: Selecting mesh points within boundaries

  4. Area Adjustment: Updating cell areas for boundary cells if needed

Boundary Methods:

  • Convex Hull: Simple, conservative boundary (default)

  • Alpha Shapes: More flexible boundaries for complex geometries

  • Custom Polygons: User-defined survey strata or management areas

Examples

Crop mesh using convex hull with 2.5 nmi buffer:

>>> krig.crop_mesh(
...     transect_df=survey_data,
...     num_nearest_transects=5,
...     mesh_buffer_distance=2.5,
...     projection='epsg:4326'
... )

Using a custom cropping function:

>>> def custom_crop(mesh, boundary_polygon, **kwargs):
...     # Custom boundary logic
...     return mesh[mesh_within_polygon(mesh, boundary_polygon)]
>>>
>>> krig.crop_mesh(
...     crop_function=custom_crop,
...     boundary_polygon=survey_stratum
... )
krige(transects: DataFrame, variable: str, extrapolate: bool = True, default_mesh_cell_area: float | None = None, adaptive_search_strategy: str = 'uniform', custom_search_kwargs: dict[str, Any] | None = None) → DataFrame#

Perform ordinary kriging interpolation and project results onto the mesh grid.

This method implements the complete kriging workflow: neighbor selection, covariance matrix construction, weight calculation, prediction, and variance estimation. The results are projected onto the mesh grid with proper area weighting for survey-level biomass estimates.

Parameters:
transectspandas.DataFrame

Georeferenced survey data containing coordinates and the target variable. Must include columns specified in coordinate_names and variable.

variablestr

Column name of the variable to interpolate (e.g., 'biomass_density').

extrapolatebool, default=True

If True, uses the full mesh grid (may extrapolate beyond data coverage). If False, uses the cropped mesh (requires prior call to echopop.geostatistics.Kriging.crop_mesh()).

default_mesh_cell_areafloat, optional

Default area (nmi²) for mesh cells when 'area' column is missing from mesh. Required if mesh lacks area information and no 'fraction' column exists.

adaptive_search_strategystr, default=’uniform’

Name of the search strategy for handling sparse data regions. Built-in strategies:

  • 'uniform': Applies uniform weights to extrapolated points (default).

Use echopop.geostatistics.Kriging.register_search_strategy() to add custom strategies.

custom_search_kwargsDict[str, Any], default={}

Additional keyword arguments passed to the adaptive search strategy function. Available parameters depend on the selected strategy but may include custom weighting schemes, distance thresholds, or algorithm-specific parameters. If the custom function incorporates coordinate_names or kriging_mesh as arguments, they will be inherited from the class instance. See echopop.geostatistics.uniform_strategy() for more details on internal argument names that can be added to the custom function call.

Returns:
pandas.DataFrame

Kriged results with columns:

  • Original mesh columns (coordinates, area, etc.)

  • {variable}: Column name associated variable

  • 'kriged_variance': Prediction variance from kriging equations

  • 'sample_variance': Coefficient of variation based variance

  • 'cell_cv': Cell-level coefficient of variation

Raises:
KeyError

If required columns are missing from transects or mesh lacks area info.

AttributeError

If extrapolate=False but mesh_cropped is None.

ValueError

If adaptive_search_strategy is not a registered strategy name.

Notes

Kriging Process:

  1. Neighbor Search: Find \(k\)-nearest neighbors within search radius

  2. Covariance Matrix: Build spatial covariance structure using variogram

  3. Weight Calculation: Solve kriging system with singular value decomposition (SVD) for stability

  4. Prediction: Compute weighted estimates and prediction variance

  5. Projection: Scale results by cell areas for survey totals

Variance Components:

  • Kriged Variance: From kriging equations, measures prediction uncertainty

  • Sample Variance: CV-based measure incorporating data variability

  • Survey CV: Overall coefficient of variation for the entire survey

Quality Indicators:

  • Negative predictions are truncated to zero with warnings

  • High kriged variance indicates uncertain predictions

  • Large survey CV suggests high spatial variability or poor model fit

Search Strategy Options:

The adaptive search handles regions with insufficient neighbors:

  • Interpolation: k_min ≤ neighbors ≤ k_max within search radius

  • Extrapolation: < k_min neighbors, uses distance-weighted nearest points

  • Full Extrapolation: No neighbors within radius, uses k_min nearest

Examples

Standard kriging with extrapolation:

>>> results = krig.krige(
...     transects=survey_data,
...     variable='biomass_density',
...     extrapolate=True,
...     default_mesh_cell_area=1.0
... )
>>> print(f"Survey CV: {krig.survey_cv:.3f}")

Conservative kriging without extrapolation:

>>> krig.crop_mesh(transect_df=survey_data, mesh_buffer_distance=2.0)
>>> results = krig.krige(
...     transects=survey_data,
...     variable='biomass_density',
...     extrapolate=False
... )

Registering and using a custom strategy:

>>> def conservative_search(**params):
...     # Custom logic for sparse regions
...     return modified_indices, weights
>>>
>>> Kriging.register_search_strategy('conservative', conservative_search)
>>> results = krig.krige(
...     transects=survey_data,
...     variable='biomass_density',
...     adaptive_search_strategy='conservative',
...     custom_search_kwargs={'threshold': 0.5}
... )
classmethod list_search_strategies() → list[str]#

List all available search strategies.

classmethod register_search_strategy(name: str, strategy: Callable) → None#

Register a custom search strategy function.

Parameters:
namestr

Name identifier for the strategy. Cannot overwrite built-in strategies.

strategyCallable

Function implementing the search strategy

Raises:
ValueError

If attempting to overwrite a built-in strategy

echopop.geostatistics.uniform_search_strategy(sparse_radii: ndarray[int], valid_distances: ndarray[int], local_points: ndarray[float], k_min: int, k_max: int, search_radius: float, wr_indices: ndarray[int], oos_indices: ndarray[number], oos_weights: ndarray[float], **kwargs) → tuple[ndarray[number], ndarray[number], ndarray[number]]#

Uniform nearest neighbors search strategy for ordinary kriging.

Uniform extrapolation search strategy for finding (and weighting) k-th nearest points (relative to a reference coordinate) required for computing the lagged semivariogram in an adaptive approach.

Parameters:
sparse_radiinumpy.ndarray[int]

Indices where there are fewer than k_min nearest neighbors.

valid_distancesnumpy.ndarray[int]

The number of masked distance matrix values where extrapolation is required.

local_pointsnumpy.ndarray[float]

An array with the sorted distances (from nearest to furthest) relative to each point.

k_minint

The minimum number of nearest neighbors required for including values for kriging within the search radius.

k_maxint

The maximum number of nearest neighbors required for including values for kriging detected within the search radius.

search_radiusfloat

The adaptive search radius that identifies the \(k\)-nearest neighbors around each georeferenced value that are subsequently kriged.

wr_indicesnumpy.ndarray[int]

Indices of within-radius (WR) (i.e. < k_max) points.

oos_indicesnumpy.ndarray[numpy.number]

Template array based on the size of the data input and k_min that will contain indices where extrapolation is required where there are fewer than k_min nearest neighbors.

oos_weightsnumpy.ndarray[float]

Weights applied to extraplolated values.

Returns:
Tuple[numpy.ndarray[numpy.number], numpy.ndarray[numpy.number], numpy.ndarray[numpy.number]]

A tuple with updated values for wr_indices, oos_indices, and oos_weights via a search strategy that applies unconstrained and uniform extrapolation to out-of-sample (OOS) points.