Coordinate Reference Systems (CRS) Explained for Python GIS

On this page

Problem statement

Coordinates are just numbers until something says what they measure: from which origin, along which axes, in which unit, on which model of the earth. That something is the coordinate reference system, and it is the reason a layer can load without an error and still be wrong — in the wrong place, in the wrong units, or measured in degrees as though they were metres.

This page explains what a CRS is, what its parts decide, and how GeoPandas carries one, using New York City's own layers: some arrive in longitude and latitude, and some in the city's State Plane system, measured in US survey feet.

Quick answer

A CRS tells GIS software how coordinate numbers relate to real locations on the earth.

In practice, there are two separate questions:

What numbers are stored in the geometry?
Example: (-73.9857, 40.7484) or (988212, 211939) — the Empire State Building in longitude/latitude and in New York's State Plane feet
What do those numbers mean?
Are they longitude/latitude in degrees, or projected coordinates in meters?

In GeoPandas the answer lives in gdf.crs. Check it on every layer, label a missing one only from what the source documents (set_crs), and convert with to_crs — the difference between those two is set_crs vs to_crs.

What a CRS does, geographic or projected, assign or reproject

An EPSG code is shorthand for all of this — which is why one number is enough to be precise.

1. Understand what a CRS does

A CRS defines how coordinates map to real-world locations. The same geometry values can represent completely different places depending on the CRS.

For example, these numbers:

POINT (12 55)

could mean:

longitude 12°, latitude 55° in a geographic CRS
x=12 meters, y=55 meters in a projected CRS

Without CRS information, the coordinates are just numbers.

2. Know the two kinds, and what a code names

A geographic CRS gives positions as angles — longitude and latitude in degrees, such as EPSG:4326, WGS 84. A projected CRS flattens the earth onto a plane and gives distances from an origin, in metres or feet, such as EPSG:32618, UTM zone 18N, or the city's EPSG:2263. Angles locate; only a plane measures — why, in detail, is projected vs geographic CRS. Every one of these has an EPSG code, a key into a public registry that holds the full definition (EPSG codes explained).

3. Inspect CRS before doing any analysis

Use the .crs attribute on a GeoDataFrame.

import geopandas as gpd

gdf = gpd.read_file("data/cities.geojson")
print(gdf.crs)
print(gdf.crs == "EPSG:4326")
print(gpd.read_file("data/roads.shp").crs.name, "|", gpd.read_file("data/roads.shp").crs.axis_info[0].unit_name)

EPSG:4326
True
NAD83 / New York Long Island (ftUS) | US survey foot

A GeoJSON reports EPSG:4326; the city's own street layer, as a shapefile, reports New York's State Plane — in US survey feet, which is the unit every distance and area from it will come out in. A missing CRS appears as None.

4. Put layers in one CRS before they meet

overlay, sjoin, clip and every measurement compare coordinates directly, so two layers in different systems are compared number to number — State Plane feet against degrees — and match nothing. Convert one to the other's crs first; which system to measure in is choosing a projected CRS.

print(gdf1.crs == gdf2.crs)

New York City layers inspected, assigned and aligned

The data used on this page

New York City's own layers. The Department of City Planning publishes, through its ArcGIS services, the borough boundaries, the community districts, the arterial and major street network, and the future floodplain for a 100-year storm in the 2080s (NPCC projections); NYC Open Data publishes the building footprints. The setup writes them in the formats and CRSs the examples read — including the streets and the footprints in New York's own State Plane Long Island zone (EPSG:2263, US survey feet), the CRS most city data comes in. The footprints and floodplain are for the southern tip of Manhattan.

from pathlib import Path

import geopandas as gpd
import pandas as pd
import requests

Path("data").mkdir(exist_ok=True)
DCP = "https://services5.arcgis.com/GfwWNkhOj9bNBqoJ/arcgis/rest/services/{}/FeatureServer/0/query?where=1%3D1&outFields=*&outSR={}&f=geojson"


def dcp(layer, srid):
    # the service returns its features a page at a time; follow exceededTransferLimit
    # with resultOffset rather than leaving the paging to GDAL
    pages, offset = [], 0
    while True:
        reply = requests.get(DCP.format(layer, srid) + f"&resultOffset={offset}", timeout=300)
        reply.raise_for_status()
        page = reply.json()
        pages.append(gpd.GeoDataFrame.from_features(page["features"], crs=srid))
        if not page.get("properties", {}).get("exceededTransferLimit") and not page.get("exceededTransferLimit"):
            return pd.concat(pages, ignore_index=True)
        offset += len(page["features"])


