Map Projections in Python: A Practical Guide

Overview

The Earth is (roughly) a sphere. A map is flat. Putting a curved surface onto a plane always introduces distortion. No projection can preserve area, shape, distance, and direction all at once.

Map projections get overlooked. The wrong projection isn’t just ugly, it misleads. Area, spatial relationships, and your actual message suffer.

This guide, ported from Dr. Dominic Royé’s R version, shows you how to pick the right projection for your purpose in Python using geopandas, cartopy, and matplotlib.

Why Projections Matter

Every projection trades one property for another:

Property	Preserved by	Use case
Area	Equal-area (equivalent)	Choropleth, density, species richness
Shape	Conformal	Navigation, large-scale topographic maps
Distance	Equidistant	Maps centered on a specific point
Direction	Azimuthal	Polar maps, great-circle routes

Question: “which distortion works for my purpose?” not “which projection is best?”

Key rule: Match the projection to the purpose. Navigation projections don’t work for choropleth maps.

Packages

Install the required packages:

pip install geopandas cartopy shapely pyproj matplotlib numpy

Package	Purpose
`geopandas`	Vector data (read, transform, manipulate)
`cartopy`	Projections, Tissot’s indicatrix, plot helpers
`matplotlib`	Plotting, subplots, colors
`shapely`	Geometry operations (buffer, intersection)
`pyproj`	CRS handling (EPSG/ESRI/WKT2 codes)

import geopandas as gpd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from shapely.geometry import Point, Polygon
from pyproj import CRS
import numpy as np

The Offenders: Stop Using These (For Most Things)

Mercator | The World’s Most Misused Projection

Created by Gerardus Mercator in 1569 for nautical navigation. It’s conformal (preserves local angles), which is exactly what sailors needed.

The problem: Mercator exaggerates areas away from the equator. At 60° latitude, areas appear about 4 times larger than reality. Greenland looks comparable to Africa, when Africa is actually ~14 times larger.

When Mercator is acceptable:

Web map tiles (OpenStreetMap), because conformal properties make tile seams invisible
Large-scale navigation charts
Conformal mapping of small areas

Never use Mercator for:

Choropleth maps where area conveys information
Global or hemispheric maps
Any map where area-dependent interpretation matters

# Download Natural Earth data
import zipfile, urllib.request, os

ne_url = "https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip"
zip_path = "/tmp/ne_110m_admin_0_countries.zip"

if not os.path.exists("/tmp/ne_110m_admin_0_countries.shp"):
    urllib.request.urlretrieve(ne_url, zip_path)
    with zipfile.ZipFile(zip_path, 'r') as zip_ref:
        zip_ref.extractall("/tmp/")

# Load world data
world = gpd.read_file("/tmp/ne_110m_admin_0_countries.shp")
world = world[world['ADMIN'] != "Antarctica"].copy()
world['area_km2'] = world.geometry.to_crs('ESRI:54035').area / 1e6

# Compare Mercator vs Equal Earth
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8),
                                 subplot_kw={'projection': ccrs.EqualEarth()})

# Equal Earth (correct)
ax1.set_global()
world.plot(ax=ax1, column='area_km2', cmap='YlGn', legend=True,
           edgecolor='white', linewidth=0.15)
ax1.set_title('Equal Earth (ESRI:54035) ✓', fontsize=11, fontweight='bold')
ax1.coastlines(linewidth=0.15)

# Mercator (shows distortion)
ax2 = plt.subplot(1, 2, 2, projection=ccrs.Mercator())
ax2.set_global()
world.plot(ax=ax2, column='area_km2', cmap='YlGn', legend=True,
           edgecolor='white', linewidth=0.15, transform=ccrs.PlateCarree())
ax2.set_title('Mercator (EPSG:3395) ✗', fontsize=11, fontweight='bold')
ax2.coastlines(linewidth=0.15)

plt.suptitle('Area Distortion Comparison', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('mercator_vs_equal_earth.png', dpi=150, bbox_inches='tight')

Geographic Coordinates

Almost as bad as Mercator: using geographic (longitude/latitude) coordinates as a projection. This is the equirectangular projection (EPSG:4326 used as a projection). It’s neither equal-area nor conformal.

In Python, if you plot geopandas data without setting a projection, you’re implicitly using this. Always set an appropriate CRS for your map.

A Taxonomy of Projections

Family	Property Preserved	Typical Use
Equal-area (equivalent)	Area	Choropleth, density, species richness
Conformal	Local shape/angles	Navigation, large-scale maps
Compromise	Neither fully, balanced distortion	General-purpose reference maps
Azimuthal	Direction from center	Polar maps, great-circle routes
Equidistant	Distance from point(s)	Maps centered on a location

Global thematic maps: use equal-area projections

Any global map where region size matters needs an equal-area projection.

Equal Earth (ESRI:54035): Modern replacement for Robinson, designed in 2018. Pleasing appearance, correct relative areas. Recommended for world thematic maps.

Mollweide (ESRI:54009): Classic pseudocylindrical equal-area projection. Used in scientific journals.

Robinson (ESRI:54030): Compromise projection, not equal-area, but modest area distortion. Works for reference maps, not thematic maps.

Natural Earth II (ESRI:54078): Compromise projection with rounded polar appearance. Alternative to Winkel Tripel.

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 12),
                                        subplot_kw={'projection': ccrs.EqualEarth()})

