Encoding Categorical Geographic Features

Categorical geographic variables are everywhere in spatial machine learning: land-cover class from a national land-use raster, soil or geology type from a survey polygon, administrative region from a boundary file, climate zone from a Köppen classification. Unlike continuous covariates such as elevation or reflectance, these variables carry no intrinsic numeric magnitude — yet almost every estimator expects a numeric feature matrix. How you convert them from labels to numbers determines whether the model learns a genuine geographic signal or silently memorises the answer.

The trap that catches most practitioners is not the encoding mechanics but where the encoding is fit. High-cardinality columns such as a municipality ID with thousands of levels invite mean or target encoding, and naive mean encoding fit on the whole dataset leaks the target straight into the features. Because geographic observations are spatially autocorrelated, that leakage is worse than in tabular data: nearby points share both a region label and a correlated target, so the encoded column becomes a near-perfect proxy for the answer. Validation looks spectacular; deployment to a new province collapses.

This guide is part of Spatial Feature Engineering for Machine Learning, and it sits alongside numeric transforms such as feature scaling for geospatial inputs. Here we cover the full categorical workflow: classifying columns by cardinality and order, one-hot and ordinal encoding for the simple cases, fold-safe target and frequency encoding for high-cardinality regions, verifying that no leakage occurred, and keeping unseen categories safe at inference time.

Prerequisites & Environment Setup

Encoding is deterministic, but the leakage guarantees depend on library behaviour that has changed across versions. Pin your stack so handle_unknown, sparse output, and the cross-validation encoders behave as documented here.

# Pinned requirements for categorical geographic encoding
scikit-learn==1.5.1
category_encoders==2.6.3
geopandas==0.14.4
pandas==2.2.2
numpy==1.26.4

Install with:

pip install "scikit-learn==1.5.1" "category_encoders==2.6.3" \
            "geopandas==0.14.4" "pandas==2.2.2" "numpy==1.26.4"

Load categorical layers and confirm dtypes

Read the vector layers that carry the categorical attributes and immediately cast the categorical columns away from any numeric dtype. Land-cover and soil surveys almost always ship their classes as integer codes, and pandas will happily treat them as continuous unless you intervene.

import geopandas as gpd
import pandas as pd

gdf = gpd.read_file("training_samples.gpkg").to_crs("EPSG:3035")

categorical_cols = ["landcover_code", "soil_type", "admin_region", "climate_zone"]
for col in categorical_cols:
    assert col in gdf.columns, f"Missing categorical column: {col}"
    gdf[col] = gdf[col].astype("string")   # force nominal, even for integer codes

print(gdf[categorical_cols].nunique())

The nunique() output is your single most important diagnostic. It tells you the cardinality of each column, which dictates the encoding strategy in the next section. A projected CRS such as EPSG:3035 is used here only so downstream geometry operations stay metric; encoding itself is CRS-agnostic, but keeping the workflow projected avoids surprises when you later join to a spatial cross-validation fold assignment.

Step-by-Step Implementation

Step 1 — Classify each column by cardinality and order

Before choosing an encoder, tabulate two properties for every categorical column: how many distinct levels it has, and whether those levels have a meaningful order. These two facts alone determine the correct method.

def profile_categoricals(df, cols, high_card_threshold=30):
    """Summarise cardinality and suggest an encoding strategy per column."""
    rows = []
    for col in cols:
        n = df[col].nunique(dropna=True)
        if n <= 2:
            strategy = "binary / one-hot"
        elif n <= high_card_threshold:
            strategy = "one-hot"
        else:
            strategy = "target or frequency (fit in CV folds)"
        rows.append({"column": col, "cardinality": n, "suggested": strategy})
    return pd.DataFrame(rows)

profile = profile_categoricals(gdf, categorical_cols)
print(profile)

A typical result: landcover_code has 12 levels (one-hot), soil_type has 8 (one-hot), climate_zone has 5 (one-hot), and admin_region has 1,400 levels (target or frequency). The threshold of 30 is a pragmatic default — one-hot beyond a few dozen columns starts to hurt tree ensembles and explode linear models.

