TopferPillewizer.py

"""
Töpfer & Pillewizer (1966) map generalisation scaling.

TN = ON * (OS/TS) ^ (x/2)

  ON  — original feature count
  OS  — original map scale denominator  (e.g. 50 000)
  TS  — target  map scale denominator  (e.g. 500 000)
  x   — 1 (points), 2 (lines), 3 (polygons)

The user supplies OS/TS as a ratio in the range (0, 1] for reduction to a
smaller scale (OS/TS < 1 means the target denominator is larger, i.e. coarser).

Iteration schedule
------------------
If the supplied increment is smaller than the target ratio the module builds a
sequence of intermediate ratios starting at `increment`, with each successive
step 1.5× the previous step size, until the target ratio is reached.
"""

import logging
import math
import os
import tempfile
import shutil
import time

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.ops import unary_union


# ---------------------------------------------------------------------------
# Core T&P calculation
# ---------------------------------------------------------------------------

def topfer_count(on, ratio, x):
    """Return integer target count TN = ON * (OS/TS)^(x/2).  Minimum 1."""
    if on <= 0 or ratio <= 0:
        return 0
    return max(1, round(on * (ratio ** (x / 2.0))))


def dipdir_to_strike(dipdir_values, original_values=None):
    """
    Convert averaged dip-direction values back to strike.

    Input strike is converted to dip direction as ``dipdir = strike - 90`` before
    the grid/Kent engines average it, so the inverse is ``strike = dipdir + 90``.
    If the original strikes all lie within 0-180 the result is kept in 0-180
    (strike is axial), otherwise in 0-360.
    """
    strike = (pd.to_numeric(dipdir_values, errors='coerce') + 90) % 360
    if original_values is not None:
        orig = pd.to_numeric(original_values, errors='coerce').dropna()
        if len(orig) and orig.max() <= 180:
            strike = strike % 180
    return strike.round(1)


def scale_iterations(target_ratio, increment):
    """
    Build the sequence of OS/TS ratios for iterative T&P scaling.

    Steps are additive: ratio_n = 1 - n * increment.
    The sequence runs from (1 - increment) down to target_ratio.

    Example: target_ratio=0.1, increment=0.1
      0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1  (9 steps)

    If increment <= 0 or the first step already reaches the target,
    a single-element list [target_ratio] is returned.
    """
    if increment <= 0:
        return [target_ratio]

    ratios = []
    current = round(1.0 - increment, 12)
    while current > target_ratio:
        ratios.append(current)
        current = round(current - increment, 12)

    ratios.append(target_ratio)
    return ratios


# ---------------------------------------------------------------------------
# Point scaling  (x = 1)
# ---------------------------------------------------------------------------

def _count_grid(engine, method, gdf, min_x, max_x, min_y, max_y, gs):
    """Return the number of output cells of a grid engine method at cell size gs."""
    # The output cell count does not depend on the outlier removal, so count the cells
    # directly rather than running the method (which can be slow) at every trial size.
    if method in ('spherical_kent', 'outlier_removal') and hasattr(engine, 'count_cells'):
        return engine.count_cells(gdf, min_x, max_x, min_y, max_y, gs, min_points=1)
    if method == 'outlier_carmichael' and hasattr(engine, 'count_cells'):
        return engine.count_cells(gdf, min_x, max_x, min_y, max_y, gs, min_points=4)
    tmpdir = tempfile.mkdtemp()
    try:
        if method == 'gridcell_average':
            stem = engine.gridcell_average(gdf, min_x, max_x, min_y, max_y, n=gs, path_out=tmpdir)
        elif method == 'spherical_kent':
            stem = engine.spherical_kent(gdf, min_x, max_x, min_y, max_y, n=gs, path_out=tmpdir)
        elif method == 'outlier_carmichael':
            stem = engine.outlier_carmichael(gdf, min_x, max_x, min_y, max_y, n=gs, path_out=tmpdir)
        else:
            stem = engine.outlier_removal(gdf, min_x, max_x, min_y, max_y, n=gs, path_out=tmpdir)
        df = pd.read_csv(os.path.join(tmpdir, stem + '.csv'))
        return int((df['DIP'] != -999).sum()) if not df.empty else 0
    finally:
        shutil.rmtree(tmpdir, ignore_errors=True)


