Skip to content

Tessellation

Tessellation divides geographic areas of interest (AOIs) into smaller rectangular tiles. These tiles are normally sized to match the sensor field-of-view for Earth-imaging satellites. This enables larger areas, ones too big to be collected in a single imaging action, to be broken down into smaller parts which can be feasibly collected.

There are infinitely many ways to tile a large area if entirely unconstrained in the tile placement. Brahe implements an orbit-geometry based tessalator that generates tiles aligned with the orbital ground-track of a satellite. This approach is particular well-suited to satellites with push-broom imaging modes such as radar imaging satellites. The OrbitGeometryTessellator uses a satellite's orbital elements and a reference epoch to determine ground-track directions at any latitude. It then tiles the target location perpendicular and parallel to the ground track. Output tiles are PolygonLocation instances with metadata properties describing the tile geometry, making them compatible with the rest of the access computation system. The tesselation configuration should be setup such that the maximum width and length remain feasible to collect in a single imaging pass.

For complete API details, see the API Reference: Tessellation.

Configuration

The OrbitGeometryTessellatorConfig controls tile dimensions, overlap, and ascending/descending pass selection. All dimensions are in meters.

Parameter Default Description
image_width 5000 m Cross-track tile width
image_length 5000 m Along-track tile length
crosstrack_overlap 200 m Cross-track overlap between adjacent strips
alongtrack_overlap 200 m Along-track overlap between adjacent tiles
asc_dsc Either Ascending/descending pass selection
min_image_length 5000 m Minimum tile length (tiles shorter than this are discarded)
max_image_length 5000 m Maximum tile length (tiles longer than this are split)
bh.initialize_eop()

# Default configuration
config = bh.OrbitGeometryTessellatorConfig()
print(f"Default image_width: {config.image_width} m")
print(f"Default image_length: {config.image_length} m")
print(f"Default crosstrack_overlap: {config.crosstrack_overlap} m")
print(f"Default alongtrack_overlap: {config.alongtrack_overlap} m")
print(f"Default min_image_length: {config.min_image_length} m")
print(f"Default max_image_length: {config.max_image_length} m")

# Custom configuration for ascending passes with larger tiles
custom_config = bh.OrbitGeometryTessellatorConfig(
    image_width=10000,
    image_length=15000,
    crosstrack_overlap=300,
    alongtrack_overlap=300,
    asc_dsc=bh.AscDsc.ASCENDING,
    min_image_length=8000,
    max_image_length=25000,
)
print(f"\nCustom image_width: {custom_config.image_width} m")
print(f"Custom image_length: {custom_config.image_length} m")
print(f"Custom asc_dsc: {custom_config.asc_dsc}")
use brahe as bh;
use bh::access::constraints::AscDsc;
use bh::access::tessellation::OrbitGeometryTessellatorConfig;

fn main() {
    bh::initialize_eop().unwrap();

    // Default configuration
    let config = OrbitGeometryTessellatorConfig::default();
    println!("Default image_width: {} m", config.image_width);
    println!("Default image_length: {} m", config.image_length);
    println!("Default crosstrack_overlap: {} m", config.crosstrack_overlap);
    println!("Default alongtrack_overlap: {} m", config.alongtrack_overlap);
    println!("Default min_image_length: {} m", config.min_image_length);
    println!("Default max_image_length: {} m", config.max_image_length);

    // Custom configuration for ascending passes with larger tiles
    let custom_config = OrbitGeometryTessellatorConfig::new(10000.0, 15000.0)
        .with_crosstrack_overlap(300.0)
        .with_alongtrack_overlap(300.0)
        .with_asc_dsc(AscDsc::Ascending)
        .with_min_image_length(8000.0)
        .with_max_image_length(25000.0);
    println!("\nCustom image_width: {} m", custom_config.image_width);
    println!("Custom image_length: {} m", custom_config.image_length);
    println!("Custom asc_dsc: {:?}", custom_config.asc_dsc);

}
Output
Default image_width: 5000.0 m
Default image_length: 5000.0 m
Default crosstrack_overlap: 200.0 m
Default alongtrack_overlap: 200.0 m
Default min_image_length: 5000.0 m
Default max_image_length: 5000.0 m