Step 2 — One-hot encode low-cardinality nominal classes

For land-cover, soil type, and climate zone, one-hot encoding is the correct default. It makes no ordering assumption and every level becomes an independent binary indicator. The critical configuration is handle_unknown="ignore", which lets inference silently zero out a level that never appeared in training rather than raising.

from sklearn.preprocessing import OneHotEncoder

low_card_cols = ["landcover_code", "soil_type", "climate_zone"]

ohe = OneHotEncoder(
    handle_unknown="ignore",     # unseen level -> all-zero row, no crash
    sparse_output=False,
    dtype="float32",
)
ohe.fit(gdf[low_card_cols])

encoded = ohe.transform(gdf[low_card_cols])
feature_names = ohe.get_feature_names_out(low_card_cols)
print(f"One-hot produced {encoded.shape[1]} columns from {len(low_card_cols)} inputs")

The number of output columns equals the total number of distinct levels across the input columns. Keep an eye on it: if a supposedly low-cardinality column contributes hundreds of columns, it belongs in Step 4, not here.

Step 3 — Ordinal encode only truly ordered variables

Ordinal encoding replaces each level with a single integer. That single column is compact and lossless only when the integers respect a real order. Soil drainage class (poor < moderate < well) or a fire-danger rating are legitimately ordinal. A region name or a nominal land-cover code is not — assigning it integers invents a false distance the model will exploit incorrectly.

from sklearn.preprocessing import OrdinalEncoder

# Explicit, domain-defined order — never rely on alphabetical default
drainage_order = [["poor", "imperfect", "moderate", "well", "excessive"]]

ord_enc = OrdinalEncoder(
    categories=drainage_order,
    handle_unknown="use_encoded_value",
    unknown_value=-1,            # unseen class -> sentinel, handled downstream
)
gdf["drainage_ord"] = ord_enc.fit_transform(gdf[["soil_drainage"]])

assert gdf["drainage_ord"].between(-1, len(drainage_order[0]) - 1).all()

Always pass an explicit categories list. The default sorts levels alphabetically, which would encode excessive as 0 and poor as 3 — the exact inverse of the intended physical order.

Step 4 — Target encode high-cardinality regions inside cross-validation folds

This is where leakage is created or prevented. Target encoding replaces each region with the mean target value of the rows in that region. Fit it on the whole dataset and every row sees its own target through the group mean; because spatial neighbours share both region and target, the encoded feature becomes a proxy for the answer. The only safe pattern is to fit the encoder inside each cross-validation fold, on the training rows only, and transform the held-out rows with statistics they never contributed to.

Use a smoothed encoder so small regions shrink toward the global prior instead of trusting a handful of samples:

import numpy as np
import pandas as pd

class SmoothedTargetEncoder:
    """Mean target encoder with additive smoothing and a global-prior fallback
    for unseen categories. Fit on the training fold ONLY."""

    def __init__(self, smoothing: float = 20.0):
        self.smoothing = smoothing
        self.mapping_ = {}
        self.prior_ = 0.0

    def fit(self, s: pd.Series, y: np.ndarray) -> "SmoothedTargetEncoder":
        df = pd.DataFrame({"cat": s.values, "y": y})
        self.prior_ = df["y"].mean()
        stats = df.groupby("cat")["y"].agg(["mean", "count"])
        # shrink group mean toward the prior by sample count
        smooth = (stats["count"] * stats["mean"] + self.smoothing * self.prior_)
        smooth /= (stats["count"] + self.smoothing)
        self.mapping_ = smooth.to_dict()
        return self

    def transform(self, s: pd.Series) -> np.ndarray:
        # unseen category -> global prior, never NaN, never the row's own target
        return s.map(self.mapping_).fillna(self.prior_).to_numpy(dtype="float32")

Now wire it through a spatial splitter so the encoder is refit on every fold. The groups array here holds spatial block IDs, exactly as produced in the spatial cross-validation strategies workflow — this is what makes the encoding fold spatially honest rather than merely random.