def _new_seed(info):
    """A random seed from the clock, noted in info so the result can be reproduced."""
    info['random_seed'] = int(time.time())
    return info['random_seed']


def scale_points_tp(gdf, ratio, engine, method, dip_col, dipdir_col,
                    grid_size=5000, firstorder_kwargs=None, info=None,
                    strike_input=False):
    """
    Reduce a point GeoDataFrame to TN = ON*(OS/TS)^0.5 using the chosen method.

    ON is always the raw point count (len(gdf)).  For stochastic, the fraction
    is drawn directly.  For grid methods, the grid cell size is found iteratively
    so that the number of occupied cells matches TN — the analytical estimate
    sqrt(total_area/TN) fails for non-uniform / sparse data, so a proportional
    seed + binary search is used instead.

    If ``info`` (a dict) is supplied it is filled with the parameters actually
    used (fraction, final grid size, first-order limits, fallbacks) for the
    T&P report.
    """
    from .FirstOrderOrientation import save_grid_to_shapefile

    if info is None:
        info = {}
    on = len(gdf)
    tn = topfer_count(on, ratio, x=1)
    info['method'] = method
    if tn >= on:
        info['note'] = 'target >= original count; all points kept'
        return gdf.copy()

    # --- First order: keep the TN best-ranked points (closest to a contact and
    # best aligned with it). Needs contact_gdf, dist_buffer and angle_tol in
    # firstorder_kwargs. If fewer than TN points pass the first-order filters,
    # all of them are returned.
    if method == 'firstorder':
        result = engine.firstorder(gdf, max_points=tn, **firstorder_kwargs)
        info['distance_limit'] = firstorder_kwargs.get('dist_buffer')
        info['angle_limit_deg'] = firstorder_kwargs.get('angle_tol')
        info['points_passing_limits'] = getattr(engine, 'last_firstorder_n_pass', None)
        info['selection'] = 'best-ranked by distance/limit + angle/limit'
        if len(result) < tn:
            info['note'] = f'only {len(result)} points pass the limits (target {tn})'
        if len(result) < tn:
            print(f"  T&P first order: only {len(result)} points pass the "
                  f"distance/angle filters, fewer than the target {tn}.")
        return result

    # --- Stochastic ---
    if method == 'stochastic':
        info['fraction_retained'] = tn / on
        tmpdir = tempfile.mkdtemp()
        try:
            result = engine.stochastic(gdf, frac=tn / on,
                                       random_state=_new_seed(info),
                                       path_out=tmpdir)
            return result.set_crs(gdf.crs, allow_override=True)
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)

    # --- Grid methods ---
    bounds = gdf.total_bounds
    min_x, max_x = int(bounds[0]), int(bounds[2]) + 1
    min_y, max_y = int(bounds[1]), int(bounds[3]) + 1

    is_geo = gdf.crs is not None and gdf.crs.is_geographic
    gs_ref = grid_size / 110000.0 if is_geo else float(grid_size)

    def count(gs):
        return _count_grid(engine, method, gdf, min_x, max_x, min_y, max_y, gs)

    # Seed with reference count to drive proportional estimate
    n_ref = count(gs_ref)
    if n_ref == 0:
        # No data at all — fall back to stochastic
        info['note'] = 'grid search found no data; fell back to stochastic'
        info['fraction_retained'] = max(0.001, tn / on)
        tmpdir = tempfile.mkdtemp()
        try:
            result = engine.stochastic(gdf, frac=max(0.001, tn / on),
                                       random_state=_new_seed(info),
                                       path_out=tmpdir)
            return result.set_crs(gdf.crs, allow_override=True)
        finally:
            shutil.rmtree(tmpdir, ignore_errors=True)

    # Proportional first estimate: N ∝ 1/gs²  →  gs ∝ sqrt(1/N)
    gs1 = gs_ref * math.sqrt(max(1, n_ref) / max(1, tn))
    n1 = count(gs1)

    tol = 0.15
    gs_n = n1
    if abs(n1 - tn) <= tol * tn:
        gs_best = gs1
    else:
        # Build binary-search brackets:
        #   fine_gs (small) → n >= tn;  coarse_gs (large) → n <= tn
        if n1 >= tn:
            fine_gs, fine_n = gs1, n1
            coarse_gs = gs1 * 2
            for _ in range(10):
                coarse_n = count(coarse_gs)
                if coarse_n <= tn:
                    break
                coarse_gs *= 2
        else:
            coarse_gs, coarse_n = gs1, n1
            fine_gs = gs1 / 2
            for _ in range(10):
                fine_n = count(fine_gs)
                if fine_n >= tn:
                    break
                fine_gs /= 2

        gs_best, best_n = gs1, n1
        for _ in range(12):
            if coarse_gs - fine_gs < 0.1:
                break
            mid = (fine_gs + coarse_gs) / 2
            mid_n = count(mid)
            if abs(mid_n - tn) < abs(best_n - tn):
                gs_best, best_n = mid, mid_n
            if abs(mid_n - tn) <= tol * tn:
                break
            if mid_n > tn:
                fine_gs = mid
            else:
                coarse_gs = mid
        gs_n = best_n

    to_m = 110000.0 if is_geo else 1.0   # report grid sizes in metres
    info['grid_size_start'] = gs_ref * to_m
    info['grid_size_final'] = gs_best * to_m
    info['cells_at_final_size'] = gs_n

    # Final run at best grid size
    tmpdir = tempfile.mkdtemp()
    try:
        if method == 'gridcell_average':
            stem = engine.gridcell_average(gdf, min_x, max_x, min_y, max_y,
                                           n=gs_best, path_out=tmpdir)
        elif method == 'spherical_kent':
            stem = engine.spherical_kent(gdf, min_x, max_x, min_y, max_y,
                                         n=gs_best, path_out=tmpdir)
        elif method == 'outlier_carmichael':
            stem = engine.outlier_carmichael(gdf, min_x, max_x, min_y, max_y,
                                             n=gs_best, path_out=tmpdir)
        else:
            stem = engine.outlier_removal(gdf, min_x, max_x, min_y, max_y,
                                          n=gs_best, path_out=tmpdir)
        result = save_grid_to_shapefile(tmpdir, stem)
    finally:
        shutil.rmtree(tmpdir, ignore_errors=True)

    if result is not None and not result.empty:
        result = result.rename(columns={'DIP': dip_col, 'DIP_DIR': dipdir_col}, errors='ignore')
        if strike_input and dipdir_col in result.columns:
            # engine averaged dip direction; restore the input's strike convention
            result[dipdir_col] = dipdir_to_strike(result[dipdir_col], gdf.get(dipdir_col))
        return result.set_crs(gdf.crs, allow_override=True)

    # Fallback to stochastic if grid method yielded nothing
    info['note'] = 'grid method produced no output; fell back to stochastic'
    info['fraction_retained'] = max(0.001, tn / on)
    tmpdir = tempfile.mkdtemp()
    try:
        result = engine.stochastic(gdf, frac=max(0.001, tn / on),
                                   random_state=_new_seed(info),
                                   path_out=tmpdir)
        return result.set_crs(gdf.crs, allow_override=True)
    finally:
        shutil.rmtree(tmpdir, ignore_errors=True)


