Cross-archive Mars: combining NASA PDS and ESA PSA over one area

Finding and working with multiple datasets over Jezero Crater

NoteWhat this tutorial shows

How to assemble data over one Mars surface area from two archives — NASA’s PDS and ESA’s PSA — and tie them together with the rest of planetarypy: per-body constants, planetary crs, geo footprints/overlaps, spicer illumination, and the (experimental) sun-direction plotting helper.

It also documents, honestly, the rough edges you hit doing this today — they are called out in callout-warning boxes so you know what works, what doesn’t yet, and the reliable path around each.

We work over Jezero Crater (Perseverance’s landing site), at about 77.4°E, 18.4°N.

LON, LAT = 77.4, 18.4   # Jezero Crater (degrees East, degrees North)
TARGET = "Mars"

1. Ground truth: the body and a local CRS

Before any data, anchor the work in the body’s own numbers and a sensible coordinate system for the area.

from planetarypy import constants as c

mars = c.Mars
mars.mean_radius   # a Constant with units + provenance

\(3389.5267 \; \mathrm{km}\)

planetarypy.crs builds IAU-2015 CRSs straight from the body — radii come from the IAU code itself, nothing hardcoded. For a feature study you usually want a local azimuthal-equidistant CRS centred on the area (distance-true near the point), plus the body-wide geographic CRS.

from planetarypy import crs

body_geographic = crs.body_crs(TARGET)          # ('IAU_2015', '49900')
local = crs.local_crs(LON, LAT, TARGET)         # azeqd centred on Jezero
body_geographic.to_authority()

('IAU_2015', '49900')

2. ESA PSA: spatially searching one area

The PSA exposes an IVOA EPN-TAP table (psa.epn_core) where every granule is a directly downloadable product. The high-level helpers (psa.missions, psa.instruments, psa.datasets) are non-spatial — but the table itself carries spatial bounds, and psa.query() lets you run raw ADQL against them.

from planetarypy import psa

psa.instruments("Mars Express").head()

	mission	instrument	products
0	Mars Express	ASPERA-3	675880
1	Mars Express	HRSC	265494
2	Mars Express	VMC	252989
3	Mars Express	MaRS	231249
4	Mars Express	SPICAM	89427

The spatial columns are c1min/c1max (longitude) and c2min/c2max (latitude).

WarningObstacle 1 — longitude is −180..180, and bounding boxes degenerate near the poles

Two things will bite you:

1. Longitude is signed (−180..180), not 0–360 — undocumented in the API. Jezero at “77.4°E” is 77.4 here, but a feature near the antimeridian needs care.
2. Orbital swaths that cross high latitudes get a degenerate bounding box that spans all longitudes (lon[-180, 180]). A naive “box contains my point” query then returns lots of false positives that don’t actually image your area.

The fix: drop the wrap-around boxes by also requiring the box to be small ((c1max - c1min) < 30). For rigorous work, use the per-product geometry index (psa.geometry_index, which carries CENTER_LAT/LON) or the s_region polygon.

hrsc = psa.query(
    "SELECT TOP 10 granule_uid, c1min, c1max, c2min, c2max "
    "FROM psa.epn_core "
    "WHERE instrument_host_name LIKE '%Mars Express%' "
    "  AND instrument_name LIKE '%HRSC%' "
    f"  AND c1min <= {LON} AND c1max >= {LON} "
    f"  AND c2min <= {LAT} AND c2max >= {LAT} "
    "  AND (c1max - c1min) < 30 AND (c2max - c2min) < 30"   # <-- kill degenerate boxes
)
[r["granule_uid"].split(":DATA:")[-1] for r in hrsc]

['HM069_0000_BL3.IMG::1.0',
 'HM069_0000_GR3.IMG::1.0',
 'HM069_0000_IR3.IMG::1.0',
 'HM069_0000_RE3.IMG::1.0',
 'HM612_0000_IR3.IMG::3.0',
 'HQ211_0000_BL3.IMG::5.0',
 'HQ211_0000_GR3.IMG::5.0',
 'HQ211_0000_IR3.IMG::5.0',
 'HQ211_0000_ND3.IMG::5.0',
 'HQ211_0000_P13.IMG::5.0']