from sklearn.model_selection import GroupKFold
from sklearn.ensemble import HistGradientBoostingRegressor

X_num = gdf[["elevation", "slope"]].to_numpy("float32")
region = gdf["admin_region"]
y = gdf["target"].to_numpy("float32")
groups = gdf["block_id"].to_numpy()      # spatial block assignment

gkf = GroupKFold(n_splits=5)
oof_pred = np.full(len(gdf), np.nan, dtype="float32")

for train_idx, val_idx in gkf.split(X_num, y, groups=groups):
    enc = SmoothedTargetEncoder(smoothing=20.0).fit(
        region.iloc[train_idx], y[train_idx]      # fit on TRAIN fold only
    )
    reg_train = enc.transform(region.iloc[train_idx]).reshape(-1, 1)
    reg_val = enc.transform(region.iloc[val_idx]).reshape(-1, 1)

    Xtr = np.hstack([X_num[train_idx], reg_train])
    Xva = np.hstack([X_num[val_idx], reg_val])

    model = HistGradientBoostingRegressor(max_iter=300, learning_rate=0.05,
                                          random_state=42)
    model.fit(Xtr, y[train_idx])
    oof_pred[val_idx] = model.predict(Xva)

rmse = np.sqrt(np.mean((oof_pred - y) ** 2))
print(f"Leakage-safe out-of-fold RMSE: {rmse:.3f}")

The encoder is instantiated inside the loop. There is no single global region_te column anywhere in the pipeline — that absence is the whole point.

Step 5 — Frequency encoding for a leakage-proof high-cardinality option

Frequency encoding replaces each region with how often it occurs. It never touches the target, so it cannot leak, and it is a strong signal when class prevalence correlates with the outcome (dense urban municipalities behave differently from sparse rural ones). It is the safest first move for high-cardinality columns and pairs well with target encoding.

def frequency_encode(train_s: pd.Series, other_s: pd.Series):
    """Map categories to their training-fold relative frequency.
    Fit on train, apply to any split; unseen -> 0."""
    freq = train_s.value_counts(normalize=True)
    return (train_s.map(freq).fillna(0.0).to_numpy("float32"),
            other_s.map(freq).fillna(0.0).to_numpy("float32"))

Compute frequencies on the training fold only, for the same reason target statistics must be: a global frequency table quietly encodes information about the held-out geography.

Step 6 — Keep unseen categories safe at inference

Production data will contain regions, soil types, and land-cover classes that never appeared in training. Fit encoders once on all training data, persist them, and guarantee that transform maps every unseen level to a defined fallback rather than crashing or producing NaN.

import joblib

# Fit final encoders on ALL training data, then persist
final_ohe = OneHotEncoder(handle_unknown="ignore", sparse_output=False,
                          dtype="float32").fit(gdf[low_card_cols])
final_te = SmoothedTargetEncoder(smoothing=20.0).fit(gdf["admin_region"], y)
joblib.dump({"ohe": final_ohe, "te": final_te, "cols": low_card_cols},
            "encoders.joblib")

def encode_for_inference(new_gdf, artifacts_path="encoders.joblib"):
    art = joblib.load(artifacts_path)
    for col in art["cols"] + ["admin_region"]:
        new_gdf[col] = new_gdf[col].astype("string")
    ohe_out = art["ohe"].transform(new_gdf[art["cols"]])          # unseen -> zeros
    te_out = art["te"].transform(new_gdf["admin_region"]).reshape(-1, 1)  # -> prior
    return np.hstack([ohe_out, te_out])

Verification & Testing

Two properties must hold before you trust the encoded features: no target leakage, and deterministic handling of unseen categories.

Assert no target leakage in target encoding. The definitive test is that fold-wise encoding underperforms naive whole-dataset encoding. If a naive global target encoder produces suspiciously better validation than the fold-safe version, the gap is leakage, not skill.

