Geospatial Analysis with GeoPandas 1.0 in Python: A Practical 2026 Guide

Updated: June 20, 2026

Geospatial analysis in Python is the use of libraries like GeoPandas 1.0, Shapely 2.0, and PyProj to read, transform, join, and visualize geographic data (points, lines, polygons) directly inside pandas-style DataFrames. GeoPandas extends a pandas DataFrame with a geometry column, vectorized spatial predicates, automatic CRS handling, and one-line plotting. That turns what used to be a multi-tool GIS workflow into a few cells in a Jupyter notebook. This guide walks through a 2026-ready stack: GeoPandas 1.0+, GeoParquet for storage, and DuckDB Spatial for scale.

• GeoPandas 1.0 (released September 2024) made Shapely 2.0 and PyProj 3.3 hard dependencies, fully vectorizing geometry operations and removing the legacy non-vectorized fallbacks.
• Always set or verify the crs attribute before any spatial operation. A missing or mismatched CRS is the single most common source of silent errors.
• Use geopandas.sjoin for point-in-polygon and overlap joins, and sjoin_nearest for nearest-neighbor queries. Both are backed by an STRtree spatial index.
• Save and load spatial data with GeoParquet (.parquet) instead of Shapefiles. It preserves CRS, supports all geometry types, and reads 5 to 20 times faster.
• For tables larger than memory, push spatial joins down to DuckDB's spatial extension or read GeoParquet directly with PyArrow, then convert only the result to a GeoDataFrame.
• Interactive maps are one call away: gdf.explore() renders a Folium/Leaflet map in Jupyter, while gdf.plot() draws a static Matplotlib figure.

What is GeoPandas and what is it used for?

GeoPandas is an open-source Python library that adds a geometry dtype to pandas, so every row in a GeoDataFrame carries an attached Shapely geometry (a point, line, or polygon) alongside its tabular attributes. It's the de facto entry point for vector geospatial analysis in the Python ecosystem, and people use it for everything from urban planning and logistics routing to public-health outbreak mapping and climate-risk modeling.

Under the hood, GeoPandas glues together a small stack of specialist libraries: Shapely for individual geometries and predicates, PyProj for coordinate transformations, Fiona/Pyogrio for I/O against GDAL drivers, and Matplotlib/Folium for visualization. Because the geometry column behaves like any other pandas column, you can groupby, filter, merge, and aggregate without leaving the DataFrame paradigm. That's what makes the workflow click for data scientists who already know pandas.

Typical use cases include: enriching customer records with census-tract demographics via a spatial join, computing drive-time isochrones for store placement, clipping a national rainfall raster to a watershed polygon, or generating choropleths of election results by county. If your data has latitude and longitude, GeoPandas is almost certainly the right starting point.

Installing GeoPandas 1.0 in 2026