# Equal Earth
ax1.set_global()
world.plot(ax=ax1, color='#3d6b8c', edgecolor='white', linewidth=0.1)
ax1.set_title('Equal Earth (ESRI:54035) ✓', fontsize=10, fontweight='bold')
ax1.coastlines(linewidth=0.1)

# Mollweide
ax2 = plt.subplot(2, 2, 2, projection=ccrs.Mollweide())
ax2.set_global()
world.plot(ax=ax2, color='#3d6b8c', edgecolor='white', linewidth=0.1, transform=ccrs.PlateCarree())
ax2.set_title('Mollweide (ESRI:54009) ✓', fontsize=10, fontweight='bold')

# Robinson
ax3 = plt.subplot(2, 2, 3, projection=ccrs.Robinson())
ax3.set_global()
world.plot(ax=ax3, color='#3d6b8c', edgecolor='white', linewidth=0.1, transform=ccrs.PlateCarree())
ax3.set_title('Robinson (ESRI:54030)', fontsize=10, fontweight='bold')

# Natural Earth II
ax4 = plt.subplot(2, 2, 4, projection=ccrs.NaturalEarth2())
ax4.set_global()
world.plot(ax=ax4, color='#3d6b8c', edgecolor='white', linewidth=0.1, transform=ccrs.PlateCarree())
ax4.set_title('Natural Earth II (ESRI:54078)', fontsize=10, fontweight='bold')

plt.suptitle('Global Projections Comparison', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('global_projections.png', dpi=150, bbox_inches='tight')

Tissot’s indicatrix: visual distortion check

Tissot’s indicatrix places identical circles at regular intervals across the globe. On equal-area projections, all circles have the same area. On conformal projections, all circles stay circular but change size.

from shapely.geometry import Point
from shapely.ops import unary_union

# Generate Tissot circles in WGS84
def generate_tissot(spacing=30, radius_deg=4):
    lons = np.arange(-180, 180, spacing)
    lats = np.arange(-60, 61, spacing)
    circles = []
    for lon in lons:
        for lat in lats:
            circle = Point(lon, lat).buffer(radius_deg, resolution=32)
            circles.append(circle)
    return gpd.GeoDataFrame({'geometry': circles}, crs='EPSG:4326')

tissot = generate_tissot(spacing=30, radius_deg=4)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8),
                                 subplot_kw={'projection': ccrs.Mercator()})

# Mercator - circles stretched at high latitudes
ax1.set_global()
tissot_merc = tissot.to_crs('EPSG:3395')
tissot_merc.plot(ax=ax1, color='#d7362d', alpha=0.75, edgecolor='none')
world.boundary.plot(ax=ax1, linewidth=0.2, color='grey', transform=ccrs.PlateCarree())
ax1.set_title('Mercator - Circles Stretched ✗', fontsize=11, fontweight='bold')

# Equal Earth - circles have equal area
ax2 = plt.subplot(1, 2, 2, projection=ccrs.EqualEarth())
ax2.set_global()
tissot_ee = tissot.to_crs('ESRI:54035')
tissot_ee.plot(ax=ax2, color='#d7362d', alpha=0.75, edgecolor='none')
world.boundary.plot(ax=ax2, linewidth=0.2, color='grey', transform=ccrs.PlateCarree())
ax2.set_title('Equal Earth - Circles Equal Area ✓', fontsize=11, fontweight='bold')

