Geostatistics

Geostatistics#

While we can infer a lot of information from the along-transect population estimates, it is sometimes desired (or even necessary) to characterize the spatial processes and interpolate over unsampled areas. This is where geostatistics can be leveraged.

Import necessary modules#

With the transect data fully procedssed, we can use the geostatistics sub-package to characterize spatial variability and use ordinary kriging to interpolate estimates. The first step requires reading in a mesh grid for the kriging interpolation. This can be done using load_mesh_data from the ingest sub-package:

from echopop.ingest import load_mesh_data

FEAT_TO_ECHOPOP_MESH_COLUMNS = {
    "centroid_latitude": "latitude",
    "centroid_longitude": "longitude",
    "fraction_cell_in_polygon": "fraction",
}


df_mesh = load_mesh_data(
    mesh_filepath=DATA_ROOT
    / "Kriging_files/Kriging_grid_files/krig_grid2_5nm_cut_centroids_2013.xlsx",
    sheet_name="krigedgrid2_5nm_forChu",
    column_name_map=FEAT_TO_ECHOPOP_MESH_COLUMNS,
)

This function has an additional argument, column_name_map, which is used to rename columns from the imported file to those compatible with Echopop. Similar to the transect and biological data, we can stratify the mesh using the join_geostrata_by_latitude function:

# INPFC (from geostrata)
df_mesh = join_geostrata_by_latitude(
    df_mesh, df_dict_geostrata["inpfc"], stratum_name="geostratum_inpfc"
)
# KS (from geostrata)
df_mesh = join_geostrata_by_latitude(df_mesh, df_dict_geostrata["ks"], stratum_name="geostratum_ks")

The next step requires transforming georeferenced values to Euclidean coordinates. The transform_coordinates function can be used to achieve this for both the transect and mesh DataFrames. For this dataset, we incorporate the 200 m isobath as an additional set of reference coordinates:

from echopop.ingest import load_isobath_data

df_isobath = load_isobath_data(
    isobath_filepath=DATA_ROOT
    / "Kriging_files/Kriging_grid_files/transformation_isobath_coordinates.xlsx",
    sheet_name="Smoothing_EasyKrig",
)

Coordinate transformation#

With this 200 m isobath read in, we can now transform the coordinates:

from echopop.geostatistics import transform_coordinates

df_nasc_no_age1_prt, delta_longitude, delta_latitude = transform_coordinates(
    data=df_nasc_all_ages,
    reference=df_isobath,
    x_offset=-124.78338,
    y_offset=45.0,
)

The arguments x_offset and y_offset act as reference coordinates used to subtract from the original coordinates. The reference argument is an optional argument that allows for a reference DataFrame to be input. In this case, we used the 200 m isobath but it can by any set of coordinates latitudinally distributed. This outputs a pandas.DataFrame with the transformed coordinates as well as the distances of the pre-transformed \(x\)- and \(y\)-axes. So our previous columns of "longitude" and "latitude" are now denoted by "x" and "y", respectively. Both sets of coordinates are retained within the DataFrame, which the first 10 rows of just these coordinate columns are represented:

display(df_nasc_no_age1_prt.filter(["longitude", "latitude", "x", "y"]).head(10))
longitude latitude x y
0 -121.143005 34.397267 -0.143604 -0.519568
1 -121.133196 34.397391 -0.141486 -0.519562
2 -121.123057 34.397435 -0.139387 -0.519560
3 -121.112871 34.397394 -0.137373 -0.519562
4 -121.102888 34.397437 -0.135307 -0.519560
5 -121.092790 34.397504 -0.133193 -0.519556
6 -121.082708 34.397426 -0.131240 -0.519560
7 -121.072621 34.397341 -0.129295 -0.519564
8 -121.062490 34.397293 -0.127299 -0.519567
9 -121.052328 34.397279 -0.125260 -0.519567

We can do the same type of transformation for the kriging mesh, but we will also use the delta_longitude and delta_latitude outputs from the transect data to further standardize the coordinates. This ensures that both the transect data and mesh share compatible coordinate systems.

df_mesh, _, _ = transform_coordinates(
    data=df_mesh,
    reference=df_isobath,
    x_offset=-124.78338,
    y_offset=45.0,
    delta_x=delta_longitude,
    delta_y=delta_latitude,
)

The Variogram#

With our transformed coordinates, we can now begin to characterize the spatial variability in our transect data. This is done by using semivariograms. In Echopop, this is done through the Variogram-class from the geostatistics sub-package:

from echopop.geostatistics import Variogram