Those are real HRSC scenes over Jezero. A granule_uid’s product-id segment is what you hand to psa.fetch_psa_product.

3. NASA PDS: registry-wide discovery with a first-class bbox

planetarypy.search wraps NASA’s registry-wide search API (80M+ products). Spatial search is a first-class bbox= parameter — no raw query string needed. bbox_from_point builds a (west, south, east, north) box around a coordinate, and count() sizes a result set without fetching any rows.

from planetarypy import search

bbox = search.bbox_from_point(LON, LAT, 1.0)     # ±1° box around Jezero
search.count(
    target="urn:nasa:pds:context:target:planet.mars",
    observationals=True,
    bbox=bbox,
)

968

That is every cart-indexed product whose footprint overlaps the box. Pull a page of them — bbox= replaces the old cart:Bounding_Coordinates query-string incantation:

mars_box = search.search_products(
    target="urn:nasa:pds:context:target:planet.mars",
    observationals=True,
    bbox=bbox,
    fields=["ref_lid_instrument", "pds:Time_Coordinates.pds:start_date_time"],
    limit=8,
)
mars_box["ref_lid_instrument"].dropna().head()

lidvid
urn:esa:psa:em16_tgo_cas:data_derived:cas_der_my35_010783_165_0_nir_558955156::1.0    urn:esa:psa:context:instrument:tgo.cassis
urn:esa:psa:em16_tgo_cas:data_derived:cas_der_my35_009577_165_0_red_556485268::1.0    urn:esa:psa:context:instrument:tgo.cassis
urn:esa:psa:em16_tgo_cas:data_derived:cas_der_my35_008085_163_0_blu_553429644::1.0    urn:esa:psa:context:instrument:tgo.cassis
urn:esa:psa:em16_tgo_cas:data_derived:cas_der_my35_007607_017_0_blu_552450116::1.0    urn:esa:psa:context:instrument:tgo.cassis
urn:esa:psa:em16_tgo_cas:data_derived:cas_der_my35_010870_164_0_pan_559133328::1.0    urn:esa:psa:context:instrument:tgo.cassis
Name: ref_lid_instrument, dtype: str

From the command line the same “what is here?” question is a one-liner:

plp search at Mars 77.4 18.4 -r 1.0          # list products overlapping the point
plp search at Mars 77.4 18.4 --count         # just the number

WarningObstacles on the NASA side (two real ones)

1. Spatial fields are not universal. The cart:Bounding_Coordinates fields are only populated when the instrument/archive team added them — common for derived/calibrated products, often absent for raw/EDR. No coverage there means no spatial hit, even if data exists.
2. target=Mars returns every Mars-targeted product, across all missions and nodes — not just orbital imagery. For example, Hayabusa2 NIRS3 cruise-phase spectra that observed Mars (remotely, during transit to Ryugu) are genuinely Mars-targeted and sort first. They’re not mislabeled; they’re just not what you want if you’re after CTX/HiRISE/CaSSIS imagery. Narrow with an instrument/lid constraint.

(Historical note: a missing-outer-parenthesis bug in search’s query builder used to make any multi-filter query fail with HTTP 400; fixed — the registry requires the whole q wrapped in (...).)

Alternative NASA spatial path: the index system

Where cart: coverage is missing (e.g. CTX EDR), planetarypy’s parsed PDS cumulative indexes carry real footprint corners. get_index reads them; filter on the corner lat/lon columns to select an area.

from planetarypy import pds

# CTX EDR index (large first download, then cached as parquet).
ctx = pds.get_index("mro.ctx.edr")
near = ctx[
    (ctx["CENTER_LONGITUDE"].between(LON - 1, LON + 1))
    & (ctx["CENTER_LATITUDE"].between(LAT - 1, LAT + 1))
]
near[["PRODUCT_ID", "CENTER_LONGITUDE", "CENTER_LATITUDE"]].head()

4. Ready-made reference rasters: windowed reads, no full download

Discovery finds products to fetch. Often you also want a ready-made global reference raster over the area — a basemap or a DTM — without pulling the whole multi-gigabyte mosaic. planetarypy.datasets is a body-namespaced registry of remote cloud-optimised GeoTIFFs that reads just the window you ask for over /vsicurl/.

from planetarypy import datasets