def leakage_probe(region, y, groups, smoothing=20.0):
    """Compare naive (leaky) vs fold-safe target-encoding signal.
    A large positive gap means the naive encoder is leaking."""
    naive = SmoothedTargetEncoder(smoothing).fit(region, y).transform(region)
    naive_corr = np.corrcoef(naive, y)[0, 1]

    gkf = GroupKFold(n_splits=5)
    oof = np.full(len(y), np.nan, dtype="float32")
    for tr, va in gkf.split(region, y, groups=groups):
        enc = SmoothedTargetEncoder(smoothing).fit(region.iloc[tr], y[tr])
        oof[va] = enc.transform(region.iloc[va])
    safe_corr = np.corrcoef(oof, y)[0, 1]

    print(f"Naive corr: {naive_corr:.3f}  |  fold-safe corr: {safe_corr:.3f}")
    assert naive_corr - safe_corr < 0.25, \
        "Naive encoding leaks strongly — always use the fold-safe encoder"

leakage_probe(region, y, groups)

Assert unseen-category handling. Feed the fitted encoders a category that is guaranteed absent and confirm the output is finite and the column count is unchanged.

import pandas as pd, numpy as np

phantom = pd.Series(["__NEVER_SEEN__"] * 3)
te_out = final_te.transform(phantom)
assert np.isfinite(te_out).all(), "Target encoder produced NaN on unseen category"
assert np.allclose(te_out, final_te.prior_), "Unseen category must map to prior"

ohe_out = final_ohe.transform(pd.DataFrame(
    {c: ["__NEVER_SEEN__"] * 3 for c in low_card_cols}))
assert ohe_out.shape[1] == len(final_ohe.get_feature_names_out(low_card_cols))
assert (ohe_out.sum(axis=1) == 0).all(), "Unseen row should be all zeros"

Troubleshooting & Common Errors

Target leakage from naive mean encoding

Symptom: Validation R² is near 0.95 during development but collapses on a new region. The admin_region target-encoded column dominates feature importance.

Cause: The encoder was fit on the whole dataset, so each row’s encoded value was computed from a group mean that included its own target. Spatial autocorrelation amplifies this because neighbours share the region label and correlated targets.

Fix: Refit the encoder inside every cross-validation fold on training rows only, as in Step 4. Confirm with the leakage_probe above. This is the same failure mode covered from the modelling side in handling spatial autocorrelation.

Category mismatch between train and inference

Symptom: ValueError: Found unknown categories at prediction time, or a feature matrix with a different column count than the model expects.

Cause: The encoder was fit with default handle_unknown="error", or a fresh encoder was fit on the inference batch instead of reusing the persisted training encoder.

Fix: Set handle_unknown="ignore" for one-hot and use_encoded_value for ordinal, persist the fitted encoder with joblib, and always transform (never fit) at inference. Assert ohe.get_feature_names_out() matches the training schema.

Integer land-cover codes wrongly treated as continuous

Symptom: A tree model splits landcover_code at 26.5, or a linear model reports a coefficient implying that code 42 is “larger” than code 11.

Cause: The column kept its integer dtype and entered the feature matrix as a magnitude, so the model inferred a false ordering and distance between classes.

Fix: Cast to string or category immediately after loading (Step 0), then one-hot encode. Verify with assert not pd.api.types.is_numeric_dtype(gdf["landcover_code"]).

Exploding dimensionality from one-hot on region IDs

Symptom: One-hot encoding a admin_region column adds 1,400 columns, training slows to a crawl, and a linear model overfits badly.

Cause: One-hot encoding scales with cardinality. High-cardinality identifiers create a sparse, wide matrix that starves tree ensembles of useful splits and gives linear models a coefficient per region to overfit.

Fix: Switch to fold-safe target or frequency encoding (Steps 4–5), or aggregate the ID to a coarser level (municipality to county) before encoding. Reserve one-hot for columns under the cardinality threshold from Step 1.

NaN encoded values after a spatial join

Symptom: Encoded columns contain NaN for a subset of rows, and the estimator raises during fit.

Cause: A prior zonal statistics or attribute join left some geometries without a matching categorical value, and the encoder mapped that missing key to NaN.

Fix: Fill missing categories with an explicit "__missing__" sentinel before encoding so they become their own level rather than propagating NaN into the model.