# ---------------------------------------------------------------------------
# Line (fault) scaling  (x = 2)
# ---------------------------------------------------------------------------

def _edge_type_rank(edge_type):
    """
    Importance rank for graph edge_type (lower rank = higher importance = kept longer).

    x-x  → 0  (most important, crosses two X-nodes)
    z-z  → 3  (least important, isolated at both ends — removed first)
    else → 1  (x-z, y-*, etc.)
    """
    if pd.isna(edge_type):
        return 2
    et = str(edge_type).lower().strip()
    if et == 'x-x':
        return 0
    if et == 'z-z':
        return 3
    return 1


def _find_offset_column(columns):
    """
    Column holding the stratigraphic offset level from FaultStratOffset.

    The tool writes ``max_off_lv`` (0 = no offset, 1..n = offset detected at
    priority field k, -1 = fault at edge of map, empty = no data). That name
    contains neither 'strat' nor 'offset', so look for it explicitly before
    falling back to any column that looks like an offset/strat field.
    """
    for preferred in ('max_off_lv', 'max_offset_level'):
        if preferred in columns:
            return preferred
    return next((c for c in columns
                 if 'strat' in c.lower() or 'offset' in c.lower()), None)


def _find_cluster_column(columns):
    """
    Orientation-cluster label column written by the cluster tool ('cluster').
    Prefer that exact name; 'cluster_ce' (cluster centre) and any user field that
    merely contains 'cluster' are only used as a fallback.
    """
    if 'cluster' in columns:
        return 'cluster'
    return next((c for c in columns if 'cluster' in c.lower()
                 and 'cent' not in c.lower() and not c.lower().endswith('_ce')), None) \
        or next((c for c in columns if 'cluster' in c.lower()), None)