plt.suptitle("Tissot's Indicatrix - Visual Distortion Check", fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('tissot_indicatrix.png', dpi=150, bbox_inches='tight')

Regional and continental maps

For continents or large sub-regions, pick based on the region’s shape (east-west vs. north-south extent).

Europe: Lambert Azimuthal Equal Area (LAEA) centered on Europe, EPSG:3035. The standard for European statistical maps.

North America: Albers Equal Area Conic (EPSG:5070). Standard for USA continental maps.

Spain/Iberian Peninsula: UTM zone 30N (EPSG:25830) for large-scale work.

# Filter Europe
europe = world[world['CONTINENT'] == 'Europe'].copy()
europe = europe[~europe['ADMIN'].isin(['Greenland', 'Svalbard', 'Russia'])]
from shapely.geometry import box
europe = europe.clip(box(-33, 26, 50, 75))

# Transform to LAEA Europe (EPSG:3035)
europe_laea = europe.to_crs("EPSG:3035")

# Plot
fig, ax = plt.subplots(1, 1, figsize=(10, 12),
                       subplot_kw={'projection': ccrs.LambertAzimuthalEqualArea(central_longitude=20, central_latitude=52)})
europe_laea.plot(ax=ax, color='#4a7fb5', edgecolor='white', linewidth=0.2)
ax.set_title('Europe - LAEA (EPSG:3035) ✓', fontsize=11, fontweight='bold')
ax.coastlines(linewidth=0.15)
plt.tight_layout()
plt.savefig('europe_laea.png', dpi=150, bbox_inches='tight')

Orthographic projection: globe view

The orthographic projection (+proj=ortho) renders Earth as viewed from space. Compression toward edges is accurate foreshortening, not a flaw.

# Simplified approach using cartopy
fig, ax = plt.subplots(1, 1, figsize=(10, 10),
                       subplot_kw={'projection': ccrs.Orthographic(central_longitude=0.5, central_latitude=51)})
ax.add_feature(cfeature.OCEAN, facecolor='#cce4f0')
ax.add_feature(cfeature.LAND, facecolor='#3d6b8c', edgecolor='white', linewidth=0.1)
ax.gridlines(draw_labels=True, linewidth=0.1, color='gray', alpha=0.5)
ax.set_title('Orthographic - Europe Centered (lat_0=51, lon_0=0.5)', fontsize=11, fontweight='bold')
plt.tight_layout()
plt.savefig('orthographic.png', dpi=150, bbox_inches='tight')

Specifying Projections in Python

Vector Data with `geopandas`

# Inspect CRS
print(gpd.get_crs(world))  # or world.crs
print(world.crs.to_epsg())  # EPSG code
print(world.crs.name)         # Human-readable name

# Transform to different CRS
world_ee = world.to_crs("ESRI:54035")   # Equal Earth
world_merc = world.to_crs("EPSG:3395")  # Mercator
world_laea = world.to_crs("EPSG:3035")  # LAEA Europe

# Get WKT2 string (best for archiving/sharing)
wkt2 = world.crs.to_wkt()
print(wkt2[:200])

CRS Code Formats

Format	Example	When to Use
EPSG	`"EPSG:4326"`	Standard geodetic CRS (preferred)
ESRI	`"ESRI:54035"`	World projections not in EPSG
WKT2	`CRS.from_wkt(...)`	Archiving, sharing, portable
PROJ4	`'+proj=ortho +lat_0=51'`	Legacy (avoid in new code)

from pyproj import CRS

# EPSG code
crs_epsg = CRS.from_epsg(4326)

# ESRI code
crs_esri = CRS.from_string("ESRI:54035")

# WKT2 string
wkt2 = crs_epsg.to_wkt()
crs_from_wkt = CRS.from_wkt(wkt2)

Quick-Reference: Projection Codes

Use Case	Projection	Code
Global thematic maps	Equal Earth	`ESRI:54035`
Global thematic maps (alternative)	Mollweide	`ESRI:54009`
World reference maps	Robinson	`ESRI:54030`
World reference maps (alternative)	Natural Earth II	`ESRI:54078`
European thematic maps	LAEA Europe	`EPSG:3035`
USA thematic maps	Albers Equal Area	`EPSG:5070`
Iberian Peninsula	UTM zone 30N	`EPSG:25830`
Arctic maps	NSIDC Polar Stereographic	`EPSG:3413`
Web/tiled maps	Web Mercator	`EPSG:3857`
Navigation charts	Mercator	`EPSG:3395`
Data storage/exchange	WGS84 geographic	`EPSG:4326`

# Visualize all common projections
projections = {
    'Equal Earth': ccrs.EqualEarth(),
    'Mollweide': ccrs.Mollweide(),
    'Robinson': ccrs.Robinson(),
    'Mercator': ccrs.Mercator(),
    'LAEA Europe': ccrs.LambertAzimuthalEqualArea(central_longitude=20, central_latitude=52),
}

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()

for idx, (name, crs) in enumerate(projections.items()):
    axes[idx].set_projection(crs)
    axes[idx].set_global()
    world.plot(ax=axes[idx], color='#3d6b8c', edgecolor='white',
               linewidth=0.1, transform=ccrs.PlateCarree())
    axes[idx].set_title(name, fontsize=9)
    axes[idx].coastlines(linewidth=0.1)

plt.suptitle('Quick-Reference: Common Projections', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('quick_reference.png', dpi=150, bbox_inches='tight')

Checklist: Common Projection Mistakes

Before submitting a map, check:

Using Mercator for thematic maps? Replace with equal-area projection
Using geographic coordinates (EPSG:4326) as projection? Reproject to suitable CRS
Wrong projection for region? Match projection to spatial extent
Choropleth without equal-area? Area-dependent maps need equal-area projections
Using st_area() on WGS84 without s2? Project to equal-area CRS first

References

Original R guide: Dr. Dominic Royé - Map Projections
GeoPandas docs: geopandas.org
Cartopy docs: scitools.org.uk/cartopy
PyProj docs: pyproj4.github.io
Natural Earth: naturalearthdata.com