Rare categories dominate after target encoding

Symptom: Regions with two or three samples get extreme encoded values that swamp the signal.

Cause: With no smoothing, a group mean over a handful of rows is high-variance and unreliable.

Fix: Increase the smoothing parameter so small groups shrink toward the global prior. Values between 10 and 50 are reasonable; tune by out-of-fold RMSE.

Performance Optimisation

Prefer native categorical support in gradient boosters

Both LightGBM and HistGradientBoostingRegressor can consume categorical columns directly, splitting on category subsets without any one-hot expansion. For low-to-moderate cardinality this is faster and often more accurate than manual one-hot. Pass the categorical column indices and skip the encoder entirely — but note that native handling still learns splits from the target, so it must live inside the same spatial folds.

from sklearn.ensemble import HistGradientBoostingRegressor

model = HistGradientBoostingRegressor(
    categorical_features=[0, 1, 2],   # indices of landcover, soil, climate
    max_iter=300, learning_rate=0.05, random_state=42,
)

Use category dtype to cut memory on wide vector tables

Casting string columns to pandas category dtype stores each value as an integer code plus a small dictionary, which slashes memory on large point tables and speeds up groupby inside target encoding. This complements the numeric memory work in feature scaling for geospatial inputs.

for col in ["landcover_code", "soil_type", "admin_region", "climate_zone"]:
    gdf[col] = gdf[col].astype("category")

Vectorise fold-wise encoding with precomputed group statistics

Refitting a groupby on every fold is wasteful when folds share most rows. Precompute per-region sums and counts once, then subtract the held-out fold’s contribution to derive the training-fold statistics in O(1) per fold — a leave-one-fold-out trick that keeps the encoding leakage-safe while avoiding repeated full scans.

Combine encoded categoricals with neighbourhood context

Categorical encodings gain power when paired with spatial context features. A region’s target mean plus the spatial lag and neighborhood statistics of surrounding observations, or a land-cover class alongside raster band math and index calculation such as NDVI, often outperforms either family alone. Encode categoricals inside the same folds you use for those numeric features so the whole feature set stays consistent.

FAQ

Should I ever one-hot encode a high-cardinality region column?

Only after aggregating it to a coarser level. Raw one-hot on a municipality ID with over a thousand levels produces a sparse matrix that slows training and invites overfitting. Roll the IDs up to counties or agro-ecological zones first, or switch to smoothed target and frequency encoding fit inside spatial folds. Reserve one-hot for genuinely low-cardinality nominal variables such as land-cover class.

How much smoothing should I use for target encoding?

Smoothing controls how strongly small groups shrink toward the global prior. Start around 20 and tune by out-of-fold RMSE. Larger values are safer when many regions have few samples; smaller values let well-sampled regions express their true mean. The wrong extreme is zero smoothing, which lets a two-sample region produce a wildly overconfident encoded value.

Can I mix encodings for the same column?

Yes, and it often helps. Frequency encoding and target encoding capture different signals — prevalence versus outcome association — and a gradient booster can use both. Compute each on the training fold, add both as columns, and let the model weight them. Keep every target-derived column inside the cross-validation folds so none of them leak.

Does the encoder have to be refit for every spatial fold, or can I fit once?

Any encoder that reads the target — target, mean, leave-one-out, or native categorical splitting — must be refit on the training rows of each fold. Encoders that ignore the target — one-hot, ordinal, frequency — can be fit once on all training data, because they cannot leak the outcome. The dividing line is simple: if the encoder looked at y, it belongs inside the fold.

Related

• Feature Scaling for Geospatial Inputs — normalising numeric covariates that sit alongside your encoded categoricals
• Raster Band Math and Index Calculation — deriving spectral indices to pair with land-cover class features
• Zonal Statistics for Polygon Aggregation — summarising raster values within the regions you are encoding
• Spatial Lag and Neighborhood Statistics — neighbourhood context features that complement categorical encodings
• Spatial Cross-Validation Strategies — the fold structure that makes target encoding leakage-safe

Part of: Spatial Feature Engineering for Machine Learning