def _offset_rank(series):
    """
    Numeric rank for sorting faults by offset level, highest = kept longest.
    Edge-of-map (-1) and no-data (NaN) faults are ranked below 'no offset' (0).
    """
    vals = pd.to_numeric(series, errors='coerce')
    return vals.where(vals >= 0, -1).fillna(-1)


def scale_lines_tp(gdf, ratio, method, info=None):
    """
    Reduce a line GeoDataFrame to TN = ON*(OS/TS)^1 features.

    A ``tp_sel_val`` field is always appended to record the criterion used:
      length        — geometry length
      graph         — edge_type value  (falls back to length if field absent)
      strat_offset  — stratigraphic-offset value  (falls back to length)
      clusters      — cluster label  (falls back to length)

    Sorting rules by method:
      length        — keep longest features
      graph         — sort by edge_type rank (x-x > x-z/y-* > z-z),
                      secondary sort by length
      strat_offset  — keep highest absolute stratigraphic-offset value;
                      falls back to length if no offset field found
      clusters      — keep members of the largest orientation clusters;
                      falls back to length if no cluster field found
    """
    if info is None:
        info = {}
    on = len(gdf)
    tn = topfer_count(on, ratio, x=2)
    info['method'] = method

    def _to_str(series):
        """Convert a Series to string, replacing NaN/None with ''."""
        return ['' if pd.isna(v) else str(v) for v in series]

    # Fast path: all features kept — add fields without sorting
    if tn >= on:
        info['note'] = 'target >= current count; all faults kept'
        result = gdf.copy()
        if 'line_len' not in result.columns:
            result['line_len'] = result.geometry.length
        if method == 'graph' and 'edge_type' in result.columns:
            result['tp_sel_val'] = _to_str(result['edge_type'])
        elif method == 'strat_offset':
            strat_col = _find_offset_column(result.columns)
            result['tp_sel_val'] = _to_str(
                pd.to_numeric(result[strat_col], errors='coerce') if strat_col
                else pd.Series([''] * len(result), index=result.index))
        elif method == 'clusters':
            cluster_col = _find_cluster_column(result.columns)
            result['tp_sel_val'] = _to_str(result[cluster_col]) if cluster_col else ''
        else:  # length or graph-without-edge_type
            result['tp_sel_val'] = ('' if method == 'graph'
                                    else _to_str(result['line_len']))
        return result

    n_keep = min(tn, on)
    work = gdf.copy()
    # Reuse an existing line_len, except for the graph method: its features are
    # re-cut/merged graph edges, and a stored line_len copied from the original
    # fault would give every piece of a split fault the whole fault's length.
    if 'line_len' in work.columns and method != 'graph':
        work['_len'] = work['line_len']
    else:
        work['_len'] = work.geometry.length
        if method == 'graph' and 'line_len' in work.columns:
            work['line_len'] = work['_len']

    if method == 'graph' and 'edge_type' in work.columns:
        work['_rank'] = work['edge_type'].apply(_edge_type_rank)
        work = work.sort_values(['_rank', '_len'], ascending=[True, False])
        work['tp_sel_val'] = _to_str(work['edge_type'])
        info['criterion'] = 'edge_type rank (x-x > x-z/y-* > z-z), then length'

    elif method == 'graph':  # graph selected but edge_type column absent
        work = work.sort_values('_len', ascending=False)
        work['tp_sel_val'] = ''
        info['criterion'] = 'length (edge_type column absent - fallback)'

    elif method == 'strat_offset':
        strat_col = _find_offset_column(work.columns)
        if strat_col:
            # Highest offset level first; no-offset (0) next; edge-of-map / no-data
            # last; ties broken by length.
            work['_sort'] = _offset_rank(work[strat_col])
            work = work.sort_values(['_sort', '_len'], ascending=False)
            work['tp_sel_val'] = _to_str(pd.to_numeric(work[strat_col], errors='coerce'))
            info['criterion'] = (f'{strat_col} (offset level, highest first; none then '
                                 f'edge-of-map last), then length')
            kept_levels = work.head(n_keep)[strat_col].value_counts(dropna=False)
            info['offset_levels_kept'] = ', '.join(
                f'{"none" if pd.isna(k) else k}:{v}' for k, v in sorted(
                    kept_levels.items(), key=lambda kv: (pd.isna(kv[0]), kv[0])))
        else:
            work = work.sort_values('_len', ascending=False)
            work['tp_sel_val'] = ''
            info['criterion'] = 'length (no strat-offset field - fallback)'

    elif method == 'clusters':
        cluster_col = _find_cluster_column(work.columns)
        if cluster_col:
            sizes = work[cluster_col].value_counts()
            work['_csz'] = work[cluster_col].map(sizes)
            work = work.sort_values(['_csz', '_len'], ascending=[False, False])
            work['tp_sel_val'] = _to_str(work[cluster_col])
            info['criterion'] = f'size of {cluster_col} cluster, then length'
        else:
            work = work.sort_values('_len', ascending=False)
            work['tp_sel_val'] = ''
            info['criterion'] = 'length (no cluster field - fallback)'

    else:   # length
        work = work.sort_values('_len', ascending=False)
        work['tp_sel_val'] = _to_str(work['_len'])
        info['criterion'] = 'length'

    if 'line_len' not in work.columns:
        work['line_len'] = work['_len']

    drop_cols = [c for c in ['_len', '_rank', '_sort', '_csz'] if c in work.columns]
    kept = work.head(n_keep)
    info['shortest_fault_kept'] = float(kept['_len'].min()) if len(kept) else None
    return kept.drop(columns=drop_cols).reset_index(drop=True)