dcp("NYC_Borough_Boundary", 4326).to_file("data/cities.geojson")
dcp("DCM_Arterials_and_Major_Streets", 2263)[["Route_Name", "Borough", "geometry"]].to_file("data/roads.shp")
dcp("NYC_Community_Districts", 4326).to_file("data/parcels.geojson")
dcp("Future_Floodplains_2080_100yr_NPCC", 4326).to_file("data/flood_zones.geojson")
footprints = requests.get("https://data.cityofnewyork.us/resource/5zhs-2jue.geojson", timeout=120, params={
    "$where": "within_box(the_geom, 40.715, -74.020, 40.700, -74.000)", "$limit": 5000})
gpd.read_file(footprints.text)[["name", "construction_year", "geometry"]].to_crs(2263).to_file("data/buildings.shp")
for f in ("cities.geojson", "roads.shp", "parcels.geojson", "buildings.shp", "flood_zones.geojson"):
    g = gpd.read_file(f"data/{f}")
    print(f"{f:20} {len(g):5} features  EPSG:{g.crs.to_epsg()}")

cities.geojson           5 features  EPSG:4326
roads.shp              740 features  EPSG:2263
parcels.geojson         71 features  EPSG:4326
buildings.shp          747 features  EPSG:2263
flood_zones.geojson      1 features  EPSG:4326

Example 1: Read a file and inspect its CRS

import geopandas as gpd

roads = gpd.read_file("data/roads.shp")
print("CRS:", roads.crs)

if roads.crs is None:
    print("This layer has no CRS metadata.")
print(f"{roads.length.sum():,.0f} {roads.crs.axis_info[0].unit_name} of arterial street")

CRS: EPSG:2263
5,883,233 US survey foot of arterial street

The unit is part of the CRS: every length from this layer is in US survey feet.

Example 2: Assign a CRS to data that has none

import geopandas as gpd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame(
    {"name": ["Empire State Building", "One World Trade Center"]},
    geometry=[Point(-73.9857, 40.7484), Point(-74.0134, 40.7127)]
)

print(gdf.crs)  # None

# Coordinates are known to be longitude/latitude in WGS84
gdf = gdf.set_crs("EPSG:4326")

print(gdf.crs)

None
EPSG:4326

This is correct only if you already know the coordinates are in WGS 84. If you guess, you may attach the wrong CRS to valid geometry.

Example 3: Reproject a layer to a projected CRS

import geopandas as gpd

parcels = gpd.read_file("data/parcels.geojson")
print("Original CRS:", parcels.crs)

parcels_utm = parcels.to_crs("EPSG:32618")
print("Projected CRS:", parcels_utm.crs)

print([round(float(v), 4) for v in parcels.total_bounds])
print([round(float(v)) for v in parcels_utm.total_bounds])

Original CRS: EPSG:4326
Projected CRS: EPSG:32618
[-74.1615, 40.5362, -73.7, 40.9155]
[570894, 4487639, 609762, 4529953]

The coordinates are now metres east and north in UTM zone 18N, the zone New York lies in.

The same distance in all three, from Example 5.

Example 4: Align two layers to the same CRS before overlay

import geopandas as gpd

buildings = gpd.read_file("data/buildings.shp")
flood_zones = gpd.read_file("data/flood_zones.geojson")

print("Buildings CRS:", buildings.crs)
print("Flood zones CRS:", flood_zones.crs)

if buildings.crs != flood_zones.crs:
    flood_zones = flood_zones.to_crs(buildings.crs)

print("Aligned:", buildings.crs == flood_zones.crs)

import shapely

raw = flood_zones.geometry.iloc[0]
print(f"floodplain: one {raw.geom_type} of {shapely.get_num_geometries(raw):,} parts, valid: {raw.is_valid}")
try:
    shapely.intersection(raw, buildings.union_all())
except shapely.errors.GEOSException as err:
    print("intersection before repair:", str(err).split(" at ")[0])

flood_zones["geometry"] = flood_zones.geometry.make_valid(method="structure")   # the published polygon is invalid
at_risk = gpd.sjoin(buildings, flood_zones, predicate="intersects").index.nunique()
print(f"{at_risk} of {len(buildings)} buildings touch the 2080 100-year floodplain")