Custom image_width: 10000.0 m
Custom image_length: 15000.0 m
Custom asc_dsc: Ascending
Default image_width: 5000 m
Default image_length: 5000 m
Default crosstrack_overlap: 200 m
Default alongtrack_overlap: 200 m
Default min_image_length: 5000 m
Default max_image_length: 5000 m

Custom image_width: 10000 m
Custom image_length: 15000 m
Custom asc_dsc: Ascending

Point Tessellation

Tessellating a PointLocation creates one tile per pass direction, centered on the point. With AscDsc.ASCENDING, a single tile is created; with AscDsc.EITHER, up to two tiles are created (one per direction). At high latitudes where ascending and descending ground tracks converge, redundant tiles may be automatically merged.

bh.initialize_eop()

# ISS TLE
line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927"
line2 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537"

# Create propagator and tessellator
prop = bh.SGPPropagator.from_tle(line1, line2, step_size=60.0)
config = bh.OrbitGeometryTessellatorConfig(
    image_width=5000,
    image_length=5000,
    asc_dsc=bh.AscDsc.ASCENDING,
)
tess = bh.OrbitGeometryTessellator(prop, prop.epoch, config, spacecraft_id="ISS")

# Tessellate a point
point = bh.PointLocation(10.0, 30.0, 0.0)
tiles = tess.tessellate_point(point)

print(f"Number of tiles: {len(tiles)}")
for i, tile in enumerate(tiles):
    props = tile.properties
    center = tile.center_geodetic()
    print(f"Tile {i}: center=({center[0]:.4f}, {center[1]:.4f})")
    print(f"  width={props['tile_width']:.0f} m, length={props['tile_length']:.0f} m")
use brahe as bh;
use bh::access::constraints::AscDsc;
use bh::access::location::{AccessibleLocation, PointLocation};
use bh::access::tessellation::{OrbitGeometryTessellator, OrbitGeometryTessellatorConfig, Tessellator};
use bh::propagators::SGPPropagator;

fn main() {
    bh::initialize_eop().unwrap();

    // ISS TLE
    let line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927";
    let line2 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537";

    // Create propagator and tessellator
    let prop = SGPPropagator::from_tle(line1, line2, 60.0).unwrap();
    let epoch = prop.epoch;
    let config = OrbitGeometryTessellatorConfig::default()
        .with_asc_dsc(AscDsc::Ascending);
    let tess = OrbitGeometryTessellator::new(
        Box::new(prop), epoch, config, Some("ISS".to_string()),
    );

    // Tessellate a point
    let point = PointLocation::new(10.0, 30.0, 0.0);
    let tiles = tess.tessellate(&point).unwrap();

    println!("Number of tiles: {}", tiles.len());
    for (i, tile) in tiles.iter().enumerate() {
        let center = tile.center_geodetic();
        let props = tile.properties();
        println!("Tile {}: center=({:.4}, {:.4})", i, center[0], center[1]);
        println!("  width={:.0} m, length={:.0} m",
            props["tile_width"].as_f64().unwrap(),
            props["tile_length"].as_f64().unwrap(),
        );
    }

}
Output
1
2
3
Number of tiles: 1
Tile 0: center=(10.0000, 29.9999)
  width=5000 m, length=5000 m
1
2
3
Number of tiles: 1
Tile 0: center=(10.0000, 29.9999)
  width=5000 m, length=5000 m

The figure below shows the difference between ascending-only and ascending+descending tessellation for a single point near San Francisco. Each tile direction produces a rectangle aligned to the satellite ground track at that latitude.