There are three required arguments for initializing a Variogram instance:

  • lag_resolution: The distance interval represented by each lag. This is a positive-only quantity that determins the spatial resolution of the variogram analysis.

  • n_lags: The number of lags used for computing the semivariance.

  • coordinate_names: The column names for the spatial coordinates in the transect and mesh DataFrames.

These values can then be used to create the object vario:

vario = Variogram(
    lag_resolution=0.002,
    n_lags=30,
    coordinate_names=("x", "y")
)

The Variogram object has a variety of attributes that can be accessed with each subsequent method and computation.

# Lag vector
vario.lags

# Lag counts
vario.lag_counts

# Lagged semivariance
vario.gamma

# Lag covariance
vario.lag_covariance

# Initial variogram model parameters
vario.variogram_params_inital

# Initial variogram model parameters fit
vario.variogram_fit_initial

# Fitted variogram model parameters
vario.variogram_params_optimized

# Fitted variogram model parameters fit
vario.variogram_fit_optimized

Nothing appears, but that is okay. These will be populated with later steps. With our Variogram object initialized, we can then calculate the empirical variogram using the Variogram.empirical_variogram method. This method has four arguments:

  • data: The transect data.

  • variable: The column in the transect DataFrame that will be used for the spatial analysis. For instance, we can use "biomass_density" for this.

  • azimuth_filter: When True, azimuth angles between coordinate pairs are used for directional filtering. This is particularly useful for constraining spatial anisotropy.

  • azimuth_angle_threshold: When azimuth_filter is True, this is the maximum azimuth angle deviation (in degrees) allowed between paired coordinates.

vario.calculate_empirical_variogram(
    data=df_nasc_no_age1_prt, 
    variable="biomass_density", 
    azimuth_filter=True, 
    azimuth_angle_threshold=180
)

Once the empirical variogram is computed, we can take another look at the attributes to see what was updated. Specifically, we can investigate the lags, lag counts, and semivariance (gamma).

vario.lags

array([0.   , 0.002, 0.004, 0.006, 0.008, 0.01 , 0.012, 0.014, 0.016,
       0.018, 0.02 , 0.022, 0.024, 0.026, 0.028, 0.03 , 0.032, 0.034,
       0.036, 0.038, 0.04 , 0.042, 0.044, 0.046, 0.048, 0.05 , 0.052,
       0.054, 0.056, 0.058])

vario.lag_counts

array([  9286,   4907,   9426,   9278,   9156,  37034,  36076,  28286,
        26568,  58624,  60056,  49082,  44484,  75502,  83456,  66558,
        60376,  88618, 105106,  82948,  76114,  93692, 129330,  99294,
        88286,  87798, 159346, 110742, 100506,  94880])

vario.gamma

array([0.        , 0.03222788, 0.27160393, 0.49486823, 0.60184138,
       0.79318917, 0.83445331, 0.85618106, 0.88999828, 0.85965762,
       0.8763788 , 0.92765322, 0.94444364, 0.92116988, 0.93427382,
       0.94469502, 0.94983059, 0.91850362, 0.90343159, 0.93291779,
       0.96267877, 0.93406458, 0.9220397 , 0.96056632, 0.97118683,
       0.97548987, 0.95498628, 0.98826505, 0.98966192, 0.99734628])

With the empirical variogram calculated, the next step is to model the variogram. There are many different semivariogram models that can be used, and which parameters one uses can vary from model-to-model. Here, we will use a simple model: exponential. The Variogram.fit_variogram_model method is used for model fitting, which has three arguments:

  • model: The name of the model.

  • model_parameters: A lmfit.Parameters-class object with the initial values and constraints for each model parameter required for optimization.

  • optimizer_kwargs: An optional dictionary for passing keyword arguments used by lmfit.minimizer.minimize for optimizing the model parameters.

Since model="exponential" in this example, we need to fit the model using the nugget, sill, and correlation_range. The parameters are set to:

  • nugget: initial = 0.0, vary=True (i.e., this parameter will be varied), minimum = 0.0

  • sill: initial = 1.0, vary=True, minimum = 0.0

  • correlation_range: 0.01, vary=True, minimum = 0.001

from lmfit import Parameters

# Set up `lmfit` parameters
# lmfit.Parameters tuples: (NAME VALUE VARY MIN MAX EXPR BRUTE_STEP)
variogram_parameters_lmfit = Parameters()
variogram_parameters_lmfit.add_many(
    ("nugget", 0.0, True, 0.0),
    ("sill", 1.0, True, 0.0),
    ("correlation_range", 0.01, True, 0.001),
)

So now we feed variogram_parameters_lmfit into vario to optimize the variogram parameters:

best_fit_parameters = vario.fit_variogram_model(
    model="exponential",
    model_parameters=variogram_parameters_lmfit,
)

First let’s investigate the initial set or parameters and their model fit as a reference.

vario.variogram_params_inital

{'nugget': 0.0, 'sill': 1.0, 'correlation_range': 0.01}

vario.variogram_fit_initial

np.float64(0.0011709552080984554)

Now we can compare to the optimized values:

vario.variogram_params_optimized

{'nugget': 0.10873628362976354,
 'sill': 0.956889261897134,
 'correlation_range': 0.007600320029644778}

vario.variogram_fit_optimized

np.float64(0.0007938833723707372)

So in this case, we see a somewhat modest improvement, although the nugget estimate is quite high. We can also set some parameters to remain fixed (and therefore not be optimized). So if we do that here:

variogram_parameters_lmfit = Parameters()
variogram_parameters_lmfit.add_many(
    ("nugget", 1e-10, False),
    ("sill", 1.0, True, 0.0),
    ("correlation_range", 0.01, True, 0.001),
)

best_fit_parameters = vario.fit_variogram_model(
    model="exponential",
    model_parameters=variogram_parameters_lmfit,
)

We can see that we get slightly different estimates for the sill and correlation range:

vario.variogram_params_optimized

{'nugget': 1e-10,
 'sill': 0.9555106907984584,
 'correlation_range': 0.007014696146941477}

vario.variogram_fit_optimized

np.float64(0.0007725528458103382)

The model fit here represents the mean absolute deviation, so in this case the improvement from the optimization was relatively modest.

Kriging#

We can leverage the semivariogram model to interpolate our transect data via ordinary kriging. In Echopop, this is done using the Kriging-class from the geostatistics sub-package. There are four required arguments for initialization:

  • mesh: The kriging mesh DataFrame.

  • kriging_params: A dictionary with keyword arguments used by the built-in ordinary kriging functions.

  • variogram_params: A dictionary with variogram parameters. In this case, we will use the best-fit parameters.

  • coordinate_names: The names of the DataFrame coordinate columns.

There are four required keys for kriging_params:

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

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

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

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

So now we initialize a Kriging object:

from echopop.geostatistics import Kriging

# Define the requisite kriging parameters
KRIGING_PARAMETERS = {
    "search_radius": best_fit_parameters["correlation_range"] * 3,
    "aspect_ratio": 0.001,
    "k_min": 3,
    "k_max": 10,
}

# Define the requisite variogram parameters and arguments
VARIOGRAM_PARAMETERS = {"model": "exponential", **best_fit_parameters}

krg = Kriging(
    mesh=df_mesh,
    kriging_params=KRIGING_PARAMETERS,
    variogram_params=VARIOGRAM_PARAMETERS,
    coordinate_names=("x", "y"),
)

Once initialized, we can use the Kriging.krige method to perform the actual ordinary kriging interpolation and then project the results on the mesh grid. This method has six arguments:

  • transects: The DataFrame containing the transect data.

  • variable: The column name from transects that will be used for kriging. It is crucial that the variable used for the variogram-fitting is the same as the one used her for kriging.

  • extrapolate: When True (default), values from transects[variable] will be interpolated over the entire mesh grid. When False, a cropped version of the grid will be defined, which requires using the Kriging.crop_mesh method prior to using Kriging.krige.

  • default_mesh_area: If the mesh DataFrame does not have the column "area", then one can be provided here (in nmi2).

  • adaptive_search_strategy: This is an optional feature that allows for custom algorithms to be created by users. Otherwise, a uniform adaptive search strategy is used for the nearest neighbor algorithm.

  • custom_search_kwargs: These are keyword arguments specific to adaptive_search_strategy is a custom one was supplied.

So first we can run the algorithm over the entire mesh grid.

extrapolated_results = krg.krige(
    transects=df_nasc_no_age1_prt,
    variable="biomass_density",
    default_mesh_cell_area=6.25,
)