Buildings CRS: EPSG:2263
Flood zones CRS: EPSG:4326
Aligned: True
floodplain: one MultiPolygon of 2,570 parts, valid: False
intersection before repair: TopologyException: side location conflict
740 of 747 buildings touch the 2080 100-year floodplain

The buildings arrive in State Plane feet and the floodplain in longitude/latitude; joined as they are, nothing would match. The floodplain also arrives invalid — one multipolygon in thousands of parts — and an intersection with it fails until make_valid repairs it.

Example 5: Compare distance results in geographic vs projected CRS

import geopandas as gpd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame(
    {"name": ["Empire State Building", "One World Trade Center"]},
    geometry=[Point(-73.9857, 40.7484), Point(-74.0134, 40.7127)],
    crs="EPSG:4326"
)

# Distance in degrees, usually not useful for analysis
distance_degrees = gdf.geometry.iloc[0].distance(gdf.geometry.iloc[1])
print("Distance in geographic CRS:", distance_degrees)

# Reproject to a local projected CRS in meters
gdf_proj = gdf.to_crs("EPSG:32618")
distance_meters = gdf_proj.geometry.iloc[0].distance(gdf_proj.geometry.iloc[1])
print("Distance in projected CRS (meters):", distance_meters)

gdf_ft = gdf.to_crs("EPSG:2263")
print("Distance in State Plane (US feet):", gdf_ft.geometry.iloc[0].distance(gdf_ft.geometry.iloc[1]))

Distance in geographic CRS: 0.04518605979724776
Distance in projected CRS (meters): 4602.085478083383
Distance in State Plane (US feet): 15103.34588461151

One distance, three numbers: degrees, which are not a length; metres in UTM; and US survey feet in State Plane. The last two describe the same 4.6 km.

set_crs vs to_crs, decided

Inspect the CRS, then label it with set_crs() or transform it with to_crs().

What the parts of a CRS decide

Nearly every wrong-number bug in GIS traces back to measuring in the left-hand column.

New York's street layer shows every part of a CRS at once. Its code, EPSG:2263, names NAD83 / New York Long Island (ftUS): a datum, NAD83, which ties the coordinates to the North American plate rather than to the WGS 84 that GPS reports; a projection, the State Plane zone designed for Long Island and the five boroughs and for nowhere else; and a unit, the US survey foot. GeoPandas does arithmetic in that unit without comment, which is why the arterial network measures 5,883,233 — feet, not metres — and why the same distance in Example 5 is 15,103 in State Plane and 4,602 in UTM.

None of those parts is visible in the coordinates themselves. Two layers of the same buildings, one in State Plane feet and one in longitude and latitude, share no numbers at all, and a join between them as they arrived would have matched none of the 740 buildings that touch the 2080s floodplain once both are in one system.

An invalid layer with a correct CRS

A right CRS does not make a layer usable. New York's 2080s floodplain arrives as one invalid multipolygon (Example 4 counts its parts); its label is correct, the alignment works, and a direct intersection still raises TopologyException: side location conflict until the geometry is repaired with make_valid() — a separate problem, covered in fixing invalid geometries.

Internal links

FAQ

What is a CRS in GIS?

A CRS is the definition that tells GIS software how coordinate values relate to real positions on the earth. Without it, coordinates have no reliable spatial meaning.

What is the difference between EPSG:4326 and a projected CRS?

EPSG:4326 is a geographic CRS using longitude and latitude in degrees. A projected CRS uses planar units such as meters or feet, which makes it more suitable for distance, area, and buffering.

Why are my distance or area results wrong in Python GIS?

This usually happens because you measured data in a geographic CRS instead of a projected CRS, or because the layer had the wrong CRS assigned.

How do I find the right UTM zone or projected CRS for my data?

Use gdf.estimate_utm_crs(), which inspects the data's extent and returns the matching UTM CRS (for example UTM zone 33N is EPSG:32633). For national projects, a published regional CRS such as EPSG:5070 (NAD83 CONUS Albers) is often a better choice than UTM.

What does `EPSG:3857` mean and should I use it for analysis?

EPSG:3857 is Web Mercator, a projected CRS in metres used by most web basemaps. Its distances and areas are heavily distorted away from the equator, so use it for display only and reproject to a local CRS like a UTM zone before measuring.