Point tessellation Point tessellation
Figure Source
        # Figure 1: Point Tessellation — ascending-only vs ascending+descending
        # ------------------------------------------------------------------

        fig, (ax1, ax2) = plt.subplots(
            1, 2, figsize=(12, 5), subplot_kw={"projection": ccrs.PlateCarree()}
        )
        if theme == "dark":
            set_dark_figure_bg(fig)

        extent = [-122.5, -122.3, 37.72, 37.88]
        for ax in (ax1, ax2):
            ax.set_extent(extent, crs=ccrs.PlateCarree())
            style_map_axis(ax, theme)
            ax.plot(
                -122.4,
                37.8,
                "*",
                color=marker_color,
                markersize=12,
                transform=ccrs.PlateCarree(),
                zorder=5,
            )

        draw_tiles(ax1, tiles_asc)
        ax1.set_title(f"Ascending only ({len(tiles_asc)} tile)", fontsize=10)

        draw_tiles(ax2, tiles_either)
        ax2.set_title(
            f"Ascending + Descending ({len(tiles_either)} tiles)", fontsize=10
        )

        for ax in (ax1, ax2):
            style_gridlines(ax, theme)

        plt.tight_layout()
        save_themed(fig, "tessellation_point", theme)

Polygon Tessellation

Tessellating a PolygonLocation divides the area into cross-track strips perpendicular to the satellite ground track, then subdivides each strip along-track into individual tiles. The algorithm handles concave polygons by detecting gaps in the along-track direction. Tiles at polygon edges may have adjusted lengths to fit the boundary.

import brahe as bh

bh.initialize_eop()

# ISS TLE
line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927"
line2 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537"

# Create propagator and tessellator
prop = bh.SGPPropagator.from_tle(line1, line2, step_size=60.0)
config = bh.OrbitGeometryTessellatorConfig(
    image_width=5000,
    image_length=5000,
    asc_dsc=bh.AscDsc.ASCENDING,
)
tess = bh.OrbitGeometryTessellator(prop, prop.epoch, config, spacecraft_id="ISS")

# Define a small polygon (approximately 0.1 deg x 0.1 deg)
vertices = np.array(
    [
        [10.0, 30.0, 0.0],
        [10.1, 30.0, 0.0],
        [10.1, 30.1, 0.0],
        [10.0, 30.1, 0.0],
    ]
)
polygon = bh.PolygonLocation(vertices)

# Tessellate the polygon
tiles = tess.tessellate_polygon(polygon)

print(f"Number of tiles: {len(tiles)}")
for i, tile in enumerate(tiles):
    props = tile.properties
    print(
        f"Tile {i}: group_id={props['tile_group_id'][:8]}... "
        f"length={props['tile_length']:.0f} m"
    )
use brahe as bh;
use bh::access::constraints::AscDsc;
use bh::access::location::{AccessibleLocation, PolygonLocation};
use bh::access::tessellation::{OrbitGeometryTessellator, OrbitGeometryTessellatorConfig, Tessellator};
use bh::propagators::SGPPropagator;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // ISS TLE
    let line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927";
    let line2 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537";

    // Create propagator and tessellator
    let prop = SGPPropagator::from_tle(line1, line2, 60.0).unwrap();
    let epoch = prop.epoch;
    let config = OrbitGeometryTessellatorConfig::default()
        .with_asc_dsc(AscDsc::Ascending);
    let tess = OrbitGeometryTessellator::new(
        Box::new(prop), epoch, config, Some("ISS".to_string()),
    );

    // Define a small polygon (approximately 0.1 deg x 0.1 deg)
    let vertices = vec![
        na::SVector::<f64, 3>::new(10.0, 30.0, 0.0),
        na::SVector::<f64, 3>::new(10.1, 30.0, 0.0),
        na::SVector::<f64, 3>::new(10.1, 30.1, 0.0),
        na::SVector::<f64, 3>::new(10.0, 30.1, 0.0),
    ];
    let polygon = PolygonLocation::new(vertices).unwrap();

    // Tessellate the polygon
    let tiles = tess.tessellate(&polygon).unwrap();

    println!("Number of tiles: {}", tiles.len());
    for (i, tile) in tiles.iter().enumerate() {
        let props = tile.properties();
        let group_id = props["tile_group_id"].as_str().unwrap();
        println!("Tile {}: group_id={}... length={:.0} m",
            i, &group_id[..8],
            props["tile_length"].as_f64().unwrap(),
        );
    }
}
Output
Number of tiles: 12
Tile 0: group_id=bc66a83b... length=5000 m
Tile 1: group_id=bc66a83b... length=5000 m
Tile 2: group_id=bc66a83b... length=5000 m
Tile 3: group_id=bc66a83b... length=5000 m
Tile 4: group_id=bc66a83b... length=5000 m
Tile 5: group_id=bc66a83b... length=5000 m
Tile 6: group_id=bc66a83b... length=5000 m
Tile 7: group_id=bc66a83b... length=5000 m
Tile 8: group_id=bc66a83b... length=5000 m
Tile 9: group_id=bc66a83b... length=5000 m
Tile 10: group_id=bc66a83b... length=5000 m
Tile 11: group_id=bc66a83b... length=5000 m
Number of tiles: 12
Tile 0: group_id=4f1a8f93... length=5000 m
Tile 1: group_id=4f1a8f93... length=5000 m
Tile 2: group_id=4f1a8f93... length=5000 m
Tile 3: group_id=4f1a8f93... length=5000 m
Tile 4: group_id=4f1a8f93... length=5000 m
Tile 5: group_id=4f1a8f93... length=5000 m
Tile 6: group_id=4f1a8f93... length=5000 m
Tile 7: group_id=4f1a8f93... length=5000 m
Tile 8: group_id=4f1a8f93... length=5000 m
Tile 9: group_id=4f1a8f93... length=5000 m
Tile 10: group_id=4f1a8f93... length=5000 m
Tile 11: group_id=4f1a8f93... length=5000 m

