Coordinate Reference Systems (CRS)

Coordinate-system reference, reviewed 2026-10-01 against PyProj/PROJ. Code with external data is illustrative. CRS names were checked against the installed PROJ database; a numerical round trip does not establish datum accuracy or grid availability.

Table of Contents

1. Fundamentals
2. Common CRS Codes
3. Projected vs Geographic
4. UTM Zones
5. Transformations
6. Best Practices

Fundamentals

What is a CRS?

A Coordinate Reference System defines how coordinates relate to positions on Earth:

• Geographic CRS: Uses angular latitude/longitude; usually degrees, but inspect axis units
• Projected CRS: Uses Cartesian coordinates (meters, feet)
• Vertical CRS: Defines height/depth (e.g., ellipsoidal heights)

Components

1. Datum: Mathematical model of Earth's shape
2. WGS 84 (EPSG:4326) - Global GPS
3. NAD 83 (EPSG:4269) - North America
4. ETRS89 (EPSG:4258) - Europe

1. Projection: Transformation from curved to flat surface
2. Cylindrical (Mercator)
3. Conic (Lambert Conformal)
4. Azimuthal (Polar Stereographic)

1. Units: Degrees, meters, feet, etc.

Common CRS Codes

Geographic CRS (Lat/Lon)

Projected CRS (Meters)

Regional Projections

United States:

• EPSG:5070 - NAD83 / Conus Albers
• EPSG:3338 - NAD83 / Alaska Albers
• ESRI:102003 - USA Contiguous Albers Equal Area
• EPSG:2227 - California Zone 3 (US Feet)

Europe:

• EPSG:3035 - Europe Equal Area 2001
• EPSG:3857 - Web Mercator (web mapping)
• EPSG:2154 - Lambert 93 (France)
• EPSG:25832-25836 - UTM zones (ETRS89)

Other:

• EPSG:3112 - GDA94 / Geoscience Australia Lambert (Australia)
• EPSG:2056 - CH1903+ / LV95 (Switzerland)
• EPSG:4326 - WGS 84 (global default)

Projected vs Geographic

When to Use Geographic (EPSG:4326)

✅ Storing data (databases, files) ✅ Global datasets ✅ Web APIs (GeoJSON, KML) ✅ Latitude/longitude queries ✅ GPS coordinates

# Bad: Distance calculation in geographic CRS
# Assume gdf already has a verified geographic CRS
distance = gdf.geometry.length  # WRONG! Returns degrees, not meters

# Good: Calculate distance in projected CRS
gdf_projected = gdf.to_crs("EPSG:32633")  # UTM Zone 33N
distance_m = gdf_projected.geometry.length  # Correct: meters

When to Use Projected

✅ Area/distance calculations ✅ Buffer operations ✅ Spatial analysis ✅ High-resolution mapping ✅ Engineering applications

import geopandas as gpd

# Project to appropriate UTM zone
gdf = gdf.to_crs(gdf.estimate_utm_crs())

# Now area and distance are accurate
area_sqm = gdf.geometry.area
buffer_1km = gdf.geometry.buffer(1000)  # 1000 meters

Web Mercator Warning

⚠️ EPSG:3857 (Web Mercator) for visualization only

# DON'T use Web Mercator for area calculations
gdf_web = gdf.to_crs("EPSG:3857")
area = gdf_web.geometry.area  # WRONG! Significant distortion

# DO use appropriate projection
gdf_utm = gdf.to_crs("EPSG:32633")  # or estimate_utm_crs()
area = gdf_utm.geometry.area  # Correct

UTM Zones

Understanding UTM Zones

Earth is divided into 60 zones (6° longitude each):

• Zones 1-60: West to East
• Each zone divided into North (326xx) and South (327xx)

Finding Your UTM Zone

def get_utm_zone(longitude, latitude):
    """Nominal 6-degree zone only; special Norway/Svalbard zones need registry lookup."""
    import math

    if not (-180 <= longitude <= 180 and -80 <= latitude <= 84):
        raise ValueError('Outside the nominal UTM domain')
    zone = min(60, math.floor((longitude + 180) / 6) + 1)

    if latitude >= 0:
        epsg = 32600 + zone  # Northern hemisphere
    else:
        epsg = 32700 + zone  # Southern hemisphere

    return f"EPSG:{epsg}"

# Example
get_utm_zone(-122.4, 37.7)  # Returns 'EPSG:32610' (Zone 10N)

Auto-Detect UTM Zone with GeoPandas

import geopandas as gpd

# Load data
gdf = gpd.read_file('data.geojson')

# Estimate best UTM zone
utm_crs = gdf.estimate_utm_crs()
print(f"Best UTM CRS: {utm_crs}")

# Reproject
gdf_projected = gdf.to_crs(utm_crs)

Special UTM Cases

UPS (Universal Polar Stereographic):

• EPSG:5041 - UPS North (Arctic)
• EPSG:5042 - UPS South (Antarctic)

