Merge pull request #273 from brandynlucca/fill_mesh_extrapolation_leaks

Changes to `kriging_analysis` and mesh cropping methods

echopop/analysis.py

@@ -39,7 +39,13 @@
     optimize_variogram,
 )
 from .statistics import stratified_transect_statistic
-from .utils.validate_dict import VariogramEmpirical
+from .utils.validate_dict import (
+    KrigingAnalysis,
+    KrigingParameters,
+    MeshCrop,
+    VariogramBase,
+    VariogramEmpirical,
+)
 
 
 def process_transect_data(
@@ -355,6 +361,26 @@ def krige(input_dict: dict, analysis_dict: dict, settings_dict: dict) -> tuple[p
         arguments and user-defined inputs.
     """
 
+    # Validate cropping method parameters
+    validated_cropping_methods = MeshCrop.create(**settings_dict["cropping_parameters"])
+    # ---- Update the dictionary
+    settings_dict["cropping_parameters"].update({**validated_cropping_methods})
+
+    # Validate the variogram parameters
+    valid_variogram_parameters = VariogramBase.create(**settings_dict["variogram_parameters"])
+    # ---- Update the dictionary
+    settings_dict["variogram_parameters"].update({**valid_variogram_parameters})
+
+    # Validate kriging parameters
+    valid_kriging_parameters = KrigingParameters.create(
+        **{**settings_dict["kriging_parameters"], **settings_dict["variogram_parameters"]}
+    )
+    # ---- Update the dictionary
+    settings_dict["kriging_parameters"].update({**valid_kriging_parameters})
+
+    # Validate the additional kriging arguments
+    _ = KrigingAnalysis(**settings_dict["variogram_parameters"])
+
     # Extract kriging mesh data
     mesh_data = input_dict["statistics"]["kriging"]["mesh_df"]
 
@@ -364,18 +390,22 @@ def krige(input_dict: dict, analysis_dict: dict, settings_dict: dict) -> tuple[p
     # Define the and prepare the processed and georeferenced transect data
     transect_data = edit_transect_columns(analysis_dict["transect"], settings_dict)
 
-    # Crop the mesh grid if the kriged data will not be extrapolated
-    if not settings_dict["extrapolate"]:
-        # ---- Compute the cropped mesh
-        mesh_full = crop_mesh(transect_data, mesh_data, settings_dict)
-        if (settings_dict["verbose"]) & (settings_dict["crop_method"] == "convex_hull"):
-            # ---- Print alert
-            print(
-                "Kriging mesh cropped to prevent extrapolation beyond the defined "
-                f"""`mesh_buffer_distance` value ({settings_dict['mesh_buffer_distance']} nmi)."""
-            )
+    # Add kriging parameters to the settings config
+    settings_dict.update(
+        {
+            "kriging_parameters": {
+                **input_dict["statistics"]["kriging"]["model_config"],
+                **valid_kriging_parameters,
+            },
+            "variogram_parameters": {
+                **settings_dict["variogram_parameters"],
+                **valid_variogram_parameters,
+            },
+        },
+    )
 
-    else:
+    # Crop the mesh grid if the kriged data will not be extrapolated
+    if settings_dict["extrapolate"]:
         # ---- Else, extract original mesh dataframe
         mesh_df = mesh_data.copy()
         # ---- Extract longitude column name
@@ -386,6 +416,17 @@ def krige(input_dict: dict, analysis_dict: dict, settings_dict: dict) -> tuple[p
         mesh_full = mesh_df.copy().rename(
             columns={f"{mesh_longitude}": "longitude", f"{mesh_latitude}": "latitude"}
         )
+    else:
+        # ---- Compute the cropped mesh
+        mesh_full = crop_mesh(transect_data, mesh_data, validated_cropping_methods)
+        if (settings_dict["verbose"]) and (
+            validated_cropping_methods["crop_method"] == "convex_hull"
+        ):
+            # ---- Print alert
+            print(
+                f"Kriging mesh cropped to prevent extrapolation beyond the defined "
+                f"`mesh_buffer_distance` value ({settings_dict['mesh_buffer_distance']} nmi)."
+            )
 
     # Standardize the x- and y-coordinates, if necessary
     if settings_dict["standardize_coordinates"]:
@@ -396,8 +437,8 @@ def krige(input_dict: dict, analysis_dict: dict, settings_dict: dict) -> tuple[p
         if settings_dict["verbose"]:
             # ---- Print alert
             print(
-                """Longitude and latitude coordinates (WGS84) converted to standardized """
-                """coordinates (x and y)."""
+                "Longitude and latitude coordinates (WGS84) converted to standardized "
+                "coordinates (x and y)."
             )
     else:
         # ---- Else, duplicate the transect longitude and latitude coordinates as 'x' and 'y'

echopop/spatial/mesh.py

@@ -10,7 +10,7 @@
 from ..spatial.transect import transect_bearing, transect_extent
 
 
-def crop_mesh(transect_data: pd.DataFrame, mesh_data: pd.DataFrame, settings_dict: dict):
+def crop_mesh(transect_data: pd.DataFrame, mesh_data: pd.DataFrame, cropping_parameters: dict):
     """
     Crop survey kriging mesh.
 
@@ -20,7 +20,7 @@ def crop_mesh(transect_data: pd.DataFrame, mesh_data: pd.DataFrame, settings_dic
         Georeferenced transect data.
     mesh_data: pd.DataFrame
         Kriging mesh.
-    settings_dict: dict
+    cropping_parameters: dict
         Dictionary containing relevant algorithm variables and arguments.
     """
     # Rename the mesh coordinate names, if necessary
@@ -35,11 +35,11 @@ def crop_mesh(transect_data: pd.DataFrame, mesh_data: pd.DataFrame, settings_dic
 
     # Select and return cropped mesh depending on the cropping method
     # ---- Interpolation
-    if settings_dict["crop_method"] == "transect_ends":
-        return transect_ends_crop_method(transect_data.copy(), mesh, settings_dict)
+    if cropping_parameters["crop_method"] == "transect_ends":
+        return transect_ends_crop_method(transect_data.copy(), mesh, cropping_parameters)
     # ---- Convex hull
-    elif settings_dict["crop_method"] == "convex_hull":
-        return hull_crop_method(transect_data.copy(), mesh, settings_dict)
+    elif cropping_parameters["crop_method"] == "convex_hull":
+        return hull_crop_method(transect_data.copy(), mesh, cropping_parameters)
 
 
 def hull_crop_method(transect_data: pd.DataFrame, mesh_data: pd.DataFrame, settings_dict: dict):
@@ -89,7 +89,7 @@ def hull_crop_method(transect_data: pd.DataFrame, mesh_data: pd.DataFrame, setti
 
 
 def transect_ends_crop_method(
-    transect_data: pd.DataFrame, mesh_data: pd.DataFrame, settings_dict: dict
+    transect_data: pd.DataFrame, mesh_data: pd.DataFrame, cropping_parameters: dict
 ):
     """
     Crop the kriging mesh by interpolating the eastern and western boundaries of the survey
@@ -101,15 +101,15 @@ def transect_ends_crop_method(
         Georeferenced transect data.
     mesh_data: pd.DataFrame
         Kriging mesh.
-    settings_dict: dict
+    cropping_parameters: dict
         Dictionary containing relevant algorithm variables and arguments.
     """
 
     # Extract the analysis settings
     # ---- Number of nearest transects
-    latitude_resolution = settings_dict["latitude_resolution"]
+    latitude_resolution = cropping_parameters["latitude_resolution"]
     # ---- Grid buffer distance (nmi)
-    bearing_tolerance = settings_dict["bearing_tolerance"]
+    bearing_tolerance = cropping_parameters["bearing_tolerance"]
 
     # Convert latitude resolution to degrees latitude
     latitude_resolution_deg = latitude_resolution / 60.0
@@ -272,7 +272,7 @@ def transect_ends_crop_method(
             )
         # -------- Append the indices
         region_2_index.append(idx[0])
-    # ---- Region 2
+    # ---- Region 3
     region_3_index = []
     # -------- Compute the change in longitude (degrees)
     delta_longitude = latitude_resolution_deg * np.cos(np.radians(region_3_latitude))
@@ -291,17 +291,24 @@ def transect_ends_crop_method(
         if np.isnan(region_3_extents[0][i]) | np.isnan(region_3_extents[1][i]):
             # -------- Compute the indices for the northern- and southernmost coordinates
             # -------- North
-            lat_w_max = np.argmax(transect_west["latitude"])
+            # lat_w_max = np.argmax(transect_west["latitude"])
+            lat_w_max = np.argmin(transect_west["latitude"])
             # -------- South
-            lat_e_max = np.argmax(transect_east["latitude"])
+            # lat_e_max = np.argmax(transect_east["latitude"])
+            lat_e_max = np.argmin(transect_east["latitude"])
             # -------- Slope
-            slope = (
+            # slope = (
+            #     transect_west["longitude"].iloc[lat_w_max]
+            #     - transect_east["longitude"].iloc[lat_e_max]
+            # ) / (transect_west["latitude"].min() - transect_east["latitude"].min())
+            slope = (transect_west["latitude"].min() - transect_east["latitude"].min()) / (
                 transect_west["longitude"].iloc[lat_w_max]
                 - transect_east["longitude"].iloc[lat_e_max]
-            ) / (transect_west["latitude"].max() - transect_east["latitude"].max())
+            )
             # -------- Set a new border threshold
             longitude_slope_i = (
-                slope * (region_3_latitude[i] - transect_east["latitude"].max())
+                # slope * (region_3_latitude[i] - transect_east["latitude"].max())
+                slope * (region_3_latitude[i] - transect_east["latitude"].min())
                 + transect_east["longitude"].iloc[lat_e_max]
             )
             if np.isnan(region_3_extents[0][i]):