The figure below shows England tessellated with 50 km tiles. Tiles are colored by tile_group_id — each color represents tiles sharing the same ground-track direction (ascending vs descending). The dashed line is the input polygon boundary.

Polygon tessellation Polygon tessellation
Figure Source
        # Figure 2: Polygon Tessellation — England with tiles colored by group
        # ------------------------------------------------------------------

        fig, ax = plt.subplots(
            1, 1, figsize=(8, 7), subplot_kw={"projection": ccrs.PlateCarree()}
        )
        if theme == "dark":
            set_dark_figure_bg(fig)

        ax.set_extent([-7, 3, 49, 56], crs=ccrs.PlateCarree())
        style_map_axis(ax, theme)
        draw_tiles(ax, tiles_eng)
        draw_polygon_outline(ax, england_verts, color=outline_color)
        ax.set_title(
            f"England — {len(tiles_eng)} tiles, colored by tile_group_id", fontsize=10
        )
        style_gridlines(ax, theme)
        plt.tight_layout()
        save_themed(fig, "tessellation_polygon", theme)

Effect of Configuration Parameters

Tile Length

Increasing image_width and image_length produces fewer, larger tiles. The left panel uses 5 km tiles and the right uses 15 km tiles for the same region near San Francisco.

Tile length comparison Tile length comparison

Overlap

Increasing crosstrack_overlap and alongtrack_overlap causes adjacent tiles to share more area, which produces more tiles for the same region. The left panel uses 0 m overlap; the right uses 1000 m overlap.

Overlap comparison Overlap comparison
Figure Source
        # Figure 3: Tile Length — 5 km vs 15 km
        # Figure 4: Overlap — 0 m vs 1000 m
        # ------------------------------------------------------------------

        # Tile Length
        fig, (ax1, ax2) = plt.subplots(
            1, 2, figsize=(14, 6), subplot_kw={"projection": ccrs.PlateCarree()}
        )
        if theme == "dark":
            set_dark_figure_bg(fig)

        for ax in (ax1, ax2):
            ax.set_extent(sf_extent, crs=ccrs.PlateCarree())
            style_map_axis(ax, theme)

        draw_tiles(ax1, tiles_5k)
        draw_polygon_outline(ax1, sf_verts, color=outline_color)
        ax1.set_title(f"5 km tiles ({len(tiles_5k)} tiles)", fontsize=10)

        draw_tiles(ax2, tiles_15k)
        draw_polygon_outline(ax2, sf_verts, color=outline_color)
        ax2.set_title(f"15 km tiles ({len(tiles_15k)} tiles)", fontsize=10)

        for ax in (ax1, ax2):
            style_gridlines(ax, theme)

        plt.tight_layout()
        save_themed(fig, "tessellation_tile_length", theme)

        # Overlap
        fig, (ax1, ax2) = plt.subplots(
            1, 2, figsize=(14, 6), subplot_kw={"projection": ccrs.PlateCarree()}
        )
        if theme == "dark":
            set_dark_figure_bg(fig)

        for ax in (ax1, ax2):
            ax.set_extent(sf_extent, crs=ccrs.PlateCarree())
            style_map_axis(ax, theme)

        draw_tiles(ax1, tiles_no_ol)
        draw_polygon_outline(ax1, sf_verts, color=outline_color)
        ax1.set_title(f"0 m overlap ({len(tiles_no_ol)} tiles)", fontsize=10)

        draw_tiles(ax2, tiles_ol)
        draw_polygon_outline(ax2, sf_verts, color=outline_color)
        ax2.set_title(f"1000 m overlap ({len(tiles_ol)} tiles)", fontsize=10)

        for ax in (ax1, ax2):
            style_gridlines(ax, theme)

        plt.tight_layout()
        save_themed(fig, "tessellation_overlap", theme)