datasets.bodies()                                   # ['mars', 'moon']
[r.short for r in datasets.list_datasets("mars")]   # what's registered for Mars

['mars', 'moon']

['hrsc_level3', 'themis_mosaics', 'ctx_dtms']

# A 0.5°-wide patch of the HRSC level-3 global mosaic, centred on Jezero.
# Only the bytes covering this window cross the network — the source COG is global.
patch = datasets.mars.hrsc_level3.read_window(LON, LAT, size=0.5)
patch                                                # georeferenced rioxarray DataArray

<xarray.DataArray 'hrsc_level3' (band: 1, y: 593, x: 593)> Size: 703kB
array([[[7235, 7300, 7577, ..., 8490, 8441, 8333],
        [7341, 7680, 8062, ..., 8134, 8352, 8277],
        [7382, 7753, 8102, ..., 7954, 8172, 8164],
        ...,
        [7610, 6515, 6225, ..., 6653, 6547, 6803],
        [7467, 6573, 6248, ..., 6824, 6777, 6980],
        [7169, 6829, 6139, ..., 7213, 7095, 7124]]],
      shape=(1, 593, 593), dtype=int16)
Coordinates:
  * band         (band) int64 8B 1
    spatial_ref  int64 8B 0
Dimensions without coordinates: y, x
Attributes:
    _FillValue:  -32768

xarray.DataArray

'hrsc_level3'

• band: 1
• y: 593
• x: 593

• 7235 7300 7577 7743 7472 7189 7326 ... 6904 7244 7316 7213 7095 7124

  array([[[7235, 7300, 7577, ..., 8490, 8441, 8333],
          [7341, 7680, 8062, ..., 8134, 8352, 8277],
          [7382, 7753, 8102, ..., 7954, 8172, 8164],
          ...,
          [7610, 6515, 6225, ..., 6653, 6547, 6803],
          [7467, 6573, 6248, ..., 6824, 6777, 6980],
          [7169, 6829, 6139, ..., 7213, 7095, 7124]]],
        shape=(1, 593, 593), dtype=int16)

