CRS & Coordinate Systems: The Reference GIS Engineers Reach For

A coordinate reference system is the contract that says what a pair of numbers means on the surface of the Earth. Get that contract right and a dashboard’s layers stack perfectly, its distances read in metres, and its area statistics are trustworthy. Get it wrong and a district polygon renders in the sea, a “distance” comes back as a fraction of a degree, and a choropleth of area-per-capita is silently off by a factor that grows with latitude. This page is the reference for making that contract explicit in a Streamlit or Panel spatial dashboard: which EPSG code to use, how to inspect and convert it in GeoPandas and pyproj, and how to avoid the axis-order and datum traps that produce plausible-looking wrong answers.

CRS handling is unusually error-prone in dashboards because the failure is often silent. There is no exception when you compute area in a geographic CRS — you simply get a number in the wrong units. There is no exception when latitude and longitude are swapped — the point just lands somewhere impossible. This reference is organised to catch those failures before they reach a user, with a numbered workflow you can lift into a data-loading function and a troubleshooting section keyed to the exact errors and symptoms you will see.

Jump to heading# Problem statement

Spatial data arrives from many sources — PostGIS, GeoJSON exports, national mapping agencies, GPS devices — and each has its own convention. One feed is in WGS 84 degrees, another in a national grid in metres, a third has no declared CRS at all. Before any of it can be layered on a map, filtered by a bounding box, or measured, it must be reconciled to a known system. The dashboard then has to reproject again at the display boundary, because web tile maps expect Web Mercator. Every one of those steps is a place to introduce a mismatch, and because the dashboard re-runs its data path on each interaction, a mismatch is not a one-time cost but a recurring one.

Jump to heading# Prerequisites

• Python 3.10+ with geopandas>=0.14, shapely>=2.0, and pyproj>=3.5
• A working understanding of the GeoDataFrame object and its .crs attribute, covered in the GeoDataFrame schema reference
• Familiarity with how a filtered layer feeds a live map, described under dynamic spatial filtering
• The parent spatial data reference for how CRS sits alongside schema and indexing in the data layer

Jump to heading# Geographic vs projected: the distinction that drives everything

A geographic CRS places a point with latitude and longitude — angles measured on a reference ellipsoid. EPSG:4326 (WGS 84) is the archetype, and it is the right choice for canonical storage and interchange because GeoJSON, GPS, and most public feeds speak it. What it is not good for is measurement: a degree of longitude is about 111 km at the equator and shrinks to zero at the poles, so a Euclidean distance between two lat/lon pairs is meaningless.

A projected CRS flattens the ellipsoid onto a plane and reports coordinates in linear units, almost always metres. EPSG:3857 (Web Mercator) is the projected system every web map uses for display. National grids such as EPSG:27700 (British National Grid) and zone systems such as EPSG:32633 (UTM zone 33N) are projected systems tuned to minimise distortion over a specific region — they are what you reproject into when you need a real distance or area.

A third, important sub-family is the equal-area projected CRS. EPSG:6933 globally, EPSG:3035 for Europe, and EPSG:5070 for the contiguous US preserve area at the cost of shape and angle. These are the systems to reproject into before any area calculation, because Web Mercator inflates area dramatically away from the equator.

Jump to heading# Reference: EPSG codes by role

This is the table to keep open while building. The units and area of use columns are the ones that prevent mistakes: a projected code is only trustworthy for measurement inside its stated extent, and any calculation on a degrees-based code needs reprojection first.

The mnemonic that collapses this table into a habit: store in 4326, display in 3857, measure in a local metric or equal-area code.

Jump to heading# Core implementation workflow

Jump to heading# Step 1 — Inspect the CRS a frame already carries

Never assume. The first line of any load path should report what the source actually declared.

import geopandas as gpd

gdf = gpd.read_file("districts.geojson")
print(gdf.crs)                 # e.g. EPSG:4326
print(gdf.crs.axis_info)       # confirms axis names and order
print(gdf.total_bounds)        # sanity-check the coordinate ranges

If total_bounds shows values in the hundreds of thousands rather than the −180 to 180 range, the data is already projected — reading it as degrees would misplace everything.

Jump to heading# Step 2 — Set a CRS when the source declares none

Shapefiles missing a .prj and some Parquet exports arrive with gdf.crs is None. You cannot reproject an unknown source, so you must assert the CRS you know the data is in. Use set_crs, which attaches metadata without moving any coordinates — the opposite of to_crs.

if gdf.crs is None:
    gdf = gdf.set_crs("EPSG:4326")   # declare, do not transform

Confusing set_crs (label) with to_crs (transform) is the most common CRS bug: calling to_crs on unlabeled data raises, and calling set_crs on already-labeled data can silently mislabel it.

Jump to heading# Step 3 — Reproject with to_crs

Once the source CRS is known, to_crs moves every coordinate into the target system.

gdf_web = gdf.to_crs("EPSG:3857")      # for the tile map
gdf_metric = gdf.to_crs("EPSG:25832")  # for distance/area in central EU

to_crs returns a new frame; it does not mutate in place unless you assign back. Each call materialises a full copy of the coordinate arrays, which is a memory consideration on large frames — see memory limit management via the schema reference for chunking strategies.

Jump to heading# Step 4 — Choose the CRS that matches the operation