Tile Metadata Properties

Each output tile is a PolygonLocation with metadata properties stored in its properties dictionary. These properties describe the tile geometry and ownership.

Property Type Description
tile_direction [x, y, z] Along-track unit vector in ECEF coordinates
tile_width float Cross-track dimension in meters
tile_length float Along-track dimension in meters
tile_area float Tile area (\(\text{width} \times \text{length}\)) in m\(^2\)
tile_group_id str UUID shared by all tiles in the same tiling direction
spacecraft_ids list[str] Spacecraft identifiers that can collect this tile
import brahe as bh

bh.initialize_eop()

# ISS TLE
line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927"
line2 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537"

# Create propagator and tessellator
prop = bh.SGPPropagator.from_tle(line1, line2, step_size=60.0)
config = bh.OrbitGeometryTessellatorConfig(
    image_width=5000,
    image_length=5000,
    asc_dsc=bh.AscDsc.ASCENDING,
)
tess = bh.OrbitGeometryTessellator(prop, prop.epoch, config, spacecraft_id="ISS")

# Tessellate a point and inspect properties
point = bh.PointLocation(10.0, 30.0, 0.0)
tiles = tess.tessellate_point(point)
tile = tiles[0]
props = tile.properties

# Along-track direction (unit vector in ECEF)
direction = np.array(props["tile_direction"])
print(f"tile_direction: [{direction[0]:.4f}, {direction[1]:.4f}, {direction[2]:.4f}]")
print(f"  magnitude: {np.linalg.norm(direction):.6f}")

# Tile dimensions
print(f"tile_width: {props['tile_width']:.0f} m")
print(f"tile_length: {props['tile_length']:.0f} m")
print(f"tile_area: {props['tile_area']:.0f} m^2")

# Group and spacecraft identifiers
print(f"tile_group_id: {props['tile_group_id'][:8]}...")
print(f"spacecraft_ids: {props['spacecraft_ids']}")
use brahe as bh;
use bh::access::constraints::AscDsc;
use bh::access::location::{AccessibleLocation, PointLocation};
use bh::access::tessellation::{OrbitGeometryTessellator, OrbitGeometryTessellatorConfig, Tessellator};
use bh::propagators::SGPPropagator;

fn main() {
    bh::initialize_eop().unwrap();

    // ISS TLE
    let line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927";
    let line2 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537";

    // Create propagator and tessellator
    let prop = SGPPropagator::from_tle(line1, line2, 60.0).unwrap();
    let epoch = prop.epoch;
    let config = OrbitGeometryTessellatorConfig::default()
        .with_asc_dsc(AscDsc::Ascending);
    let tess = OrbitGeometryTessellator::new(
        Box::new(prop), epoch, config, Some("ISS".to_string()),
    );

    // Tessellate a point and inspect properties
    let point = PointLocation::new(10.0, 30.0, 0.0);
    let tiles = tess.tessellate(&point).unwrap();
    let props = tiles[0].properties();

    // Along-track direction (unit vector in ECEF)
    let dir = props["tile_direction"].as_array().unwrap();
    let dx = dir[0].as_f64().unwrap();
    let dy = dir[1].as_f64().unwrap();
    let dz = dir[2].as_f64().unwrap();
    let mag = (dx * dx + dy * dy + dz * dz).sqrt();
    println!("tile_direction: [{:.4}, {:.4}, {:.4}]", dx, dy, dz);
    println!("  magnitude: {:.6}", mag);

    // Tile dimensions
    println!("tile_width: {:.0} m", props["tile_width"].as_f64().unwrap());
    println!("tile_length: {:.0} m", props["tile_length"].as_f64().unwrap());
    println!("tile_area: {:.0} m^2", props["tile_area"].as_f64().unwrap());

    // Group and spacecraft identifiers
    let group_id = props["tile_group_id"].as_str().unwrap();
    println!("tile_group_id: {}...", &group_id[..8]);
    println!("spacecraft_ids: {}", props["spacecraft_ids"]);
}
Output
1
2
3
4
5
6
7
tile_direction: [-0.4755, 0.6135, 0.6305]
  magnitude: 1.000000
