Cartopy Plotting

Cartopy / Matplotlib Plotting Skills

Use Cartopy/Matplotlib for Python-based geographic maps when the user asks for a general map, scatter map, station map, or cross-section and does not specify GMT.

Do not use this skill when the user explicitly asks for GMT, gmt, Generic Mapping Tools, gmt begin, gmt coast, or run_gmt; use gmt_plotting instead.

1. Basic earthquake location map (scatter)

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np

# --- lon, lat, depth (and optionally magnitude) are numpy arrays ---
pad = max((lon.max()-lon.min())*0.1, (lat.max()-lat.min())*0.1, 0.5)
extent = [lon.min()-pad, lon.max()+pad, lat.min()-pad, lat.max()+pad]

fig, ax = plt.subplots(figsize=(10, 8),
                        subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_extent(extent, crs=ccrs.PlateCarree())

# Background features
ax.add_feature(cfeature.LAND,      facecolor='#f0ede5', zorder=0)
ax.add_feature(cfeature.OCEAN,     facecolor='#d6eaf8', zorder=0)
ax.add_feature(cfeature.COASTLINE, linewidth=0.8, color='#444', zorder=2)
ax.add_feature(cfeature.BORDERS,   linewidth=0.5, linestyle=':', color='#777', zorder=2)
ax.add_feature(cfeature.RIVERS,    linewidth=0.3, color='#90caf9', zorder=1)

# Grid labels
gl = ax.gridlines(draw_labels=True, linewidth=0.4, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False

# Plot earthquakes coloured by depth
sc = ax.scatter(lon, lat, c=depth, cmap='plasma_r', s=20,
                transform=ccrs.PlateCarree(), zorder=5,
                alpha=0.85, edgecolors='none',
                vmin=0, vmax=depth.max() if depth.max() > 0 else 1)
plt.colorbar(sc, ax=ax, label='Depth (km)', shrink=0.65, pad=0.02)
ax.set_title('Seismicity Map', fontsize=13, pad=10)

savefig('seismicity_map.png')
plt.close()

2. Magnitude-scaled scatter (size ∝ magnitude)

# mag is a numpy array
sizes = (2 ** mag) * 3          # exponential scaling — looks natural
sc = ax.scatter(lon, lat, c=depth, cmap='plasma_r',
                s=sizes, transform=ccrs.PlateCarree(),
                alpha=0.7, edgecolors='black', linewidths=0.3, zorder=5)

3. Station + event map (two symbol types)

# Events
ax.scatter(ev_lon, ev_lat, c='red', s=20, marker='o',
           transform=ccrs.PlateCarree(), zorder=5, label='Events')
# Stations
ax.scatter(st_lon, st_lat, c='blue', s=60, marker='^',
           transform=ccrs.PlateCarree(), zorder=6, label='Stations')
ax.legend(loc='lower right', fontsize=9)

4. Coloured by time (seismicity evolution)

import matplotlib.dates as mdates
import pandas as pd

# times is a list of datetime or timestamp strings → convert to float for colourmap
t_num = mdates.date2num(pd.to_datetime(times))
sc = ax.scatter(lon, lat, c=t_num, cmap='viridis',
                s=15, transform=ccrs.PlateCarree(), zorder=5, alpha=0.8)
cb = plt.colorbar(sc, ax=ax, shrink=0.65, pad=0.02, label='Time')
cb.ax.yaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))

5. SRTM-like hill shading with cartopy (no internet required)

# Uses NaturalEarth built-in shaded relief at ~10m resolution
ax.stock_img()                  # NASA Blue Marble background (requires cartopy data)
# OR use a simple terrain colour from OSM tiles (requires internet):
# import cartopy.io.img_tiles as cimgt
# ax.add_image(cimgt.Stamen('terrain-background'), 8)  # zoom level 8

6. Add fault lines from a file (lon/lat pairs separated by >)

import numpy as np

# fault_file has lon lat rows; segments separated by blank lines or '>'
faults = np.loadtxt(fault_file, comments='>')
# OR read segment by segment:
with open(fault_file) as f:
    seg, segs = [], []
    for line in f:
        line = line.strip()
        if not line or line.startswith('>'):
            if seg:
                segs.append(np.array(seg))
                seg = []
        else:
            seg.append([float(x) for x in line.split()])
    if seg:
        segs.append(np.array(seg))

for s in segs:
    ax.plot(s[:, 0], s[:, 1], 'k-', linewidth=0.5,
            transform=ccrs.PlateCarree(), zorder=3)

7. Cross-section: depth vs epicentral distance

from obspy.geodetics import gps2dist_azimuth

# reference point
ref_lon, ref_lat = 110.0, 35.0
dist_km = np.array([gps2dist_azimuth(ref_lat, ref_lon, la, lo)[0]/1000
                    for lo, la in zip(lon, lat)])

fig, ax = plt.subplots(figsize=(10, 5))
sc = ax.scatter(dist_km, depth, c=mag, cmap='hot_r', s=20, alpha=0.7,
                edgecolors='none', vmin=mag.min(), vmax=mag.max())
plt.colorbar(sc, ax=ax, label='Magnitude', shrink=0.8)
ax.set_xlabel('Distance (km)')
ax.set_ylabel('Depth (km)')
ax.invert_yaxis()
ax.set_title('Vertical Cross-Section')
ax.grid(True, linewidth=0.3, alpha=0.5)
savefig('cross_section.png')
plt.close()

8. Multi-panel: map + cross-section side by side

fig = plt.figure(figsize=(14, 6))
ax_map = fig.add_subplot(1, 2, 1, projection=ccrs.PlateCarree())
ax_xsec = fig.add_subplot(1, 2, 2)

# --- Map panel ---
ax_map.set_extent(extent, crs=ccrs.PlateCarree())
ax_map.add_feature(cfeature.LAND,      facecolor='#f0ede5')
ax_map.add_feature(cfeature.OCEAN,     facecolor='#d6eaf8')
ax_map.add_feature(cfeature.COASTLINE, linewidth=0.7)
sc1 = ax_map.scatter(lon, lat, c=depth, cmap='plasma_r', s=15,
                      transform=ccrs.PlateCarree(), alpha=0.8)
fig.colorbar(sc1, ax=ax_map, label='Depth (km)', shrink=0.6)
ax_map.set_title('Epicentre Map', fontsize=11)

# --- Cross-section panel ---
ax_xsec.scatter(dist_km, depth, c=mag, cmap='hot_r', s=15, alpha=0.7)
ax_xsec.invert_yaxis()
ax_xsec.set_xlabel('Distance (km)')
ax_xsec.set_ylabel('Depth (km)')
ax_xsec.set_title('Cross-Section', fontsize=11)
ax_xsec.grid(True, linewidth=0.3, alpha=0.4)

plt.tight_layout()
savefig('map_and_xsec.png')
plt.close()

9. Projection variants

# Lambert Conformal — good for mid-latitude regions
ccrs.LambertConformal(central_longitude=lon.mean(), central_latitude=lat.mean())

# Albers Equal Area — good for large countries
ccrs.AlbersEqualArea(central_longitude=lon.mean(), standard_parallels=(25, 47))

# PlateCarree — simple, global or regional
ccrs.PlateCarree()

# Mercator — avoid for high latitudes
ccrs.Mercator()

Common mistakes to avoid

Tip: saving multiple figures

for i, (title, data) in enumerate(panels):
    fig, ax = plt.subplots(figsize=(8, 6))
    ax.plot(data['x'], data['y'])
    ax.set_title(title)
    savefig(f'panel_{i:02d}.png')   # each call registers one [FIGURE]
    plt.close()