get_bearing_summary_stats produces inconsistent results with repeated aggregation

[dhimmel]

import ski_bearings.analyze

from ski_bearings.openskimap_utils import get_ski_area_to_runs, load_runs
from ski_bearings.osmnx_utils import create_networkx_with_metadata
from ski_bearings.analyze import aggregate_ski_areas_pl

ski_area_ids =  ['54bceeb62ccd2209e50b8dc7e06b8daa1727ce34', '74e0060a96e0399ace1b1e5ef5af1e5197a19752']

# create a pseudo ski area with all runs
runs = load_runs()
ski_area_to_runs = get_ski_area_to_runs(runs)
combined_runs = ski_area_to_runs[ski_area_ids[0]] + ski_area_to_runs[ski_area_ids[1]]
combined_graph = create_networkx_with_metadata(combined_runs, ski_area_metadata={})
combined_graph.graph
# mean_bearing_strength 0.7113152

ski_areas_pl = ski_bearings.analyze.load_ski_areas_pl()
ski_areas_pl.filter(pl.col("ski_area_id").is_in(ski_area_ids))

# aggregate the separate ski areas
aggregate_ski_areas_pl(group_by=["hemisphere"], ski_area_filters=[pl.col("ski_area_id").is_in(ski_area_ids)])
# mean_bearing_strength 0.9843988

[dhimmel]

Once this issue is fixed, a follow up question is whether we can take a weighted average of poleward/eastward affinities across ski areas.