EPSG Codes Explained: How to Choose the Right CRS in Python
On this page
Problem statement
EPSG:4326, EPSG:3857, EPSG:32633: an EPSG code is a number in a public registry that stands for a complete coordinate reference system — datum, ellipsoid, projection, units, axis order. Code reads and writes them constantly, so two practical questions come up again and again: what does a given code actually define, and which code should a dataset for a given place be in?
This page answers both for Vienna. It asks the EPSG registry, through pyproj, which projected systems cover the city, and then measures the city's 23 districts in four of them against the city's own published areas — one of which reproduces the official figures to the square metre.
Quick answer
An EPSG code is a key into the EPSG registry; pyproj.CRS.from_epsg(code) gives the full definition behind it, and gdf.crs.to_epsg() goes the other way. To choose one, ask the registry what covers your area, then prefer the system your data's publisher uses: measured on Vienna, that one matched the official areas exactly, a UTM zone was 0.06% off, and Web Mercator more than doubled them.
Check the code, read its definition, ask the registry, measure
Step 1: Check the current CRS of your data
Always inspect the CRS first.
import geopandas as gpd
gdf = gpd.read_file("data/city_parks.geojson")
print(gdf.crs)EPSG:4326
A code, a full WKT definition with no code, or None — three different situations, the last of which is cannot transform naive geometries.
Step 2: Read what the code defines
The code is a short name; the definition behind it is what PROJ uses:
import pyproj
for code in (4326, 3857, 32633, 31256):
crs = pyproj.CRS.from_epsg(code)
axes = ", ".join(f"{a.abbrev} {a.direction} ({a.unit_name})" for a in crs.axis_info)
print(f"EPSG:{code:<6} {crs.name:28} {crs.datum.name:42} {axes}")EPSG:4326 WGS 84 World Geodetic System 1984 ensemble Lat north (degree), Lon east (degree)
EPSG:3857 WGS 84 / Pseudo-Mercator World Geodetic System 1984 ensemble X east (metre), Y north (metre)
EPSG:32633 WGS 84 / UTM zone 33N World Geodetic System 1984 ensemble E east (metre), N north (metre)
EPSG:31256 MGI / Austria GK East Militar-Geographische Institut X north (metre), Y east (metre)
Two of the four put north first. EPSG:4326 is officially latitude then longitude, and MGI / Austria GK East calls its northing X and its easting Y. GeoPandas and Shapely always treat x as east, so this rarely shows in Python; a tool that follows the registry literally, such as a WMS 1.3 server, does it the registry's way. The datum line matters too: three of the four share WGS 84, and 31256 sits on the Austrian MGI datum, so converting to or from it is a datum shift as well as a projection.
Step 3: Ask the registry what covers your area
The EPSG database knows every system's area of use, so the candidates for a place can be listed rather than remembered:
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_crs_info
from pyproj.enums import PJType
districts = gpd.read_file("data/districts.geojson")
w, s, e, n = districts.total_bounds
found = query_crs_info(auth_name="EPSG", pj_types=PJType.PROJECTED_CRS, contains=True,
area_of_interest=AreaOfInterest(w, s, e, n))
print(len(found), "projected EPSG systems cover all of Vienna")
austrian = [c for c in found if "Austria" in c.name]
print(len(austrian), "of them Austrian, e.g.", [f"{c.code} {c.name}" for c in austrian if c.code in ("31256", "3416")])
print([f"{c.code} {c.name}" for c in found if c.code in ("32633", "3035", "3857")])50 projected EPSG systems cover all of Vienna
9 of them Austrian, e.g. ['3416 ETRS89 / Austria Lambert', '31256 MGI / Austria GK East']
['3035 ETRS89-extended / LAEA Europe', '3857 WGS 84 / Pseudo-Mercator', '32633 WGS 84 / UTM zone 33N']
Fifty projected systems in the registry cover the whole city. Nine are Austrian, among them Austria Lambert (EPSG:3416) and the national grid the city itself publishes in, MGI / Austria GK East (EPSG:31256); the rest are wider systems such as LAEA Europe, UTM zone 33N and Web Mercator, whose area of use is the world. That is a list of candidates, not an answer; Example 5 measures four of them against the city's own figures.
Step 4: Reproject data when needed
Measure the districts in UTM zone 33N against the city's own area column, FLAECHE:
import geopandas as gpd
districts = gpd.read_file("data/districts.geojson")
print(districts.crs)
districts_m = districts.to_crs("EPSG:32633")
districts_m["area_sq_m"] = districts_m.area
print(districts_m[["NAMEK", "area_sq_m", "FLAECHE"]].head(3).round(0).to_string(index=False))EPSG:4326
NAMEK area_sq_m FLAECHE
Innere Stadt 2867180.0 2868772.0
Simmering 23243928.0 23256174.0
Neubau 1607341.0 1608246.0
FLAECHE is the city's own figure for each district. UTM zone 33N lands within a few hundredths of a percent of it. If you calculate area directly in EPSG:4326, the results are in square degrees — all 23 districts of Vienna add up to 0.05 — and in Web Mercator they come out 2.25 times too large (Example 5).
Step 5: Verify the result before analysis
After assigning or transforming CRS, check that the result makes sense.
import geopandas as gpd
roads = gpd.read_file("data/roads.shp")
parcels = gpd.read_file("data/census_districts.shp")
print(roads.crs)
print(parcels.crs)
parcels = parcels.to_crs(roads.crs)
print(roads.total_bounds.round(4))
print(parcels.total_bounds.round(4))EPSG:4326
EPSG:31256
[16.1874 48.121 16.5501 48.3214]
[16.1818 48.1179 16.5775 48.3227]
The bounds of both layers now fall around 16.2–16.6°E and 48.1–48.3°N, which is Vienna.
Vienna's districts measured in four CRSs
The data used on this page
Vienna, from the city's open data WFS: the 23 district boundaries with the city's own area for each, the 1,054 park locations, the 2,308 public drinking fountains, the 29,658 edges of the street graph, and the 250 census districts (Zählbezirke). Vienna sits in UTM zone 33N, which is why EPSG:32633 appears throughout. The census districts are saved in MGI / Austria GK East (EPSG:31256), the national CRS the city itself uses; the fountains are also saved once as a shapefile with its .prj deleted.
import os
import geopandas as gpd
import requests
WFS = "https://data.wien.gv.at/daten/geo"
def wien(layer):
reply = requests.get(WFS, timeout=180, params={
"service": "WFS", "request": "GetFeature", "version": "1.1.0",
"typeName": f"ogdwien:{layer}", "srsName": "EPSG:4326", "outputFormat": "json"})
reply.raise_for_status()
return gpd.read_file(reply.text)
os.makedirs("data", exist_ok=True)
if not os.path.exists("data/districts.geojson"):
wien("BEZIRKSGRENZEOGD")[["NAMEK", "BEZNR", "FLAECHE", "geometry"]].to_file("data/districts.geojson")
wien("PARKINFOOGD")[["ANL_NAME", "BEZIRK", "geometry"]].to_file("data/city_parks.geojson")
fountains = wien("TRINKBRUNNENOGD")[["BASIS_TYP_TXT", "geometry"]]
fountains.to_file("data/gps_points.geojson")
fountains.to_file("data/fountains_noprj.shp")
os.remove("data/fountains_noprj.prj")
wien("STRASSENGRAPHOGD")[["FEATURENAME", "EDGECATEGORY", "geometry"]].to_file("data/roads.shp")
wien("ZAEHLBEZIRKOGD")[["ZBEZ", "BEZNR", "FLAECHE", "geometry"]].to_crs(31256).to_file(
"data/census_districts.shp")
for name in sorted(os.listdir("data")):
if name.endswith((".geojson", ".shp")):
print(f"{name:24} {len(gpd.read_file('data/' + name)):>6,} features")census_districts.shp 250 features
city_parks.geojson 1,054 features
districts.geojson 23 features
fountains_noprj.shp 2,308 features
gps_points.geojson 2,308 features
roads.shp 29,658 features
Example 1: Read a file and inspect its CRS
import geopandas as gpd
buildings = gpd.read_file("data/census_districts.shp")
print("CRS:", buildings.crs, "-", buildings.crs.name)CRS: EPSG:31256 - MGI / Austria GK East
This confirms whether CRS metadata exists — here a national projected CRS rather than a UTM zone or EPSG:4326.
Example 2: Assign a missing EPSG code to data
import geopandas as gpd
points = gpd.read_file("data/fountains_noprj.shp")
print("as read:", points.crs, points.total_bounds.round(3))
if points.crs is None:
points = points.set_crs("EPSG:4326")
print(points.crs)
print("the GeoJSON copy:", gpd.read_file("data/gps_points.geojson").crs)as read: None [16.184 48.123 16.545 48.317]
EPSG:4326
the GeoJSON copy: EPSG:4326
Use this only when the coordinates are already in that CRS and the metadata is missing — here the bounds are plainly degrees around Vienna, 16.2–16.6°E and 48.1–48.3°N. The GeoJSON copy never needs it: GeoJSON is defined as WGS 84, so GDAL gives it EPSG:4326 whatever the file says, and a GeoJSON written in another CRS without its legacy crs member comes back wrongly labelled rather than unlabelled.
Example 3: Reproject data to a projected CRS for analysis
import geopandas as gpd
points = gpd.read_file("data/gps_points.geojson")
if points.crs is None:
points = points.set_crs("EPSG:4326")
points_utm = points.to_crs("EPSG:32633")
points_utm["buffer_500m"] = points_utm.buffer(500)
print(points_utm.crs, f"one buffer: {points_utm['buffer_500m'].iloc[0].area:,.0f} m²")EPSG:32633 one buffer: 784,137 m²
A circle of radius 500 m has an area of 785,398 m²; the buffer is a 64-sided polygon, so it comes out slightly smaller. The same buffer(500) in EPSG:4326 would be a 500-degree circle.
This is the correct pattern when you need meter-based analysis.
Example 4: Confirm layer alignment after reprojection
import geopandas as gpd
rivers = gpd.read_file("data/roads.shp")
catchments = gpd.read_file("data/census_districts.shp")
print("before:", catchments.total_bounds.round(0))
catchments_aligned = catchments.to_crs(rivers.crs)
print("Roads bounds: ", rivers.total_bounds.round(4))
print("Census districts bounds: ", catchments_aligned.total_bounds.round(4))
joined = gpd.sjoin(rivers, catchments_aligned, predicate="intersects")
print(f"{joined.index.nunique():,} of {len(rivers):,} road edges touch a census district")before: [-11180. 331051. 18258. 353820.]
Roads bounds: [16.1874 48.121 16.5501 48.3214]
Census districts bounds: [16.1818 48.1179 16.5775 48.3227]
29,650 of 29,658 road edges touch a census district
If the bounds now fall in the same region and a join finds matches, the reprojection likely worked. Before it, the census districts' bounds are in metres east and north of the Austrian grid origin, which no amount of looking at a map in degrees would line up with.
Example 5: The same districts measured in four CRSs
import geopandas as gpd
districts = gpd.read_file("data/districts.geojson")
official = districts["FLAECHE"].sum()
print(f"official area of Vienna's 23 districts: {official:,.0f} m²")
for code in ("EPSG:4326", "EPSG:3857", "EPSG:32633", "EPSG:31256"):
area = districts.to_crs(code).area
worst = (area / districts["FLAECHE"] - 1).abs().max()
print(f"{code:10} total {area.sum():>16,.4f} {area.sum() / official:8.5f} × official"
+ (f" worst district off by {worst:.3%}" if code != "EPSG:4326" else " (square degrees)"))
print("estimate_utm_crs():", districts.estimate_utm_crs())official area of Vienna's 23 districts: 414,871,121 m²
EPSG:4326 total 0.0502 0.00000 × official (square degrees)
EPSG:3857 total 933,444,696.6628 2.24996 × official worst district off by 125.633%
EPSG:32633 total 414,642,193.6182 0.99945 × official worst district off by 0.060%
EPSG:31256 total 414,871,121.5160 1.00000 × official worst district off by 0.000%
estimate_utm_crs(): EPSG:32633
EPSG:31256 reproduces the city's figures to the square metre because it is the CRS they were computed in. UTM zone 33N is within 0.06% for every district — the scale factor of a transverse Mercator zone, not an error. Web Mercator more than doubles every area at Vienna's latitude.
Why the publisher's code wins
The EPSG registry, maintained by the oil and gas industry's geomatics committee (IOGP), assigns a number to each datum, ellipsoid, projection and complete CRS, and PROJ ships a copy of it. A code is therefore not a description but a lookup: two files that both say EPSG:31256 agree on every parameter, and a file that carries only WKT may describe the same system without saying so — to_epsg() recovers the code when it can.
For choosing among the candidates, Example 5 is the argument. Vienna's district areas in FLAECHE were computed by the city in MGI / Austria GK East, EPSG:31256, and in that system GeoPandas reproduces them to the square metre. UTM zone 33N, a perfectly good metric system, is 0.06% off for every district: not an error, the zone's scale factor. Web Mercator is 2.25 times too large. When the data's publisher used a national system, measuring in it is the only way to match their numbers exactly; when no such system applies, a UTM zone is the reasonable default (choosing a projected CRS).
WKT-only files, ESRI codes and Web Mercator
- A file with WKT and no code is not missing a CRS.
gdf.crs.to_epsg()finds the matching code when one exists;to_epsg(min_confidence=…)controls how close the match must be. ESRI:codes are a separate registry;pyproj.CRS("ESRI:102003")reads them, and they have no EPSG number.- EPSG:3857 has a world-wide area of use and appears for every query; it is for web tiles, not measurement.
Internal links
To understand the broader concept, see Coordinate Reference Systems (CRS) Explained for Python GIS.
For related tasks, read GeoPandas Basics: Working with Spatial Data in Python and Python for GIS: What It Is and When to Use It.
FAQ
What is the difference between an EPSG code and a CRS?
A CRS is the full coordinate reference system definition. An EPSG code is a standard identifier for one specific CRS.
When should I use EPSG:4326 in Python GIS?
Use EPSG:4326 for latitude/longitude data, GPS data, GeoJSON, and many API or storage workflows. It is usually not the best choice for distance or area analysis.
How do I find the right UTM EPSG code for my area?
Pick the UTM zone that covers your data's longitude and hemisphere; for example EPSG:32633 is UTM zone 33N (northern hemisphere) and EPSG:32733 is the matching zone 33S. You can also let GeoPandas suggest one with gdf.estimate_utm_crs(), which returns the appropriate UTM CRS for the layer's extent.
What does EPSG:3857 actually mean and when should I avoid it?
EPSG:3857 is Web Mercator, a projected CRS in metres used by most web basemaps like OpenStreetMap and Google Maps. Avoid it for distance and area work, because its scale distortion grows sharply toward the poles, so measurements are unreliable away from the equator.
Which EPSG code should I use for the US or the UK?
For continental US analysis, EPSG:5070 (NAD83 CONUS Albers, metres) is a good equal-area choice for national distance and area work. For Great Britain, use EPSG:27700 (British National Grid, metres), the standard projected CRS for UK data.
Does changing the EPSG code with set_crs move my data on the map?
No. set_crs() only attaches metadata and assumes the coordinates already match that CRS, so it never moves the geometry. If you assign the wrong code with set_crs(), the data will plot in the wrong place; use to_crs() when you actually want to transform coordinates.