tile_width: 5000 m
tile_length: 5000 m
tile_area: 25000000 m^2
tile_group_id: b6e590c4...
spacecraft_ids: ['ISS']
1
2
3
4
5
6
7
tile_direction: [-0.4755, 0.6135, 0.6305]
  magnitude: 1.000000
tile_width: 5000 m
tile_length: 5000 m
tile_area: 25000000 m^2
tile_group_id: 19790c42...
spacecraft_ids: ["ISS"]

Merging Tiles from Multiple Spacecraft

When multiple spacecraft have similar orbital planes, their ground-track directions at a given latitude will be similar. The tile_merge_orbit_geometry function clusters tiles by direction and merges groups whose directions fall within a configurable angular threshold. Rather than creating duplicate tiles, it adds the additional spacecraft's ID to the base tile's spacecraft_ids list.

bh.initialize_eop()

# SC-1 and SC-2 TLEs with slightly different inclinations (~1.4 degree offset)
line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927"
line2_sc1 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537"
line2_sc2 = "2 25544  53.0000 247.4627 0006703 130.5360 325.0288 15.72125391563532"

# Create two tessellators with different spacecraft IDs
config = bh.OrbitGeometryTessellatorConfig(
    image_width=5000,
    image_length=5000,
    asc_dsc=bh.AscDsc.ASCENDING,
)

prop1 = bh.SGPPropagator.from_tle(line1, line2_sc1, step_size=60.0)
tess1 = bh.OrbitGeometryTessellator(prop1, prop1.epoch, config, spacecraft_id="SC-1")

prop2 = bh.SGPPropagator.from_tle(line1, line2_sc2, step_size=60.0)
tess2 = bh.OrbitGeometryTessellator(prop2, prop2.epoch, config, spacecraft_id="SC-2")

# Tessellate the same point with both spacecraft
point = bh.PointLocation(10.0, 30.0, 0.0)
tiles_sc1 = tess1.tessellate_point(point)
tiles_sc2 = tess2.tessellate_point(point)
all_tiles = tiles_sc1 + tiles_sc2

print(f"Before merge: {len(all_tiles)} tiles")

# Merge tiles with similar directions
merged = bh.tile_merge_orbit_geometry(all_tiles, 200.0, 200.0, 2.0)

print(f"After merge: {len(merged)} tiles")
for tile in merged:
    print(f"  spacecraft_ids: {tile.properties['spacecraft_ids']}")
use brahe as bh;
use bh::access::constraints::AscDsc;
use bh::access::location::{AccessibleLocation, PointLocation};
use bh::access::tessellation::{
    OrbitGeometryTessellator, OrbitGeometryTessellatorConfig, Tessellator,
    tile_merge_orbit_geometry,
};
use bh::propagators::SGPPropagator;

