From Shapefile to Feature Vector: A Practical Bridge Between GeoPandas and Scikit-Learn

The Gap Nobody Talks About

Most GIS tutorials stop at the map. You load your shapefile, you visualize your polygons, you run a dissolve or a spatial join, and you call it done. The spatial analysis world and the machine learning world both have mature, well-documented Python ecosystems — but the bridge between them is left to you to figure out, usually under deadline.

That bridge is what this article is about.

The core problem is structural. GeoPandas operates on geometries: Shapely objects living inside a geometry column, CRS-aware, topologically meaningful. Scikit-learn operates on numerical arrays: (n_samples, n_features) matrices where every value must be a float and spatial awareness is completely absent. Moving from one to the other is not just a data type conversion. It requires decisions about encoding, coordinate systems, feature selection, and what to do with the geometry column itself.

Get it wrong and your model silently trains on meaningless features — or fails with a cryptic dtype error.

Step 1: Load and Inspect Your GeoDataFrame

Start with something concrete. Say you have a shapefile of informal settlements across Nairobi, with polygon geometries and attributes like area_sqm, pop_estimate, distance_to_road, and a land-use classification label you want to predict.

import geopandas as gpd
import pandas as pd
import numpy as np
gdf = gpd.read_file("nairobi_settlements.shp")
print(gdf.dtypes)
print(gdf.crs)

Two things to check immediately: the dtypes of your attribute columns (are there mixed types, objects, or nulls?), and the CRS. The CRS matters more than people think at this stage — if you plan to derive any geometric features (area, perimeter, centroid coordinates), you need a projected CRS, not WGS84. A polygon’s .area in EPSG:4326 gives you degrees-squared, which is functionally useless.

# Project to a metric CRS before extracting geometric features
gdf_proj = gdf.to_crs(epsg=32737) # UTM Zone 37S - appropriate for Nairobi
gdf_proj["area_m2"] = gdf_proj.geometry.area
gdf_proj["perimeter_m"] = gdf_proj.geometry.length
gdf_proj["compactness"] = (4 * np.pi * gdf_proj["area_m2"]) / (gdf_proj["perimeter_m"] ** 2)

The compactness ratio (the Polsby-Popper score, if you want the formal name) is a good illustrative example of what “feature engineering from geometry” actually looks like. You are not feeding Shapely objects to your model. You are extracting numerical signals from them first.

Step 2: Extract Geometric Features Intentionally

The geometry column contains far more information than most practitioners extract. Beyond area and perimeter, consider:

Centroid coordinates — useful if spatial position itself is a signal (it often is in land-use classification).

gdf_proj["centroid_x"] = gdf_proj.geometry.centroid.x
gdf_proj["centroid_y"] = gdf_proj.geometry.centroid.y

Bounding box dimensions — the width-to-height ratio of a polygon’s envelope can distinguish organic settlement boundaries from planned plot layouts.

bounds = gdf_proj.geometry.bounds # returns minx, miny, maxx, maxy
gdf_proj["bbox_width"] = bounds["maxx"] - bounds["minx"]
gdf_proj["bbox_height"] = bounds["maxy"] - bounds["miny"]
gdf_proj["elongation"] = gdf_proj["bbox_width"] / gdf_proj["bbox_height"].replace(0, np.nan)

Vertex count — a proxy for polygon complexity. Regularized administrative boundaries have few vertices; digitized informal areas from satellite imagery can have hundreds.

gdf_proj["vertex_count"] = gdf_proj.geometry.apply(lambda geom: len(geom.exterior.coords))

The principle here is that geometry is data. Treat it that way. Before you drop the geometry column, extract everything from it that might carry predictive signal.

Step 3: Handle Categorical Attributes

Shapefiles often carry categorical string columns — land-use class, administrative zone, building type. Scikit-learn will not accept strings. You have two sensible options: one-hot encoding for low-cardinality nominals, and ordinal encoding when the category has a meaningful order.

from sklearn.preprocessing import LabelEncoder, OneHotEncoder

import pandas as pd
# One-hot encode a nominal categorical
gdf_encoded = pd.get_dummies(gdf_proj, columns=["zone_type"], drop_first=True)
# Label-encode the target variable
le = LabelEncoder()
gdf_encoded["label"] = le.fit_transform(gdf_proj["land_use_class"])