C:\Users\Brandyn\GitHub\echopop\echopop\geostatistics\kriging.py:1451: UserWarning: 142 invalid kriged estimates (< 0.0) found. These have been replaced with '0.0'; however, this truncation may distort error distributions. If this is problematic, adjust the choice of variogram and kriging algorithm parameterization.
  kriged_results, _ = project_kriging_results(
et_id area (km^2) latitude longitude fraction cell cut by polygon? geostratum_inpfc geostratum_ks x y area biomass_density kriged_variance sample_variance cell_cv
0 85670 21.4369 49.099727 -126.024144 1.0 N 6 7 0.125747 0.200900 6.25 0.000000 0.500427 NaN 0.044646
1 85427 21.4369 49.057959 -126.024127 1.0 N 6 7 0.119084 0.198853 6.25 0.000000 0.131102 NaN 0.022852
2 85184 21.4369 49.016196 -126.024110 1.0 N 6 7 0.109859 0.196807 6.25 0.000000 0.414271 NaN 0.040621
3 84941 21.4369 48.974438 -126.024093 1.0 N 6 7 0.099092 0.194760 6.25 0.000000 0.646779 NaN 0.050756
4 84698 21.4369 48.932686 -126.024076 1.0 N 6 7 0.088307 0.192714 6.25 45543.654721 0.732735 1.373507 0.054024
5 84455 21.4369 48.890939 -126.024060 1.0 N 6 8 0.077506 0.190669 6.25 0.000000 0.757142 NaN 0.054916
6 84212 21.4369 48.849198 -126.024043 1.0 N 6 8 0.066959 0.188623 6.25 54015.686837 0.718963 1.812516 0.053514
7 83969 21.4369 48.807461 -126.024026 1.0 N 6 8 0.056558 0.186578 6.25 0.000000 0.809815 NaN 0.056795
8 83726 21.4369 48.765730 -126.024009 1.0 N 6 8 0.046140 0.184533 6.25 258.880035 0.325678 6.237652 0.036017
9 83483 21.4369 48.724004 -126.023992 1.0 N 6 8 0.035707 0.182488 6.25 7458.941696 0.085441 0.242226 0.018448

This can then be converted to biomass since we have defined cell areas:

extrapolated_results["biomass"] = extrapolated_results["biomass_density"] * extrapolated_results["area"]
extrapolated_results["biomass"].sum()

np.float64(1734285929.5831578)

However, depending on the survey’s degree of coverage, it may not make sense to extrapolate over the entire grid region. In that case, we can use the Kriging.crop_mesh method, which defaults to using a convex hull around the survey to delimit the mesh grid. The crop_function can be parameterized with any user-defined cropping function with the only requirements that they 1) accept mesh as a keyword argument and 2) output a pandas.DataFrame with the coordinates columns intact. For the default case which uses geostatistics.hull_crop, it requires us to supply the transect DataFrame as a keyword argument. Additional keyword arguments for this function include:

  • num_nearest_transects: 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. This defaults to 3.

  • mesh_buffer_distance: Buffer distance in nautical miles applied to the survey polygon before filtering mesh cells. This ensures adequate coverage around the survey boundary. This defaults to 2.5 (nmi).

  • projection: EPSG projection code for the input coordinate system. Default is WGS84 ("epsg:4326").

  • coordinate_names: Names of the coordinate columns when using DataFrames. Expected format: (x_col, y_col), which defaults to ("longitude", "latitude").

So if we crop our mesh, we see that our mesh grid has been cut down substantially:

krg.crop_mesh(transects=df_nasc_no_age1_prt)
print(
    f"""
    Original mesh shape: {krg.mesh.shape}
    Cropped mesh shape: {krg.mesh_cropped.shape}
    """
)      

    Original mesh shape: (19843, 10)
    Cropped mesh shape: (11233, 10)

So now when we do the interpolation, we can see the total biomass when we do and do not extrapolate beyond the survey extentm, which in this example did not amount to a very large difference:

unextrapolated_results = krg.krige(
    transects=df_nasc_no_age1_prt,
    variable="biomass_density",
    extrapolate=False,
    default_mesh_cell_area=6.25,
)

# Get biomass
unextrapolated_results["biomass"] = unextrapolated_results["biomass_density"] * unextrapolated_results["area"]

# Compare
print(
    f"""
    Extrapolated total biomass: {extrapolated_results["biomass"].sum() * 1e-6} kmt
    Unextrapolated total biomass: {unextrapolated_results["biomass"].sum() * 1e-6} kmt
    Difference: {extrapolated_results["biomass"].sum() * 1e-6 - unextrapolated_results["biomass"].sum() * 1e-6} kmt
    """    
)

    Extrapolated total biomass: 1734.2859295831577 kmt
    Unextrapolated total biomass: 1713.0707058823502 kmt
    Difference: 21.215223700807428 kmt
    

C:\Users\Brandyn\GitHub\echopop\echopop\geostatistics\kriging.py:1451: UserWarning: 142 invalid kriged estimates (< 0.0) found. These have been replaced with '0.0'; however, this truncation may distort error distributions. If this is problematic, adjust the choice of variogram and kriging algorithm parameterization.
  kriged_results, _ = project_kriging_results(