Planetary Coordinate Reference Systems

planetarypy.crs builds planetary CRS from IAU codes with pyproj. The body’s ellipsoid/radii come straight from the IAU code — nothing is looked up or hardcoded. It also builds feature-centered local projections for distance-true work near a point (craters, landing sites). This module was adapted from Christian Tai Udovicic’s craterpy.

from planetarypy.crs import body_crs, local_crs, get_crs

A body’s CRS

Pass a body name (resolved via planetarypy.constants) or its NAIF id. The radii come from the code itself:

mars = body_crs("mars")          # or body_crs(499)
print(mars.name)
print("geographic?", mars.is_geographic)
print("semi-major axis (m):", mars.ellipsoid.semi_major_metre)

Mars (2015) - Sphere / Ocentric
geographic? True
semi-major axis (m): 3396190.0

On a fresh machine the name path triggers a one-time NSSDC download (via planetarypy.constants). Pass the NAIF id integer to skip that:

assert body_crs(499).to_wkt() == body_crs("mars").to_wkt()
print("499 == 'mars' ✓")

499 == 'mars' ✓

Ocentric vs ographic

system="ocentric" (the spherical IAU definition) exists for every body; "ographic" only for some. Asking for one a body lacks raises a clear error:

print(body_crs(499, system="ographic").name)
try:
    body_crs(301, system="ographic")   # the Moon has no ographic code
except ValueError as e:
    print("Moon ographic ->", e)

Mars (2015) / Ographic
Moon ographic -> No IAU_2015 'ographic' CRS for body 301 (code 30101). The body may only define an 'ocentric' system, or be absent from the IAU_2015 edition.

Local, feature-centered projection

local_crs(lon, lat, body) returns an Azimuthal-Equidistant CRS centered on a point, built on the body’s IAU geodetic CRS (so the sphere is the body’s). Distances and small buffers near the center are true — ideal for crater annuli or ROI buffers. Here: a 50 km circle around Gale crater.

from pyproj import Transformer
from shapely.geometry import Point as SPoint
from shapely.ops import transform as shp_transform

lon, lat = 137.4, -4.6
lcrs = local_crs(lon, lat, "mars")
print(lcrs.name, "| projected?", lcrs.is_projected)

# A 50 km circle is trivial in the local metric CRS (center is 0,0)...
circle_local = SPoint(0, 0).buffer(50_000)   # 50 km, in metres
# ...transform it back to lon/lat to place it on the globe:
to_lonlat = Transformer.from_crs(
    lcrs, body_crs("mars").geodetic_crs, always_xy=True).transform
circle_lonlat = shp_transform(to_lonlat, circle_local)
print("50 km circle lon/lat bounds:",
      [round(v, 3) for v in circle_lonlat.bounds])

AzimuthalEquidistant(-4.6000, 137.4000) on mars | projected? True
50 km circle lon/lat bounds: [136.554, -5.444, 138.246, -3.756]

Footprints across the anti-meridian

Polygons that cross ±180° longitude break naive shapely geometry. planetarypy.geo.split_at_antimeridian splits them into hemisphere pieces (and handles pole-containing caps):

from planetarypy.geo import split_at_antimeridian

# a footprint straddling the antimeridian (170°E .. -170°E)
corners = [(170, -10), (-170, -10), (-170, 10), (170, 10)]
parts = split_at_antimeridian(corners)
print(f"{len(parts)} polygon(s):")
for p in parts:
    lons = [c[0] for c in p.exterior.coords]
    print(f"  lon {min(lons):.0f}..{max(lons):.0f}")

2 polygon(s):
  lon 170..180
  lon -180..-170

craterpy compatibility

get_crs(body, system="default") mirrors craterpy’s entry point ("default" → ocentric):

get_crs("mars", "default").name

'Mars (2015) - Sphere / Ocentric'