One-hot encoding inflates your feature count if you have many unique categories. If a column like admin_ward has 30+ unique values, think carefully before dummying it — you may want to encode it as a hierarchical numeric index, aggregate by it, or simply drop it if it carries too much cardinality relative to your sample size.

Step 4: Build the Feature Matrix

Now you assemble your X and y arrays. This is the moment you drop the geometry column. It cannot go into Scikit-learn — not because GeoPandas and Scikit-learn are incompatible in spirit, but because a Shapely geometry is not a float.

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
feature_cols = [
"area_m2", "perimeter_m", "compactness", "elongation",
"vertex_count", "centroid_x", "centroid_y",
"distance_to_road", "pop_estimate",
"zone_type_residential", "zone_type_commercial" # one-hot encoded
]
X = gdf_encoded[feature_cols].values
y = gdf_encoded["label"].values
# Scale - especially important if mixing geometric features (large values) with ratios (0–1)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
X_train, X_test, y_train, y_test = train_test_split(
X_scaled, y, test_size=0.2, random_state=42, stratify=y
)

The stratify=y argument is worth highlighting. With spatial data, class imbalance is common — especially in land-use scenarios where one category dominates the landscape. Stratified splitting preserves the class distribution in both train and test splits, which gives you more honest evaluation metrics.

Step 5: Train, Evaluate, and Bring Results Back to the Map

This is where most ML tutorials end — with a confusion matrix floating in a notebook with no spatial context. The pattern worth building into your workflow is to always bring predictions back into a GeoDataFrame so you can look at them geographically.

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report
clf = RandomForestClassifier(n_estimators=150, random_state=42, class_weight="balanced")
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
print(classification_report(y_test, y_pred, target_names=le.classes_))

Now reattach predictions to the original geometry for spatial visualization:

Keep the original index alignment by predicting on the full dataset

gdf_proj["predicted_class"] = le.inverse_transform(clf.predict(X_scaled))
gdf_proj["prediction_confidence"] = clf.predict_proba(X_scaled).max(axis=1)
# Export back to file - predictions are now spatially anchored
gdf_proj[["geometry", "predicted_class", "prediction_confidence"]].to_file(
"nairobi_settlements_predicted.gpkg", driver="GPKG"
)

Loading this in QGIS or visualizing it with gdf_proj.plot(column=”predicted_class”, legend=True) immediately tells you something no confusion matrix can: where in space is your model failing? Misclassifications that cluster spatially point to missing spatial features — perhaps distance to a transit corridor, or elevation, or proximity to a water body that you have not encoded yet.

That spatial pattern of errors is information. Use it.

What to Watch Out For

Spatial autocorrelation in your train/test split. Random splitting assumes observations are independent. Spatial data violates this — polygons that are geographically close share environmental context. A random split can leak spatial information from training to test, inflating your accuracy. For serious work, use spatial cross-validation (sklearn-evaluation, libpysal, or a manual block split by administrative boundary).

Null geometry rows. GeoPandas will sometimes read shapefiles with null geometries, which will break your .area and .centroid extractions silently. Always run gdf = gdf[gdf.geometry.notna()] before feature extraction.

CRS consistency across joined layers. If you are enriching your GeoDataFrame with features from a second layer (road network, elevation raster, census polygons), verify that both are in the same CRS before the spatial join. The error you get when they are not is rarely informative.

Conclusion

The GeoPandas-to-Scikit-learn pipeline is not complicated, but it requires deliberate choices at every step: which CRS to project into before extracting geometric features, how to encode categoricals without introducing cardinality problems, and how to reassemble predictions back into spatial context after modeling.

What makes spatial ML genuinely different from tabular ML is not the geometry column — it is the obligation to think geographically about your results. A confusion matrix tells you your model is 78% accurate. A map of your predictions tells you it is 95% accurate in Westlands and completely lost in Kibera. That spatial specificity is the entire value proposition of bringing GIS thinking into your modeling pipeline.

The feature vector is not the end. It is the middle. The map is where you start and where you finish.

The code examples in this article use GeoPandas 0.14, Scikit-learn 1.4, and NumPy 1.26. Projections use UTM Zone 37S (EPSG:32737) as a metric reference for the Nairobi region.