Other national grids (not UTM):

• EPSG:31466-31469 - German Gauss-Krüger zones
• EPSG:2056 - Swiss LV95 (oblique Mercator)

Transformations

Basic Transformation

from pyproj import Transformer

# Create transformer
transformer = Transformer.from_crs(
    "EPSG:4326",  # WGS 84 (lat/lon)
    "EPSG:32610", # UTM Zone 10N for this California point
    always_xy=True  # Input: x=lon, y=lat (not y=lat, x=lon)
)

# Transform single point
lon, lat = -122.4, 37.7
x, y = transformer.transform(lon, lat)
print(f"Easting: {x:.2f}, Northing: {y:.2f}")

Batch Transformation

import numpy as np
from pyproj import Transformer

# Arrays of coordinates
lon_array = [-122.4, -122.3]
lat_array = [37.7, 37.8]

transformer = Transformer.from_crs("EPSG:4326", "EPSG:32610", always_xy=True)
xs, ys = transformer.transform(lon_array, lat_array)

Transformation with PyProj CRS

from pyproj import CRS

# Get CRS information
crs = CRS.from_epsg(32633)

print(f"Name: {crs.name}")
print(f"Type: {crs.type_name}")
print(f"Area of use: {crs.area_of_use.name}")
print(f"Datum: {crs.datum.name}")
print(f"Ellipsoid: {crs.ellipsoid.name}")

Best Practices

1. Always Know Your CRS

import geopandas as gpd

gdf = gpd.read_file('data.geojson')

# Check CRS immediately
print(f"CRS: {gdf.crs}")  # Should never be None!

# If absent, obtain the source CRS from metadata; never guess it from bounds.
if gdf.crs is None:
    raise ValueError("Source CRS is unknown")

2. Verify CRS Before Operations

def ensure_same_crs(gdf1, gdf2):
    """Ensure two GeoDataFrames have same CRS."""
    if gdf1.crs is None or gdf2.crs is None:
        raise ValueError("Source CRS is unknown")
    if gdf1.crs != gdf2.crs:
        gdf2 = gdf2.to_crs(gdf1.crs)
        print(f"Reprojected gdf2 to {gdf1.crs}")
    return gdf1, gdf2

# Use before spatial operations
zones, points = ensure_same_crs(zones_gdf, points_gdf)
result = gpd.sjoin(points, zones, predicate='within')

3. Use Appropriate Projections

# For local analysis (< 500km extent)
gdf_local = gdf.to_crs(gdf.estimate_utm_crs())

# For national/regional analysis
gdf_us = gdf.to_crs("EPSG:5070")  # NAD83 / Conus Albers
gdf_eu = gdf.to_crs("EPSG:3035")  # Europe Equal Area

# For web visualization
gdf_web = gdf.to_crs("EPSG:3857")  # Web Mercator

4. Preserve Original CRS

# Keep original as backup
gdf_original = gdf.copy()
original_crs = gdf.crs

# Do analysis in projected CRS
gdf_projected = gdf.to_crs(gdf.estimate_utm_crs())
result = gdf_projected.geometry.buffer(1000)

# Convert back if needed
result = result.to_crs(original_crs)

Common Pitfalls

Mistake 1: Area in Degrees

# WRONG: Area in square degrees
gdf = gpd.read_file('data.geojson')
area = gdf.geometry.area  # Wrong!

# CORRECT: Use projected CRS
gdf_proj = gdf.to_crs(gdf.estimate_utm_crs())
area_sqm = gdf_proj.geometry.area
area_sqkm = area_sqm / 1_000_000

Mistake 2: Buffer in Geographic CRS

# WRONG: Buffer of 1000 degrees
gdf['buffer'] = gdf.geometry.buffer(1000)

# CORRECT: Project first
gdf_proj = gdf.to_crs("EPSG:32610")
gdf_proj['buffer_km'] = gdf_proj.geometry.buffer(1000)  # 1000 meters

Mistake 3: Mixing CRS

# WRONG: Spatial join without checking CRS
result = gpd.sjoin(gdf1, gdf2, predicate='intersects')

# CORRECT: Ensure same CRS
if gdf1.crs != gdf2.crs:
    gdf2 = gdf2.to_crs(gdf1.crs)
result = gpd.sjoin(gdf1, gdf2, predicate='intersects')

Quick Reference

# Common operations

# Check CRS
gdf.crs
rasterio.open('file.tif').crs

# Reproject
gdf.to_crs("EPSG:32633")

# Auto-detect UTM
gdf.estimate_utm_crs()

# Transform single point
from pyproj import Transformer
tx = Transformer.from_crs("EPSG:4326", "EPSG:32610", always_xy=True)
x, y = tx.transform(lon, lat)

# Create custom CRS
from pyproj import CRS
custom_crs = CRS.from_proj4(
    "+proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
)

For more information, see:

• EPSG Registry
• PROJ documentation
• pyproj documentation