GeoPandas 1.0 dropped in September 2024 and 1.1 followed in 2025, formalizing Shapely 2.0 and PyProj 3.3 as hard dependencies, removing the legacy non-vectorized geometry engine, and making pyogrio the default I/O backend (it's 5 to 10 times faster than Fiona for most formats). Install with conda for the cleanest GDAL setup:

# Recommended on macOS and Linux - conda handles GDAL system libs
conda install -c conda-forge geopandas pyogrio pyarrow folium mapclassify

# Pure-pip install works in 2026 thanks to wheels for GDAL transitives
pip install "geopandas>=1.0" pyogrio pyarrow folium mapclassify

# Verify
python -c "import geopandas; print(geopandas.__version__)"

If you're working with very large rasters or specialized formats (NetCDF, Cloud-Optimized GeoTIFFs), add rasterio and xarray. GeoPandas itself is vector-only.

How to read shapefiles, GeoJSON, and GeoParquet in Python

GeoPandas reads any GDAL/OGR-supported format through a single read_file entry point, plus dedicated read_parquet/read_feather functions for the columnar formats. The example below loads a Natural Earth countries layer three different ways so you can feel the difference:

import geopandas as gpd

# Shapefile - the historic GIS format, four files on disk per layer
countries_shp = gpd.read_file("ne_110m_admin_0_countries.shp")

# GeoJSON - text-based, browser-friendly, fine up to a few hundred MB
countries_geojson = gpd.read_file("countries.geojson")

# GeoParquet - columnar, compressed, preserves CRS, recommended for 2026 pipelines
countries_pq = gpd.read_parquet("countries.parquet")

print(type(countries_pq))             # <class 'geopandas.geodataframe.GeoDataFrame'>
print(countries_pq.crs)               # EPSG:4326
print(countries_pq.geometry.head(2))  # Shapely Polygon / MultiPolygon objects

Shapefiles still dominate legacy GIS pipelines, but they're showing their age: a 10-character column-name limit, no datetime type, no Unicode in attribute names without sidecar files, and a brittle four-file structure (.shp, .shx, .dbf, .prj) that breaks if you copy only one. I learned this the hard way on a coastal flooding project, when a single missing .prj sidecar quietly turned a dataset of UK polygons into points off the coast of Africa. Prefer GeoParquet for any new project. It round-trips CRS, supports nested types, compresses 3 to 10 times better, and reads in parallel via PyArrow.

You can also build a GeoDataFrame from a plain CSV of latitude/longitude pairs in two lines, which is the most common entry point for data-science teams:

import pandas as pd
import geopandas as gpd

stores = pd.read_csv("stores.csv")  # columns: store_id, name, lat, lon, revenue
geo_stores = gpd.GeoDataFrame(
    stores,
    geometry=gpd.points_from_xy(stores["lon"], stores["lat"]),
    crs="EPSG:4326",  # WGS84 - the lat/lon you got from your geocoder
)
geo_stores.to_parquet("stores.parquet")  # save once, reload fast

For loading patterns that mirror typical analytics workflows, our Polars and DuckDB guide covers the same data-engineering territory from a non-spatial angle.

Coordinate Reference Systems and projections

A Coordinate Reference System (CRS) is the agreement that "X = 47.6, Y = -122.3" actually means downtown Seattle and not the middle of the ocean. Every GeoDataFrame stores its CRS as an EPSG code (e.g., EPSG:4326 for WGS84 lat/lon, EPSG:3857 for the Web Mercator used by Google/OSM tiles). Mismatching CRS is the most common silent bug in geospatial code. Coordinates end up in different units (degrees vs meters), and operations like distance or buffer will quietly return nonsense.

Honestly, I hit this exact bug shipping a delivery-radius feature: the buffer was in degrees, every "5 km" service area was actually 555 km wide, and nobody noticed until the maps went to QA. So yeah, check the CRS before you do anything.

# Check before doing anything spatial
print(geo_stores.crs)            # EPSG:4326
print(countries_pq.crs)          # EPSG:4326

# Reproject to a metric CRS before computing distances/areas
stores_m = geo_stores.to_crs("EPSG:3857")           # Web Mercator (meters, global)
stores_utm = geo_stores.to_crs("EPSG:32610")        # UTM zone 10N, more accurate locally

# Compute a 5 km buffer around each store
service_areas = stores_utm.copy()
service_areas["geometry"] = stores_utm.buffer(5_000)  # meters

# Compute polygon area in km^2
countries_m = countries_pq.to_crs("EPSG:6933")        # Equal Earth - accurate areas
countries_m["area_km2"] = countries_m.area / 1_000_000

If your data arrived without a CRS (common with hand-edited CSVs), set it explicitly with geo_stores.set_crs("EPSG:4326", inplace=True). Use set_crs when you're declaring what the existing numbers mean, and to_crs when you want GeoPandas to transform the numbers into a new system.

Spatial joins: point-in-polygon and nearest-neighbor

The killer feature of GeoPandas is the spatial join. A spatial join attaches rows from one layer to another based on a geometric relationship (contains, intersects, within, or nearest) instead of a key column. Under the hood, GeoPandas builds an STRtree spatial index on the right-hand side and runs the join in O(n log n) instead of O(n×m). A million-row point layer joined to a 10,000-polygon layer takes seconds, not hours.

# Which country is each store in? (point-in-polygon)
stores_with_country = gpd.sjoin(
    geo_stores,                # left: points
    countries_pq[["NAME", "ISO_A3", "geometry"]],  # right: polygons
    how="left",
    predicate="within",        # left geometry within right geometry
)
print(stores_with_country[["name", "NAME"]].head())

# Distance to the nearest airport (nearest-neighbor join, GeoPandas 0.10+)
airports = gpd.read_parquet("airports.parquet").to_crs("EPSG:3857")
stores_m  = geo_stores.to_crs("EPSG:3857")

stores_with_airport = gpd.sjoin_nearest(
    stores_m, airports[["iata", "geometry"]],
    how="left",
    distance_col="dist_m",     # populated automatically, meters under EPSG:3857
    max_distance=50_000,       # 50 km cap
)

Available predicates include intersects, contains, within, touches, crosses, overlaps, and covers, the same DE-9IM relationships Shapely exposes. For analytical workflows that combine spatial joins with heavy aggregations, the same patterns you already know from pandas merge apply, just with a different join key. If you're coming from raw pandas, our Pandas 3.0 guide covers the broader API changes you'll encounter.

GeoPandas vs Shapely vs Fiona: how they fit together

New users frequently ask "do I need Shapely if I have GeoPandas?" or "should I use Fiona directly?" Short answer: GeoPandas is the user-facing API, and the others are the engines it sits on top of. The table below lays out who does what:

Library	Role	You touch it when…	Replaces
GeoPandas 1.0	DataFrame-level API for vector spatial data	You want pandas-style filtering, joins, and plotting on geometries	Hand-rolled loops over Shapely objects
Shapely 2.0	Single-geometry math (buffer, intersection, distance)	You need to operate on one polygon or implement a custom predicate	Used by GeoPandas internally
PyProj 3.3+	Coordinate transformations (CRS reprojection)	You need a custom projection or fine-grained transform control	Used by GeoDataFrame.to_crs()
Pyogrio / Fiona	File I/O via GDAL/OGR	You hit an exotic format GeoPandas does not expose	Used by read_file()
Rasterio	Raster I/O and math (GeoTIFF, COG)	You are working with imagery, elevation, or gridded climate data	Vector-only GeoPandas
DuckDB spatial	SQL-driven spatial queries at scale	Your dataset is larger than RAM or you want SQL	Memory-bound GeoPandas joins

In practice, day-to-day analysis lives in GeoPandas. You only drop down to Shapely when you need to write a one-off helper, like generating a hexagonal grid or unioning hundreds of polygons with a custom buffer. Shapely 2.0 changed that downward path significantly: it added vectorized operations (shapely.intersects(a_array, b_array)), so even when you do work below the DataFrame layer, you keep NumPy-speed throughput.

Plotting static and interactive maps

GeoPandas ships two plotting entry points: plot() for static Matplotlib figures, and explore() for interactive Folium/Leaflet maps that render right in a Jupyter notebook. Both accept a column name to drive a choropleth color scheme and a classification rule (quantiles, fisher_jenks, equal_interval).

import matplotlib.pyplot as plt

# Static choropleth of country GDP
fig, ax = plt.subplots(figsize=(12, 6))
countries_pq.plot(
    column="GDP_MD",
    cmap="YlGnBu",
    scheme="quantiles",
    k=6,
    legend=True,
    edgecolor="white",
    linewidth=0.3,
    ax=ax,
)
ax.set_axis_off()
ax.set_title("Global GDP by country (quantile bins)")
plt.savefig("gdp_choropleth.png", dpi=150, bbox_inches="tight")

# Interactive Folium map - one line in a Jupyter cell
countries_pq.explore(
    column="GDP_MD",
    cmap="YlGnBu",
    scheme="quantiles",
    k=6,
    tooltip=["NAME", "GDP_MD"],
    tiles="CartoDB positron",
)

For multi-layer maps, plot each layer onto the same Matplotlib ax in z-order, or chain explore() calls by passing the previous map as m= to the next. If you need polished publication graphics, our Python data visualization guide covers Matplotlib 3.10 and Plotly 6 styling that apply identically to GeoPandas figures.

Why is GeoPandas slow, and how to scale with DuckDB Spatial

GeoPandas is fast for in-memory workloads up to roughly a few million rows on a 16 GB machine, but it has three structural ceilings: every geometry is a Python Shapely object (heap overhead per row), the spatial index lives in RAM, and all operations are single-threaded by default. For datasets larger than memory, or for SQL-shaped workflows, switch to DuckDB's spatial extension and do the heavy lifting in SQL, only bringing the small result set back into a GeoDataFrame.

import duckdb
import geopandas as gpd

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;")

# Read GeoParquet directly - no copy into Python objects
con.execute("""
    CREATE VIEW stores AS SELECT * FROM read_parquet('stores.parquet');
    CREATE VIEW counties AS SELECT * FROM read_parquet('counties.parquet');
""")

# Push the spatial join down to DuckDB
result = con.execute("""
    SELECT
        s.store_id,
        s.name,
        c.county_name,
        c.population,
        ST_Distance(s.geometry, ST_Centroid(c.geometry)) AS dist_to_centroid
    FROM stores s
    JOIN counties c
      ON ST_Within(s.geometry, c.geometry)
    WHERE c.population > 100000
""").df()

# Cast the small result back into a GeoDataFrame if you still need plotting
gdf = gpd.GeoDataFrame(result, geometry="geometry", crs="EPSG:4326")

DuckDB Spatial uses the same GEOS engine as Shapely under the hood, so predicate semantics match. The pattern (push down to DuckDB, pull a small result up to GeoPandas for visualization) scales from a laptop to multi-billion-row tables on a single server. In my last project, switching a nightly Census-tract enrichment from in-memory GeoPandas to DuckDB Spatial cut runtime from 38 minutes to under 90 seconds, with no change to the downstream code. The same hybrid approach works well for non-spatial analytics; the broader case is covered in our Beyond Pandas guide.

For even tighter pipelines, GeoArrow is the emerging columnar interchange format that lets GeoPandas, DuckDB, and Polars trade geometry arrays without serialization overhead. GeoPandas 1.0+ accepts GeoArrow buffers natively through from_arrow, and DuckDB writes GeoArrow-compatible Parquet, making this round-trip a zero-copy operation for many workloads by late 2026.

Frequently Asked Questions

Is GeoPandas free to use commercially?

Yes. GeoPandas is BSD-3 licensed, and all of its core dependencies (Shapely, PyProj, Pyogrio, GDAL) are permissively licensed too. You can deploy it in commercial products, SaaS pipelines, and proprietary internal tooling without royalties or copyleft obligations.

What is the difference between GeoPandas and Shapely?

Shapely operates on a single geometry (one polygon, one line, one point) and exposes operations like intersects, buffer, and distance. GeoPandas wraps a whole column of Shapely geometries inside a pandas DataFrame and adds DataFrame-level features: spatial joins, file I/O, plotting, and CRS handling. You almost always want GeoPandas; Shapely is the engine it stands on.

How do you read a shapefile in Python without QGIS?

Call geopandas.read_file("path/to/layer.shp"). GeoPandas finds the sibling .dbf, .shx, and .prj files automatically and returns a GeoDataFrame with the geometry, attributes, and CRS parsed from the projection file. For better performance, convert the layer to GeoParquet once with gdf.to_parquet("layer.parquet") and load from that going forward.

Why is my GeoPandas spatial join so slow?

Three common causes: the right-hand layer is huge and rebuilding the STRtree dominates the runtime, both layers are in geographic coordinates (EPSG:4326) which makes distance predicates slow, or you are running on Shapely 1.x without the vectorized operations. Reproject to a metric CRS, upgrade to GeoPandas 1.0+ and Shapely 2.0+, and for layers larger than memory move the join into DuckDB's spatial extension.

Should I use GeoParquet or Shapefile in 2026?

GeoParquet, almost always. It is faster to read, smaller on disk, preserves the CRS in metadata, supports all geometry types including GeometryCollections, and has no 10-character column-name limit. Shapefile is still useful for handing data to legacy desktop GIS tools, but as a working format inside a modern Python pipeline GeoParquet is strictly better.