The single decision that prevents most CRS errors is matching the system to the task. Display, distance, and area each want a different CRS.

# distance between two points, correctly, in metres
metric = gdf.to_crs("EPSG:25832")
dist_m = metric.geometry.iloc[0].distance(metric.geometry.iloc[1])

# area, correctly, in an equal-area CRS
gdf["area_km2"] = gdf.to_crs("EPSG:3035").area / 1_000_000

Jump to heading# Step 5 — Verify axis order and units before trusting the result

After any reprojection, confirm the coordinate ranges and units are what you expect. A quick bounds check catches an axis swap immediately: valid Web Mercator eastings are roughly ±20 million metres, valid WGS 84 longitudes are ±180 degrees.

minx, miny, maxx, maxy = gdf_web.total_bounds
assert abs(maxx) < 2.1e7, "x out of Web Mercator range — check axis order"

Jump to heading# Advanced patterns

Jump to heading# Custom CRS from a PROJ string

When a dataset uses a projection with no EPSG code — a local mine grid, a historical survey — build the CRS from its PROJ definition and attach it explicitly.

from pyproj import CRS

custom = CRS.from_proj4(
    "+proj=laea +lat_0=52 +lon_0=10 +x_0=0 +y_0=0 "
    "+ellps=GRS80 +units=m +no_defs"
)
gdf_local = gdf.to_crs(custom)

Jump to heading# Area-weighted calculations in an equal-area CRS

Population density, land-cover fractions, and any per-area statistic must be computed in an equal-area system or the result skews with latitude. Reproject once, compute, and join the result back to the display frame.

ea = gdf.to_crs("EPSG:6933")
gdf["density"] = gdf["population"] / (ea.area / 1_000_000)   # people per km²

Jump to heading# Geodesic distance without projecting

For point-to-point distances that span many degrees — flight paths, long baselines — a single projected CRS distorts too much. Use pyproj’s Geod to measure on the ellipsoid directly.

from pyproj import Geod

geod = Geod(ellps="WGS84")
# Berlin (13.40, 52.52) to Rome (12.50, 41.90)
_, _, dist_m = geod.inv(13.40, 52.52, 12.50, 41.90)
print(f"{dist_m / 1000:.0f} km")   # ~1183 km

Jump to heading# Verification and testing

Round-tripping is the fastest correctness check: reproject out and back, and confirm coordinates return to within floating-point tolerance.

import geopandas as gpd

gdf = gpd.read_file("districts.geojson").to_crs("EPSG:4326")
roundtrip = gdf.to_crs("EPSG:3857").to_crs("EPSG:4326")

assert gdf.geometry.geom_equals_exact(roundtrip.geometry, tolerance=1e-6).all(), \
    "Round-trip drifted — check for a datum shift or clipping at high latitude"
assert gdf.crs.to_epsg() == 4326

For measurement code, assert the units are plausible: a UK district area should be on the order of tens to hundreds of square kilometres, never a fraction of one (the signature of measuring in degrees).

Jump to heading# Troubleshooting

pyproj.exceptions.CRSError: Invalid projection : The string passed to set_crs or to_crs is not a recognised CRS. Use the canonical form "EPSG:4326" rather than "WGS84" or "4326", or build the CRS with CRS.from_user_input(...) which accepts more formats.

ValueError: Cannot transform naive geometries. Please set a crs on the object first : You called to_crs on a frame whose .crs is None. Declare the known source CRS with set_crs first (Step 2), then reproject.

Coordinates land in the wrong hemisphere after transform : An axis-order swap. EPSG:4326 formally orders axes latitude-then-longitude, but GeoJSON and most tools expect longitude-then-latitude. Pass always_xy=True when building a pyproj.Transformer so it consistently treats input as (lon, lat).

Area or distance values are tiny fractions : You measured in a geographic CRS. gdf.area on a 4326 frame returns square degrees. Reproject to a projected or equal-area CRS (Step 4) before calling .area or .distance.

Layers align at the equator but drift toward the poles : You are treating Web Mercator as equal-area, or overlaying a 4326 layer on a 3857 base without reprojecting. Reproject the overlay to match the base and never measure area in 3857.

Jump to heading# Performance considerations

Reprojection is not free: each to_crs copies the full coordinate array, so on a large frame it is worth doing once at load time and caching the result rather than per interaction. If you only need a projected copy for measurement, drop it after computing the attribute column so two copies of the geometry do not linger in memory. When many small reprojections are unavoidable — for example per-feature geodesic distances — prefer a vectorised pyproj.Transformer over per-geometry Shapely transforms, which are an order of magnitude slower.

For the two tasks engineers repeat most, this section drills deeper: looking up EPSG codes for your region walks through finding the right code programmatically, and reprojecting between EPSG:4326 and EPSG:3857 covers the display-boundary conversion without distortion.

Jump to heading# Frequently asked questions

Back to Spatial Data Reference

Related

• Looking Up EPSG Codes for Your Region — find the right code with pyproj’s UTM and area-of-interest queries
• Reprojecting Between EPSG:4326 and EPSG:3857 — the display-boundary conversion without distortion or clipping errors
• Dynamic Spatial Filtering — how reprojected coordinates feed a live, filterable map
• Spatial Data Reference — the parent section linking CRS, schema, and indexing