fn main() {
    bh::initialize_eop().unwrap();

    // SC-1 and SC-2 TLEs with slightly different inclinations (~1.4 degree offset)
    let line1 = "1 25544U 98067A   08264.51782528 -.00002182  00000-0 -11606-4 0  2927";
    let line2_sc1 = "2 25544  51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537";
    let line2_sc2 = "2 25544  53.0000 247.4627 0006703 130.5360 325.0288 15.72125391563532";

    // Create two tessellators with different spacecraft IDs
    let config = OrbitGeometryTessellatorConfig::default()
        .with_asc_dsc(AscDsc::Ascending);

    let prop1 = SGPPropagator::from_tle(line1, line2_sc1, 60.0).unwrap();
    let epoch1 = prop1.epoch;
    let tess1 = OrbitGeometryTessellator::new(
        Box::new(prop1), epoch1, config.clone(), Some("SC-1".to_string()),
    );

    let prop2 = SGPPropagator::from_tle(line1, line2_sc2, 60.0).unwrap();
    let epoch2 = prop2.epoch;
    let tess2 = OrbitGeometryTessellator::new(
        Box::new(prop2), epoch2, config, Some("SC-2".to_string()),
    );

    // Tessellate the same point with both spacecraft
    let point = PointLocation::new(10.0, 30.0, 0.0);
    let tiles_sc1 = tess1.tessellate(&point).unwrap();
    let tiles_sc2 = tess2.tessellate(&point).unwrap();
    let mut all_tiles = tiles_sc1;
    all_tiles.extend(tiles_sc2);

    println!("Before merge: {} tiles", all_tiles.len());

    // Merge tiles with similar directions
    let merged = tile_merge_orbit_geometry(&all_tiles, 200.0, 200.0, 2.0);

    println!("After merge: {} tiles", merged.len());
    for tile in &merged {
        println!("  spacecraft_ids: {}", tile.properties()["spacecraft_ids"]);
    }
}
Output
1
2
3
Before merge: 2 tiles
After merge: 1 tiles
  spacecraft_ids: ['SC-1', 'SC-2']
1
2
3
Before merge: 2 tiles
After merge: 1 tiles
  spacecraft_ids: ["SC-1","SC-2"]

The figure below shows tiles from two spacecraft with slightly different inclinations (~1.4° offset). Before merging, the tiles from SC-1 and SC-2 are visibly offset; after merging with a 2° angular threshold, overlapping tiles are combined with both spacecraft IDs in the spacecraft_ids list.

Merging before/after Merging before/after
Figure Source
        # Figure 5: Merging — before/after merge from two spacecraft
        # SC-2 has ~1.4° inclination offset so tiles visibly differ
        # ------------------------------------------------------------------

        sc1_color = (0.12, 0.47, 0.71)  # blue
        sc2_color = (1.0, 0.50, 0.05)  # orange
        merged_color = (0.17, 0.63, 0.17)  # green

        fig, (ax1, ax2) = plt.subplots(
            1, 2, figsize=(14, 6), subplot_kw={"projection": ccrs.PlateCarree()}
        )
        if theme == "dark":
            set_dark_figure_bg(fig)

        for ax in (ax1, ax2):
            ax.set_extent(sf_extent, crs=ccrs.PlateCarree())
            style_map_axis(ax, theme)

        # Before: SC-1 in blue, SC-2 in orange
        draw_tiles(ax1, tiles_sc1, color_by_group=False, color_cycle=[sc1_color])
        draw_tiles(ax1, tiles_sc2, color_by_group=False, color_cycle=[sc2_color])
        draw_polygon_outline(ax1, sf_verts, color=outline_color)
        ax1.set_title(f"Before merge ({len(tiles_before)} tiles)", fontsize=10)
        ax1.plot([], [], "s", color=sc1_color, label="SC-1", markersize=8)
        ax1.plot([], [], "s", color=sc2_color, label="SC-2", markersize=8)
        ax1.legend(loc="lower right", fontsize=8)

        # After: merged tiles in green
        draw_tiles(ax2, tiles_after, color_by_group=False, color_cycle=[merged_color])
        draw_polygon_outline(ax2, sf_verts, color=outline_color)
        ax2.set_title(f"After merge ({len(tiles_after)} tiles)", fontsize=10)
        multi_count = sum(
            1 for t in tiles_after if len(t.properties.get("spacecraft_ids", [])) > 1
        )
        ax2.plot(
            [],
            [],
            "s",
            color=merged_color,
            label=f"Merged (SC-1 + SC-2): {multi_count}",
            markersize=8,
        )
        ax2.legend(loc="lower right", fontsize=8)

        for ax in (ax1, ax2):
            style_gridlines(ax, theme)

        plt.tight_layout()
        save_themed(fig, "tessellation_merging", theme)

See Also