Coordinate Reference Systems (CRS)

Complete guide to coordinate systems, projections, and transformations for geospatial data.

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 latitude/longitude (degrees)
• 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

   • WGS 84 (EPSG:4326) - Global GPS
   • NAD 83 (EPSG:4269) - North America
   • ETRS89 (EPSG:4258) - Europe

2. Projection: Transformation from curved to flat surface

   • Cylindrical (Mercator)
   • Conic (Lambert Conformal)
   • Azimuthal (Polar Stereographic)

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

Common CRS Codes

Geographic CRS (Lat/Lon)

Projected CRS (Meters)

Regional Projections

United States:

• EPSG:5070 - US National Atlas Equal Area (CONUS)
• EPSG:6350 - US National Atlas (Alaska)
• EPSG: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 / MGA zone 52 (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
gpd.geographic_crs = "EPSG:4326"
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

# Project to appropriate UTM zone
gdf = gpd.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):
    """Get UTM zone EPSG code from coordinates."""

    zone = 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

# 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)

UTM Non-standard:

• EPSG:31466-31469 - German Gauss-Krüger zones
• EPSG:2056 - Swiss LV95 (based on UTM principles)

Transformations

Basic Transformation

from pyproj import Transformer

# Create transformer
transformer = Transformer.from_crs(
    "EPSG:4326",  # WGS 84 (lat/lon)
    "EPSG:32633", # UTM Zone 33N (meters)
    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

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

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

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

# If None, set it
if gdf.crs is None:
    gdf.set_crs("EPSG:4326", inplace=True)

2. Verify CRS Before Operations

def ensure_same_crs(gdf1, gdf2):
    """Ensure two GeoDataFrames have same CRS."""
    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")  # US National Atlas Equal Area
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