# ---------------------------------------------------------------------------
# Polygon scaling helpers
# ---------------------------------------------------------------------------

def _vals_match(a, b):
    """True if both values are non-null and equal."""
    try:
        if pd.isna(a) or pd.isna(b):
            return False
    except (TypeError, ValueError):
        pass
    return a == b


def _best_neighbor_idx(gdf, si, neighbors, priority_cols):
    """
    Return the positional index of the best merge target for the polygon at
    position *si*, honouring the attribute priority hierarchy.

    Parameters
    ----------
    gdf          : current working GeoDataFrame (reset_index applied)
    si           : positional index of the polygon to absorb
    neighbors    : list of (positional_index, shared_boundary_length) tuples
    priority_cols: column names tried in order (strat1..4, then lithoname)

    Within each priority tier the candidate with the *longest* shared boundary
    is preferred.  Falls back to longest boundary when no attribute match exists.
    """
    small_row = gdf.iloc[si]

    for col in priority_cols:
        matching = [
            (ci, sl) for ci, sl in neighbors
            if _vals_match(small_row[col], gdf.iloc[ci][col])
        ]
        if matching:
            return max(matching, key=lambda x: x[1])[0]

    return max(neighbors, key=lambda x: x[1])[0]


def _merge_to_count(gdf, target_count, priority_cols=None):
    """
    Count-controlled merge: repeatedly dissolve the smallest polygon into its
    best stratigraphically-matching neighbour until len(gdf) == target_count.

    Every polygon is *merged* (dissolved via union), never deleted, so no
    holes are created in the map.  The priority hierarchy is:

        strat1 match > strat2 match > strat3 match > strat4 match
        > lithoname match > longest shared boundary (fallback)  (as Rescale Map)

    Parameters
    ----------
    gdf           : input GeoDataFrame
    target_count  : desired number of polygons
    priority_cols : list of column names to check in priority order
    """
    if priority_cols is None:
        priority_cols = []

    result = gdf.copy().reset_index(drop=True)
    geom_col = result.geometry.name
    crs = result.crs

    while len(result) > target_count:
        result = result.reset_index(drop=True)
        areas = result.geometry.area.values
        si = int(np.argmin(areas))
        small_geom = result.geometry.iloc[si]

        # Bounding-box candidates via spatial index
        try:
            candidates = list(result.sindex.intersection(small_geom.bounds))
        except Exception:
            candidates = list(range(len(result)))
        candidates = [c for c in candidates if c != si]

        # Measure actual shared boundary length for touching candidates
        neighbors = []
        for ci in candidates:
            try:
                ni_geom = result.geometry.iloc[ci]
                if ni_geom.touches(small_geom) or ni_geom.intersects(small_geom):
                    shared = small_geom.intersection(ni_geom)
                    sl = shared.length if not shared.is_empty else 0.0
                    if sl > 0:
                        neighbors.append((ci, sl))
            except Exception:
                logging.getLogger(__name__).debug("Ignored exception", exc_info=True)

        if not neighbors:
            # Isolated polygon — nearest centroid as emergency fallback
            if candidates:
                best_ci = min(
                    candidates,
                    key=lambda c: small_geom.centroid.distance(
                        result.geometry.iloc[c].centroid)
                )
                neighbors = [(best_ci, 0.0)]
            else:
                # No candidates at all — remove to avoid infinite loop
                result = result.drop(index=si).reset_index(drop=True)
                continue

        # Select best neighbour honouring the attribute priority hierarchy
        best_ci = _best_neighbor_idx(result, si, neighbors, priority_cols)

        # Union geometries, drop the absorbed polygon
        geoms = result.geometry.tolist()
        geoms[best_ci] = unary_union([geoms[si], geoms[best_ci]])
        geoms.pop(si)   # si removed; if si < best_ci the updated geom shifts left

        result = (result
                  .drop(columns=[geom_col])
                  .drop(index=si)
                  .reset_index(drop=True))
        result = gpd.GeoDataFrame(result, geometry=geoms, crs=crs)

    return result


