dataset.core.base

base

Abstract base for Mars HiRISE dataset variants (RDR and DTM).

ProductMeta dataclass

ProductMeta(scaling_factor: float = _DEFAULT_SCALING_FACTOR, offset: float = _DEFAULT_OFFSET, sample_bits: int = _DEFAULT_SAMPLE_BITS, effective_max_dn: int = (1 << _EFFECTIVE_BIT_DEPTH) - 1, filter_names: list[str] = (lambda: ['NEAR-INFRARED', 'RED', 'BLUE-GREEN'])(), bands: int = 3, incidence_angle: float = 45.0, solar_azimuth: float = 270.0, emission_angle: float | None = None, phase_angle: float | None = None, local_time: float | None = None, solar_longitude: float | None = None, sub_solar_azimuth: float | None = None, north_azimuth: float | None = None, observation_id: str | None = None, product_id: str | None = None, product_version_id: str | None = None, target_name: str | None = None, mission_phase_name: str | None = None, orbit_number: int | None = None, rationale_desc: str | None = None, start_time: str | None = None, stop_time: str | None = None, product_creation_time: str | None = None, map_scale: float | None = None, map_resolution: float | None = None, map_projection_type: str | None = None, center_latitude: float | None = None, center_longitude: float | None = None, center_filter_wavelength: float | None = None)

Radiometric metadata parsed from a per-product PDS3 LBL file.

Shared across RDR and DTM ortho products. The defaults correspond to representative COLOR RDR label values; DTM ortho images should parse their own LBL which typically has 8-bit samples and different offsets.

from_lbl classmethod

from_lbl(lbl_path: Path) -> ProductMeta

Parse a PDS3 LBL and return a populated instance.

Source code in src/dataset/core/base.py