• Coordinates: (2)
  • band
    (band)
    int64
    1

    array([1])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["Mars (2015) - Sphere / Ocentric / Equirectangular, clon = 0",GEOGCS["Mars (2015) - Sphere / Ocentric",DATUM["Mars (2015) - Sphere",SPHEROID["Mars (2015) - Sphere",3396190,0]],PRIMEM["Reference Meridian",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Equirectangular"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    spatial_ref :

    PROJCS["Mars (2015) - Sphere / Ocentric / Equirectangular, clon = 0",GEOGCS["Mars (2015) - Sphere / Ocentric",DATUM["Mars (2015) - Sphere",SPHEROID["Mars (2015) - Sphere",3396190,0]],PRIMEM["Reference Meridian",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Equirectangular"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    GeoTransform :

    4573000.0 50.0 0.0 1105500.0 0.0 -50.0

    array(0)

• Attributes: (1)

  _FillValue :

  -32768

The result is a rioxarray DataArray in the mosaic’s own IAU CRS (49910, equirectangular), so it composes with crs/geo and plots directly. The STAC-backed collections (THEMIS mosaics, CTX/LOLA DTMs) read the same way; pass out="patch.tif" to also write a GeoTIFF alongside the array.

5. geo: footprints, anti-meridian, and overlaps

Turn the PSA bounding boxes into polygons and reason about them with geo. This is also where the planetarypy anti-meridian helpers earn their keep — exactly the degeneracy from Obstacle 1.

from shapely.geometry import box
import geopandas as gpd

polys = [box(r["c1min"], r["c2min"], r["c1max"], r["c2max"]) for r in hrsc]
# NOTE: geo.overlaps expects the id column to be named exactly "id"
# (it errors with a cryptic KeyError: 'id_left' otherwise).
gdf = gpd.GeoDataFrame(
    {"id": [r["granule_uid"].split(":DATA:")[-1] for r in hrsc]},
    geometry=polys, crs="EPSG:4326",
)
gdf.head()

	id	geometry
0	HM069_0000_BL3.IMG::1.0	POLYGON ((78.38914 15.28595, 78.38914 20.95471...
1	HM069_0000_GR3.IMG::1.0	POLYGON ((78.36745 15.29424, 78.36745 20.94038...
2	HM069_0000_IR3.IMG::1.0	POLYGON ((78.33126 15.23663, 78.33126 20.9729,...
3	HM069_0000_RE3.IMG::1.0	POLYGON ((78.43064 15.2684, 78.43064 21.01134,...
4	HM612_0000_IR3.IMG::3.0	POLYGON ((85.59466 -4.92346, 85.59466 25.0253,...

from planetarypy import geo

# Which scenes overlap each other over the area? -> id_1, id_2, geometry, area
geo.overlaps(gdf).head()

	id_1	id_2	geometry	area
0	HM069_0000_BL3.IMG::1.0	HM612_0000_IR3.IMG::3.0	POLYGON ((76.473 15.28595, 76.473 20.95471, 78...	10.862121
1	HM069_0000_BL3.IMG::1.0	HQ211_0000_IR3.IMG::5.0	POLYGON ((76.473 15.28595, 76.473 20.95471, 78...	10.862121
2	HM069_0000_BL3.IMG::1.0	HQ211_0000_P13.IMG::5.0	POLYGON ((76.473 15.28595, 76.473 20.95471, 78...	10.862121
3	HM069_0000_BL3.IMG::1.0	HQ211_0000_GR3.IMG::5.0	POLYGON ((76.473 15.28595, 76.473 20.95471, 78...	10.862121
4	HM069_0000_BL3.IMG::1.0	HQ211_0000_BL3.IMG::5.0	POLYGON ((76.473 15.28595, 76.473 20.95471, 78...	10.862121

6. spicer: illumination at the area

How was Jezero lit at a given time? Spicer computes full illumination geometry (needs the [spice] extra).

from planetarypy.spice.spicer import Spicer

s = Spicer("Mars")
illum = s.illumination(LON, LAT, "2021-02-18T20:55:00")   # Perseverance landing
sun_az_from_north = s.solar_azimuth_at(LON, LAT, "2021-02-18T20:55:00")
illum, sun_az_from_north

WarningObstacle 3 — solar azimuth conventions don’t match between modules

Spicer.solar_azimuth_at returns azimuth clockwise from north (geographic), but the plotting helper add_sun_indicator (below) expects clockwise from image top (the PDS SUB_SOLAR_AZIMUTH convention). They agree only when image-north points up. For a rotated/projected image you must convert (subtract the image’s north azimuth) before plotting. This mismatch is why the sun helper is experimental.

7. Fetch and open a product

Both archives download into {storage_root} and open through the same planetarypy.open().

import planetarypy as plp

# ESA PSA — direct FTP fetch by product id:
paths = plp.psa.fetch_psa_product("H5270_0000_ND3")
ds = plp.open(paths[0])     # rioxarray/pdr-backed reader
ds

8. Plotting with the (experimental) sun indicator

import planetarypy as plp

# image = ds.squeeze().values   # a 2D array from the opened product
# CW-from-top azimuth required (see Obstacle 3 — convert from Spicer if needed).
plp.plotting.imshow_with_sun(image, sun_azimuth_deg=135, title="HRSC over Jezero")

WarningObstacle 4 — add_sun_indicator / imshow_with_sun are experimental

The placement was recently fixed (it no longer rescales the image), but the azimuth-convention handoff (Obstacle 3) is not yet validated end-to-end. Treat the sun arrow as indicative until that’s resolved.

Summary

• Both archives support spatial discovery. NASA is now a first-class bbox= parameter on search.search_products / search.count (plus the plp search at CLI); ESA is still an EPN-TAP bbox in a raw psa.query ADQL string. Both carry data-quality traps (NASA: missing cart coverage + cross-mission products that legitimately match target=Mars; PSA: degenerate polar bounding boxes).
• Where registry spatial metadata is missing (e.g. CTX EDR), planetarypy’s PDS index system (pds.get_index + corner lat/lon) is the fallback.
• For a ready-made reference raster over the area, planetarypy.datasets streams a windowed read over /vsicurl/ (datasets.mars.hrsc_level3.read_window(...)) — no full-mosaic download.
• constants, crs, geo, spicer, open, and plotting compose cleanly across both archives — the cross-archive friction is in discovery, not in the analysis layer.

See the obstacles called out above; several are tracked for fixes.