# ---------------------------------------------------------------------------
# Polygon scaling  (x = 3)
# ---------------------------------------------------------------------------

def scale_polygons_tp(gdf, ratio, polygon_subsampler_class,
                      distance_threshold=1e-6,
                      lithoname=None, strat1=None, strat2=None,
                      strat3=None, strat4=None,
                      dyke_field=None, dyke_codes=None,
                      triangulator_class=None, info=None):
    """
    Reduce a polygon GeoDataFrame to TN = ON*(OS/TS)^1.5 features.

    Strategy
    --------
    1. Run ``clean_small_polygons_and_holes_new`` once with a threshold sized
       to cover all the intended-to-remove polygons.  This efficiently handles
       the bulk of merges with full stratigraphic hierarchy.
    2. If the count is still above TN (the threshold-based algorithm stalls when
       small polygons are only adjacent to other small polygons), fall through to
       ``_merge_to_count``, which iteratively absorbs the smallest remaining
       polygon into its best stratigraphically-matching neighbour until the T&P
       target is reached.  Every polygon is merged (dissolved), never deleted,
       so no holes appear in the output map.
    """
    def _or_none(v):
        return v if v else None

    if info is None:
        info = {}
    on = len(gdf)
    tn = topfer_count(on, ratio, x=3)
    info['distance_threshold'] = distance_threshold
    info['priority_fields'] = ', '.join(
        c for c in [strat1, strat2, strat3, strat4, lithoname] if c) or 'none'
    if tn >= on:
        info['note'] = 'target >= current count; all polygons kept'
        return gdf.copy()

    work = gdf.copy()
    helper_columns = []

    # Optional dyke pre-processing (once, before merging)
    if (dyke_field and dyke_codes and triangulator_class is not None
            and dyke_field in work.columns):
        # A polygon is only triangulated if its area is below the threshold, so
        # the threshold must be the area below which this step will remove
        # polygons (the same cut the bulk merge below uses). With a threshold of
        # 0 nothing is triangulated and an ignore-code polygon such as a dyke is
        # merged whole into the neighbour with the longest boundary, painting that
        # unit across the dyke even where it crosses other units.
        _areas0 = np.sort(work.geometry.area.values)
        _n_remove0 = max(1, on - tn)
        tri_threshold = float(_areas0[min(_n_remove0 - 1, len(_areas0) - 1)]) * (1.0 + 1e-9) + 1e-6
        info['dyke_triangulation_threshold'] = tri_threshold
        trig = triangulator_class(
            gdf=work,
            id_column=dyke_field,
            min_area_threshold=tri_threshold,
            distance_threshold=distance_threshold,
            strat1=strat1, strat2=strat2, strat3=strat3, strat4=strat4,
            lithoname=lithoname,
        )
        work = trig.triangulate_polygons(target_ids=dyke_codes)
        helper_columns = getattr(trig, 'helper_columns', [])
        info['dyke_field'] = dyke_field
        info['dyke_codes'] = ', '.join(str(c) for c in dyke_codes)
        s = getattr(trig, 'last_summary', None)
        if s:
            info['dyke_polygons_matched'] = s['matched']
            info['dyke_polygons_triangulated'] = s['triangulated']
            if s['codes_not_found']:
                info['dyke_codes_not_found'] = ', '.join(s['codes_not_found'])

    merge_kwargs = dict(
        distance_threshold=distance_threshold,
        lithoname=_or_none(lithoname),
        strat1=_or_none(strat1),
        strat2=_or_none(strat2),
        strat3=_or_none(strat3),
        strat4=_or_none(strat4),
    )

    # Priority column list for count-controlled merge: Priority 1-4, then the lithology
    # (Priority 5), the same order as clean_small_polygons_and_holes_new (Rescale Map)
    priority_cols = [
        c for c in [strat1, strat2, strat3, strat4, lithoname]
        if c and c in work.columns
    ]

    # Pass 1: threshold-based bulk merge (efficient)
    n_current = len(work)
    if n_current > tn:
        n_remove = n_current - tn
        areas = work.geometry.area.sort_values().values
        # The n_remove smallest polygons are areas[0 .. n_remove-1]; the threshold
        # must sit just above the largest of those so exactly n_remove (not
        # n_remove + 1) polygons fall below it.
        idx = min(n_remove - 1, len(areas) - 1)
        threshold = float(areas[idx]) * (1.0 + 1e-9) + 1e-6
        info['min_area_threshold'] = threshold
        try:
            sub = polygon_subsampler_class(work)
            merged = sub.clean_small_polygons_and_holes_new(
                work, min_area_threshold=threshold, **merge_kwargs)
            if len(merged) < n_current:
                work = merged
        except Exception as exc:
            print(f"T&P polygon bulk-merge error: {exc}")
            info['note'] = f'bulk merge failed: {exc}'
        info['count_after_area_merge'] = len(work)

    # Pass 2: count-controlled merge for any remainder — no holes, priorities respected
    if len(work) > tn:
        info['count_merge_polygons'] = len(work) - tn
        work = _merge_to_count(work, tn, priority_cols=priority_cols)

    # remove the triangulation's bookkeeping columns once merging is done
    drop = [c for c in helper_columns if c in work.columns]
    if drop:
        work = work.drop(columns=drop)
    return work