@classmethod
def from_lbl(cls, lbl_path: pathlib.Path) -> ProductMeta:
    """Parse a PDS3 LBL and return a populated instance."""
    obj = cls()
    if not lbl_path.exists():
        logger.debug("LBL not found, using defaults: %s", lbl_path)
        return obj
    try:
        text = lbl_path.read_text(errors="replace")
    except OSError as exc:
        logger.warning("Could not read LBL %s: %s", lbl_path, exc)
        return obj

    def _float(pat: str) -> float | None:
        m = re.search(pat, text, re.MULTILINE)
        return float(m.group(1)) if m else None

    def _int(pat: str) -> int | None:
        m = re.search(pat, text, re.MULTILINE)
        return int(m.group(1)) if m else None

    def _str(pat: str) -> str | None:
        m = re.search(pat, text, re.MULTILINE)
        return m.group(1).strip().strip('"').strip("'") if m else None

    if (v := _float(r"^\s*INCIDENCE_ANGLE\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.incidence_angle = v
    # Match SOLAR_AZIMUTH but not SUB_SOLAR_AZIMUTH
    if (v := _float(r"^\s*(?<!SUB_)SOLAR_AZIMUTH\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.solar_azimuth = v
    if (v := _float(r"^\s*EMISSION_ANGLE\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.emission_angle = v
    if (v := _float(r"^\s*PHASE_ANGLE\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.phase_angle = v
    if (v := _float(r"^\s*LOCAL_TIME\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.local_time = v
    if (v := _float(r"^\s*SOLAR_LONGITUDE\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.solar_longitude = v
    if (v := _float(r"^\s*SUB_SOLAR_AZIMUTH\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.sub_solar_azimuth = v
        # If solar_azimuth wasn't set from SOLAR_AZIMUTH, fall back here.
        if not re.search(r"^\s*(?<!SUB_)SOLAR_AZIMUTH\s*=", text, re.MULTILINE):
            obj.solar_azimuth = v
    if (v := _float(r"^\s*NORTH_AZIMUTH\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.north_azimuth = v

    # Identification / timing
    obj.observation_id = _str(r"^\s*OBSERVATION_ID\s*=\s*\"?([^\"\n]+?)\"?\s*$")
    obj.product_id = _str(r"^\s*PRODUCT_ID\s*=\s*\"?([^\"\n]+?)\"?\s*$")
    obj.product_version_id = _str(r"^\s*PRODUCT_VERSION_ID\s*=\s*\"?([^\"\n]+?)\"?\s*$")
    obj.target_name = _str(r"^\s*TARGET_NAME\s*=\s*\"?([^\"\n]+?)\"?\s*$")
    obj.mission_phase_name = _str(r"^\s*MISSION_PHASE_NAME\s*=\s*\"?([^\"\n]+?)\"?\s*$")
    obj.rationale_desc = _str(r"^\s*RATIONALE_DESC\s*=\s*\"?([^\"\n]+?)\"?\s*$")
    obj.start_time = _str(r"^\s*START_TIME\s*=\s*([^\s\n]+)")
    obj.stop_time = _str(r"^\s*STOP_TIME\s*=\s*([^\s\n]+)")
    obj.product_creation_time = _str(
        r"^\s*PRODUCT_CREATION_TIME\s*=\s*([^\s\n]+)"
    )
    if (v := _int(r"^\s*ORBIT_NUMBER\s*=\s*(\d+)")) is not None:
        obj.orbit_number = v

    # Map projection
    if (v := _float(r"^\s*MAP_SCALE\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.map_scale = v
    if (v := _float(r"^\s*MAP_RESOLUTION\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.map_resolution = v
    obj.map_projection_type = _str(
        r"^\s*MAP_PROJECTION_TYPE\s*=\s*\"?([^\"\n]+?)\"?\s*$"
    )
    if (v := _float(r"^\s*CENTER_LATITUDE\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.center_latitude = v
    if (v := _float(r"^\s*CENTER_LONGITUDE\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.center_longitude = v
    if (v := _float(r"^\s*CENTER_FILTER_WAVELENGTH\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.center_filter_wavelength = v

    if (v := _float(r"^\s*SCALING_FACTOR\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.scaling_factor = v
    if (v := _float(r"^\s*OFFSET\s*=\s*([\d.eE+\-]+)")) is not None:
        obj.offset = v
    if (v := _int(r"^\s*SAMPLE_BITS\s*=\s*(\d+)")) is not None:
        obj.sample_bits = v
    if (v := _int(r"^\s*BANDS\s*=\s*(\d+)")) is not None:
        obj.bands = v
    if m := re.search(r"SAMPLE_BIT_MASK\s*=\s*2#([01]+)#", text):
        obj.effective_max_dn = (1 << m.group(1).count("1")) - 1
    if m := re.search(r"FILTER_NAME\s*=\s*\(([^)]+)\)", text, re.DOTALL):
        names = [s.strip().strip('"').strip("'") for s in m.group(1).split(",")]
        if names:
            obj.filter_names = names
    elif m := re.search(r'FILTER_NAME\s*=\s*"?(\w[\w-]*)"?', text):
        obj.filter_names = [m.group(1)]

    return obj

MarsHiRISEBase

MarsHiRISEBase(root: Path, *, split: str = 'train', target: str | None = None, transforms: Callable[[Sample], Sample] | None = None, download: bool = False, bbox: tuple[float, float, float, float] | None = None, checksum: bool = False, reuse_cache: bool = True, res: float = 1.0 / 118502.26464032)

Bases: GeoDataset

Abstract base for Mars HiRISE-derived datasets.

Provides the common skeleton for index downloading/parsing, spatial filtering, parallel footprint extraction, async image downloading, coverage visualisation, and the verify/cache lifecycle.

Subclasses must implement the hooks marked @abstractmethod.

Source code in src/dataset/core/base.py

def __init__(
        self,
        root: Path,
        *,
        split: str = "train",
        target: str | None = None,
        transforms: Callable[[Sample], Sample] | None = None,
        download: bool = False,
        bbox: tuple[float, float, float, float] | None = None,
        checksum: bool = False,
        reuse_cache: bool = True,
        res: float = 1.0 / 118_502.26464032,
) -> None:
    super(GeoDataset, self).__init__()

    self.root = pathlib.Path(root)
    self.bbox = bbox
    self.split = split
    self.target = target
    self.transforms = transforms
    self.download = download
    self.checksum = checksum
    self.reuse_cache = reuse_cache

    self.res: float = res
    self._crs = MARS_GEOGRAPHIC_CRS

    self.index: gpd.GeoDataFrame | None = None
    self._raw_index: pd.DataFrame | None = None

    self._verify()

spatial_index_cache property

spatial_index_cache: Path

Path to the GeoPackage cache file.

Incorporates target, bbox, and subclass-specific suffix parts to avoid stale-cache collisions when configuration changes.

prefer_cog staticmethod

prefer_cog(path: Path | None) -> pathlib.Path | None

Return COG sidecar (.tif) for path if it exists, else path.

Source code in src/dataset/core/base.py

@staticmethod
def prefer_cog(path: pathlib.Path | None) -> pathlib.Path | None:
    """Return COG sidecar (.tif) for *path* if it exists, else *path*."""
    if path is None:
        return None
    cog = path.with_suffix(".tif")
    return cog if cog.exists() else path

merge_tiles staticmethod

merge_tiles(tiles: list[Tensor]) -> torch.Tensor

Mosaic co-registered tiles with a first-non-zero-wins strategy.

Source code in src/dataset/core/base.py

@staticmethod
def merge_tiles(tiles: list[torch.Tensor]) -> torch.Tensor:
    """Mosaic co-registered tiles with a first-non-zero-wins strategy."""
    if len(tiles) == 1:
        return tiles[0]

    max_h = max(t.shape[1] for t in tiles)
    max_w = max(t.shape[2] for t in tiles)
    n_ch = max(t.shape[0] for t in tiles)

    if not all(t.shape[0] == n_ch for t in tiles):
        logger.warning(
            "merge_tiles: channel mismatch (%d tiles, max %d ch) — padding.",
            len(tiles), n_ch,
        )
        padded = []
        for t in tiles:
            if t.shape[0] < n_ch:
                pad = torch.zeros(
                    n_ch - t.shape[0], t.shape[1], t.shape[2], dtype=t.dtype
                )
                t = torch.cat([t, pad])
            padded.append(t)
        tiles = padded

    merged = torch.zeros((n_ch, max_h, max_w), dtype=torch.float32)
    for tile in tiles:
        h, w = tile.shape[1], tile.shape[2]
        empty = merged[:, :h, :w] == 0.0
        merged[:, :h, :w][empty] = tile[empty]
    return merged

plot_coverage

plot_coverage(resolution: float | None = None, show_count: bool = True, suptitle: str | None = None) -> Figure

Visualise the spatial coverage of all entries in the index.

The plot is cropped to the actual extent of the entries, with a count heatmap plus per-entry bounding boxes.

Source code in src/dataset/core/base.py

def plot_coverage(
        self,
        resolution: float | None = None,
        show_count: bool = True,
        suptitle: str | None = None,
) -> Figure:
    """Visualise the spatial coverage of all entries in the index.

    The plot is cropped to the actual extent of the entries, with
    a count heatmap plus per-entry bounding boxes.
    """
    if self.index is None or len(self.index) == 0:
        raise RuntimeError("Spatial index is empty — nothing to plot.")

    all_bounds = self.index.geometry.bounds
    lon_min = float(all_bounds["minx"].min())
    lon_max = float(all_bounds["maxx"].max())
    lat_min = float(all_bounds["miny"].min())
    lat_max = float(all_bounds["maxy"].max())

    lon_span = lon_max - lon_min or 1.0
    lat_span = lat_max - lat_min or 1.0
    margin_lon = lon_span * 0.05
    margin_lat = lat_span * 0.05
    lon_min -= margin_lon
    lon_max += margin_lon
    lat_min -= margin_lat
    lat_max += margin_lat

    if resolution is None:
        resolution = min(lon_span, lat_span) / 20.0

    coverage, lon_edges, lat_edges = self._coverage_grid(
        resolution=resolution,
        lon_min=lon_min, lon_max=lon_max,
        lat_min=lat_min, lat_max=lat_max,
    )
    lon_centers = (lon_edges[:-1] + lon_edges[1:]) / 2
    lat_centers = (lat_edges[:-1] + lat_edges[1:]) / 2

    n_obs = len(self.index)
    covered_cells = int(np.count_nonzero(coverage))
    total_cells = coverage.size
    pct_covered = 100.0 * covered_cells / total_cells if total_cells else 0.0
    median_overlap = (
        float(np.median(coverage[coverage > 0])) if covered_cells else 0.0
    )
    max_overlap = int(coverage.max())

    fig = plt.figure(figsize=(16, 10))
    gs = fig.add_gridspec(
        2, 2,
        width_ratios=[4, 1], height_ratios=[1, 4],
        hspace=0.05, wspace=0.05,
    )
    ax_main = fig.add_subplot(gs[1, 0])
    ax_top = fig.add_subplot(gs[0, 0], sharex=ax_main)
    ax_right = fig.add_subplot(gs[1, 1], sharey=ax_main)

    display = np.ma.masked_equal(coverage, 0)
    norm = LogNorm(vmin=1, vmax=max(max_overlap, 1)) if show_count else None
    im = ax_main.pcolormesh(
        lon_edges, lat_edges, display,
        cmap="YlOrRd", norm=norm, shading="flat",
        rasterized=True, zorder=1,
    )
    ax_main.pcolormesh(
        lon_edges, lat_edges,
        np.ma.masked_not_equal(coverage, 0),
        cmap="Greys", vmin=0, vmax=1,
        shading="flat", rasterized=True, zorder=0,
    )

    cmap_boxes = matplotlib.colormaps.get_cmap("tab20")
    rects = []
    for geom in self.index.geometry:
        b = geom.bounds
        rects.append(Rectangle((b[0], b[1]), b[2] - b[0], b[3] - b[1]))

    box_colors = [cmap_boxes(i % 20) for i in range(n_obs)]
    pc = PatchCollection(
        rects,
        facecolors=[(r, g, b, 0.08) for r, g, b, _ in box_colors],
        edgecolors=[(r, g, b, 0.85) for r, g, b, _ in box_colors],
        linewidths=0.6, zorder=2,
    )
    ax_main.add_collection(pc)

    ax_main.set_xlim(lon_min, lon_max)
    ax_main.set_ylim(lat_min, lat_max)
    ax_main.set_xlabel("Longitude (°E, normalised to [−180, 180])")
    ax_main.set_ylabel("Latitude (°)")

    def _nice_ticks(lo: float, hi: float, n: int = 6) -> np.ndarray:
        span = hi - lo
        raw_step = span / n
        magnitude = 10 ** np.floor(np.log10(raw_step))
        for candidate in [1, 2, 2.5, 5, 10]:
            step = candidate * magnitude
            if span / step <= n + 1:
                break
        start = np.ceil(lo / step) * step
        return np.arange(start, hi + step * 0.5, step)

    ax_main.set_xticks(_nice_ticks(lon_min, lon_max))
    ax_main.set_yticks(_nice_ticks(lat_min, lat_max))
    ax_main.grid(color="white", linewidth=0.3, alpha=0.5, zorder=3)

    if show_count:
        cbar = fig.colorbar(im, ax=ax_main, pad=0.01, fraction=0.025)
        cbar.set_label("Entries per cell (log scale)", fontsize=8)

    legend_handles = [
        Patch(facecolor="#bbbbbb", label="No coverage"),
        Patch(facecolor=plt.cm.YlOrRd(0.15), label="Low overlap"),
        Patch(facecolor=plt.cm.YlOrRd(0.85), label="High overlap"),
    ]
    ax_main.legend(
        handles=legend_handles, loc="lower left",
        fontsize=7, framealpha=0.75,
    )

    ax_top.bar(
        lon_centers, coverage.sum(axis=0),
        width=resolution, color="steelblue", alpha=0.8, linewidth=0,
    )
    ax_top.set_ylabel("Sum", fontsize=8)
    ax_top.tick_params(labelbottom=False, labelsize=7)
    ax_top.grid(axis="y", linewidth=0.4, alpha=0.5)

    ax_right.barh(
        lat_centers, coverage.sum(axis=1),
        height=resolution, color="steelblue", alpha=0.8, linewidth=0,
    )
    ax_right.set_xlabel("Sum", fontsize=8)
    ax_right.tick_params(labelleft=False, labelsize=7)
    ax_right.grid(axis="x", linewidth=0.4, alpha=0.5)

    summary = (
        f"Entries : {n_obs:,}\n"
        f"Cells covered : {pct_covered:.1f}%  "
        f"({covered_cells:,} / {total_cells:,}  @  {resolution:.4f}°/cell)\n"
        f"Overlap — median : {median_overlap:.1f}×   max : {max_overlap}×"
    )
    ax_main.text(
        0.01, 0.99, summary,
        transform=ax_main.transAxes,
        va="top", ha="left", fontsize=8,
        bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.80),
        zorder=4,
    )

    title = suptitle or (
        f"{type(self).__name__} coverage — {n_obs:,} entries"
    )
    if self.target:
        title += f"  (filter: '{self.target}')"
    fig.suptitle(title, fontsize=12, y=1.005)

    return fig

plot_global_coverage

plot_global_coverage(output_pdf: str | Path, input_svg: str | Path = 'MarsTopography.svg', lat_bounds: tuple[float, float] = (-57.0, 57.0), target_id: str = 'topopgraphy') -> None

Injects dataset coordinates directly into an SVG and exports a vector PDF.

Handles SVGs where the map is embedded as a and instantiated via a tag.

Source code in src/dataset/core/base.py

def plot_global_coverage(self,  # pragma: no cover
                         output_pdf: str | pathlib.Path,
                         input_svg: str | pathlib.Path = "MarsTopography.svg",
                         lat_bounds: tuple[float, float] = (-57.0, 57.0),
                         target_id: str = "topopgraphy"  # Note: matching your snippet's spelling
                         ) -> None:
    """Injects dataset coordinates directly into an SVG and exports a vector PDF.

    Handles SVGs where the map is embedded as a <def> and instantiated via a <use> tag.
    """
    import xml.etree.ElementTree as ET
    import math
    import cairosvg

    # 1. Parse the SVG and register namespaces
    # Keeping namespaces empty prevents ugly ns0: prefixes in the output
    ET.register_namespace('', "http://www.w3.org/2000/svg")
    ET.register_namespace('xlink', "http://www.w3.org/1999/xlink")
    tree = ET.parse(input_svg)
    root = tree.getroot()

    # Build a parent map so we can easily inject nodes side-by-side
    parent_map = {c: p for p in tree.iter() for c in p}

    # 2. Locate the <use> node
    use_node = tree.find(f".//*[@id='{target_id}']")
    if use_node is None:
        raise ValueError(f"Could not find an SVG node with id='{target_id}'")

    # 3. Find the referenced definition
    # It could be under 'href' (SVG 2) or 'xlink:href' (SVG 1.1)
    href = use_node.attrib.get("{http://www.w3.org/1999/xlink}href") or use_node.attrib.get("href")
    if not href or not href.startswith("#"):
        raise ValueError(f"The node '{target_id}' does not have a valid href pointing to a def.")

    source_id = href[1:]  # Strip the '#'
    source_node = tree.find(f".//*[@id='{source_id}']")
    if source_node is None:
        raise ValueError(f"Could not find the referenced <image> node with id='{source_id}'")

    # 4. Extract raw dimensions from the source definition
    img_w = float(source_node.attrib['width'])
    img_h = float(source_node.attrib['height'])

    # Sometimes <use> elements have their own x/y offsets in addition to the transform
    use_x = float(use_node.attrib.get('x', 0.0))
    use_y = float(use_node.attrib.get('y', 0.0))

    # 5. Setup Projection Math
    min_lat, max_lat = lat_bounds

    def lat_to_mercator(lat_deg: float) -> float:
        lat_deg = max(min(lat_deg, 89.9), -89.9)
        return math.log(math.tan(math.pi / 4.0 + math.radians(lat_deg) / 2.0))

    merc_top = lat_to_mercator(max_lat)
    merc_bottom = lat_to_mercator(min_lat)
    merc_range = merc_top - merc_bottom

    def get_svg_coords(lon: float, lat: float) -> tuple[float, float]:
        """Map lon/lat to the raw image pixel coordinates."""
        norm_x = (lon + 180.0) / 360.0
        cx = use_x + (norm_x * img_w)

        merc_y = lat_to_mercator(lat)
        norm_y = (merc_top - merc_y) / merc_range
        cy = use_y + (norm_y * img_h)

        return cx, cy

    # 6. Create SVG Groups
    data_group = ET.Element('g', id="hirise_data_layers")

    # CRITICAL STEP: Copy the transform matrix from the <use> node to our data group
    # This forces the SVG renderer to align our dots with the scaled/moved image
    if 'transform' in use_node.attrib:
        data_group.set('transform', use_node.attrib['transform'])

    # Note: Because the transform matrix scales everything down (e.g., ~0.47x),
    # we need to increase the circle radius so they don't become microscopic.

    # Use "fill-opacity" instead of "opacity" to prevent CairoSVG from rasterizing the vectors.
    raw_group = ET.SubElement(data_group, 'g', id="layer_raw_index", fill="blue", stroke="black",
                              **{"fill-opacity": "0.7"})
    filtered_group = ET.SubElement(data_group, 'g', id="layer_filtered_index", fill="red", stroke="black",
                                   **{"stroke-width": "1.0"})

    df = self._load_cum_index()

    # 7. Inject the raw cumulative index
    if df is not None and not df.empty:
        raw_lons = ((df["MINIMUM_LONGITUDE"].astype(float) + 180.0) % 360.0) - 180.0
        raw_lats = df["MINIMUM_LATITUDE"].astype(float)

        for lon, lat in zip(raw_lons, raw_lats):
            if not (min_lat <= lat <= max_lat):
                continue
            cx, cy = get_svg_coords(lon, lat)
            ET.SubElement(raw_group, 'circle', cx=f"{cx:.2f}", cy=f"{cy:.2f}", r="2.5")

    # 8. Inject the precise filtered index points
    if self.index is not None and not self.index.empty:
        centroids = self.index.to_crs(MARS_MERCATOR_CRS).geometry.centroid.to_crs(MARS_GEOGRAPHIC_CRS)
        for lon, lat in zip(centroids.x, centroids.y):
            if not (min_lat <= lat <= max_lat):
                continue
            cx, cy = get_svg_coords(lon, lat)
            ET.SubElement(filtered_group, 'circle', cx=f"{cx:.2f}", cy=f"{cy:.2f}", r="8")

    # 9. Append the new elements right after the <use> node
    parent = parent_map[use_node]
    use_index = list(parent).index(use_node)
    parent.insert(use_index + 1, data_group)

    # 10. Save out to PDF
    import logging
    logger = logging.getLogger(__name__)

    if isinstance(output_pdf, str):
        output_pdf = Path(output_pdf)

    logger.info(f"Writing vectorized PDF to {output_pdf}")
    svg_bytes = ET.tostring(root, encoding='utf-8', method='xml')

    with open(output_pdf.with_suffix(".svg"), "wb") as f:
        f.write(svg_bytes)

    cairosvg.svg2pdf(bytestring=svg_bytes, write_to=str(output_pdf), dpi=600)

corners_to_polygon

corners_to_polygon(row: Series) -> Polygon | None

Build a Shapely Polygon from CORNER1–4 lat/lon columns in row.

Longitudes are normalised from PDS [0°, 360°] to [−180°, 180°]. Returns None on missing/degenerate data.

Source code in src/dataset/core/base.py

def corners_to_polygon(row: pd.Series) -> Polygon | None:
    """Build a Shapely Polygon from CORNER1–4 lat/lon columns in *row*.

    Longitudes are normalised from PDS [0°, 360°] to [−180°, 180°].
    Returns ``None`` on missing/degenerate data.
    """
    try:
        coords = [
            (
                ((float(row[f"CORNER{i}_LONGITUDE"]) + 180.0) % 360.0) - 180.0,
                float(row[f"CORNER{i}_LATITUDE"]),
            )
            for i in (1, 2, 3, 4)
        ]
    except (KeyError, ValueError, TypeError):
        return None

    if any(math.isnan(lon) or math.isnan(lat) for lon, lat in coords):
        return None

    poly = Polygon(coords)
    if not poly.is_valid:
        poly = poly.buffer(0)
    return poly if (poly.is_valid and not poly.is_empty) else None

extract_footprint

extract_footprint(file_path: str | PathLike | None, mars_crs: CRS = MARS_GEOGRAPHIC_CRS, nodata_test: Callable | None = None) -> tuple[list[tuple[float, float]] | None, tuple[float, float, float, float] | None]

Extract the convex hull of valid pixels from a raster.

Reads band 1 at the coarsest available overview level so even multi-GB images resolve to a few hundred pixels.

Thread-safe: each call opens its own file handle.

Parameters:

Name	Type	Description	Default
file_path	str | PathLike | None	Path to the raster file (JP2, IMG, or COG).	required
mars_crs	CRS	Target geographic CRS to reproject vertices into.	MARS_GEOGRAPHIC_CRS
nodata_test	Callable | None	Optional callable (data_array) -> bool_mask that returns a boolean mask of valid pixels. Defaults to data > 0.	None

Returns:

Type	Description
list[tuple[float, float]] | None	(hull_coords, file_bounds) — hull_coords is a list of
tuple[float, float, float, float] | None	(lon, lat) vertices, or None on failure. file_bounds is a
tuple[list[tuple[float, float]] | None, tuple[float, float, float, float] | None]	(west, south, east, north) fallback.

Source code in src/dataset/core/base.py

def extract_footprint(
        file_path: str | PathLike | None,
        mars_crs: rasterio.crs.CRS = MARS_GEOGRAPHIC_CRS,
        nodata_test: Callable | None = None,
) -> tuple[list[tuple[float, float]] | None, tuple[float, float, float, float] | None]:
    """Extract the convex hull of valid pixels from a raster.

    Reads band 1 at the coarsest available overview level so even multi-GB
    images resolve to a few hundred pixels.

    Thread-safe: each call opens its own file handle.

    Args:
        file_path: Path to the raster file (JP2, IMG, or COG).
        mars_crs: Target geographic CRS to reproject vertices into.
        nodata_test: Optional callable ``(data_array) -> bool_mask`` that
            returns a boolean mask of *valid* pixels.  Defaults to
            ``data > 0``.

    Returns:
        ``(hull_coords, file_bounds)`` — hull_coords is a list of
        ``(lon, lat)`` vertices, or ``None`` on failure.  file_bounds is a
        ``(west, south, east, north)`` fallback.
    """
    from rasterio.warp import transform as warp_transform

    if file_path is None:
        return None, None

    if isinstance(file_path, str):
        path = pathlib.Path(file_path)
    else:
        path = file_path

    cog = path.with_suffix(".tif")
    actual = cog if cog.exists() else path
    if not actual.exists():
        return None, None

    try:
        with rasterio.open(actual) as src:
            src_crs = src.crs
            if src_crs is None:
                return None, None

            # File bounds (always — cheap fallback)
            try:
                fl, fb, fr, ft = transform_bounds(src_crs, mars_crs, *src.bounds)
                fl = ((fl + 180.0) % 360.0) - 180.0
                fr = ((fr + 180.0) % 360.0) - 180.0
                if not (-180.0 <= fl < fr <= 180.0 and -90.0 <= fb < ft <= 90.0):
                    file_bounds = None
                else:
                    file_bounds = (fl, fb, fr, ft)
            except Exception:
                file_bounds = None

            # Read band 1 at coarsest overview
            ovrs = src.overviews(1)
            factor = max(ovrs) if ovrs else max(1, min(src.height, src.width) // 500)
            oh = max(1, src.height // factor)
            ow = max(1, src.width // factor)
            data = src.read(1, out_shape=(oh, ow))

            if nodata_test is not None:
                valid_mask = nodata_test(data)
            else:
                valid_mask = data > 0

            ys, xs = np.where(valid_mask)
            if len(xs) < 3:
                return None, file_bounds

            step = max(1, len(xs) // 4000)
            hull = MultiPoint(
                list(zip(xs[::step].tolist(), ys[::step].tolist()))
            ).convex_hull
            if hull.is_empty:
                return None, file_bounds

            hull_px = np.array(hull.exterior.coords)
            ovr_tf = rasterio.transform.from_bounds(*src.bounds, ow, oh)
            src_xs, src_ys = rasterio.transform.xy(
                ovr_tf, hull_px[:, 1].tolist(), hull_px[:, 0].tolist()
            )
            geo_xs, geo_ys = warp_transform(
                src_crs, mars_crs, list(src_xs), list(src_ys)
            )
            geo_xs = [((x + 180.0) % 360.0) - 180.0 for x in geo_xs]

            return list(zip(geo_xs, geo_ys)), file_bounds

    except Exception:
        return None, None

reproject_band

reproject_band(src_dataset: DatasetReader, band_idx: int, dst_crs: CRS, dst_transform: Affine, out_h: int, out_w: int, *, src_nodata: float | None = None, dst_nodata: float = 0.0, resampling: Resampling = Resampling.bilinear) -> np.ndarray

Reproject a single band from an open rasterio dataset.

Returns an (out_h, out_w) float32 array.

Source code in src/dataset/core/base.py

def reproject_band(
        src_dataset: rasterio.DatasetReader,
        band_idx: int,
        dst_crs: rasterio.crs.CRS,
        dst_transform: rasterio.Affine,
        out_h: int,
        out_w: int,
        *,
        src_nodata: float | None = None,
        dst_nodata: float = 0.0,
        resampling: Resampling = Resampling.bilinear,
) -> np.ndarray:
    """Reproject a single band from an open rasterio dataset.

    Returns an ``(out_h, out_w)`` float32 array.
    """
    src_crs = src_dataset.crs
    if src_crs is None:
        src_crs = dst_crs

    dest = np.empty((out_h, out_w), dtype=np.float32)

    kwargs: dict = dict(
        source=rasterio.band(src_dataset, band_idx),
        destination=dest,
        src_transform=src_dataset.transform,
        src_crs=src_crs,
        dst_transform=dst_transform,
        dst_crs=dst_crs,
        resampling=resampling,
        dst_nodata=dst_nodata,
    )
    if src_nodata is not None:
        kwargs["src_nodata"] = src_nodata

    reproject(**kwargs)
    return dest

check_overlap

check_overlap(src_dataset: DatasetReader, dst_crs: CRS, x: slice, y: slice, tol: float = _SPATIAL_TOL) -> bool

Return True if the file bounds overlap the query window.

Absorbs floating-point imprecision and normalises longitudes.

Source code in src/dataset/core/base.py

def check_overlap(
        src_dataset: rasterio.DatasetReader,
        dst_crs: rasterio.crs.CRS,
        x: slice,
        y: slice,
        tol: float = _SPATIAL_TOL,
) -> bool:
    """Return True if the file bounds overlap the query window.

    Absorbs floating-point imprecision and normalises longitudes.
    """
    src_crs = src_dataset.crs
    if src_crs is None:
        return True  # can't check — assume overlap
    try:
        fl, fb, fr, ft = transform_bounds(src_crs, dst_crs, *src_dataset.bounds)
        fl = ((fl + 180.0) % 360.0) - 180.0
        fr = ((fr + 180.0) % 360.0) - 180.0
        if fl <= fr:
            if (
                    fr + tol < x.start
                    or fl - tol > x.stop
                    or ft + tol < y.start
                    or fb - tol > y.stop
            ):
                return False
    except Exception:
        pass
    return True