Gdzie na świecie najbliżej do największej liczby państw?
Zastanawiałem się czy da się znaleźć punkt na Ziemi, z którego do kilku różnych państw jest w linii prostej naprawdę blisko? Aby to sprawdzić, przetestowałem algorytm oparty na danych kartograficznych i analizie przestrzennej. Nie szukałem odległości do stolic czy geometrycznych środków krajów – kluczowa była rzeczywista odległość od terytorium państwowego.
Metodologia: Jak to zostało policzone?
Do obliczeń wykorzystałem oficjalny zestaw danych granic państwowych Natural Earth (Admin 0 – Sovereignty), który zawiera geometrię suwerennych państw odwzorowującą stan faktyczny (de facto) przy użyciu wektorów 1:10m.
Pomysł na algorytm wyszukiwania podrzuciłem LLMom ChatGTP i Gemini. Proces opiera się na hierarchicznym przeszukiwaniu siatki globalnej:
- Wstępna gruboziarnista siatka – podział globu na komórki i szybkie wstępne filtrowanie obszarów za pomocą indeksów przestrzennych.
- Automatyczne zawężanie obszaru – rekurencyjne zagęszczanie siatki w miejscach o najwyższym potencjale (największej liczbie nakładających się buforów odległościowych).
- Lokalna optymalizacja – precyzyjne wyznaczanie minimalnego promienia i współrzędnych geograficznych za pomocą transformacji układów współrzędnych i obliczeń geodezyjnych.
Do implementacji użyłem Pythona oraz bibliotek GeoPandas (przetwarzanie danych GIS), Shapely (operacje geometryczne) i pyproj (transformacje układów odniesienia).
Wyniki obliczeń:
4 państwa są w promieniu około 250 m
- Współrzędne:
−17.7923, 25.2584 - Mapy: [OpenStreetMap] | [Google Maps]
- Państwa: Namibia, Zambia, Botswana, Zimbabwe
- Lokalizacja: Okolice Kazunguli nad rzekami Chobe i Zambezi.
- Kontekst: To jeden z najbardziej unikalnych punktów geopolitycznych świata. Wbrew obiegowym opiniom nie jest to idealny czwórstyk (punkt, w którym zbiegają się cztery granice), lecz dwa bardzo bliskie trójstyki połączone krótkim, kilkusetmetrowym odcinkiem granicy wodnej między Botswaną a Zambią, nad którym przebiega nowoczesny most Kazungula Bridge.
5 państw jest w promieniu około 35 km
- Współrzędne:
47.0223, 10.0599 - Mapy: [OpenStreetMap] | [Google Maps]
- Państwa: Austria, Szwajcaria, Niemcy, Włochy, Liechtenstein
- Lokalizacja: Alpy Centralne.
- Kontekst: Krajobraz drastically się zmienia – zamiast nizinnych rzek mamy wysokogórski teren. Odległości liczone są „z lotu ptaka”, co w praktyce oznacza, że 35 km w linii prostej na mapie przekłada się na znacznie dłuższą i trudniejszą trasę w rzeczywistości ze względu na ukształtowanie terenu.
6 państw jest w promieniu około 93 km
- Współrzędne:
47.1333, 8.3496 - Mapy: [OpenStreetMap] | [Google Maps]
- Państwa: Szwajcaria, Niemcy, Włochy, Francja, Liechtenstein, Austria
- Lokalizacja: Region alpejski w Europie Środkowej.
- Kontekst: Promień 93 km pozwala objąć zasięgiem terytoria aż sześciu państw, co wynika z zagęszczenia małych krajów i federacji w tej części Europy.
7 państw jest w promieniu około 139 km
- Współrzędne:
43.4801, 21.1570 - Mapy: [OpenStreetMap] | [Google Maps]
- Państwa: Serbia, Czarnogóra, Bułgaria, Bośnia i Hercegowina, Albania, Rumunia, Macedonia Północna
- Lokalizacja: Półwysep Bałkański.
- Kontekst: Największa liczba państw uchwycona w jednym algorithmicznym zapytaniu. Zbiega się tu skomplikowana historia podziałów terytorialnych regionu bałkańskiego.
Ograniczenia: Granica na mapie a rzeczywistość
Przedstawione wyniki należy traktować jako eksperyment matematyczno-geograficzny. Program operuje na czystej geometrii i mierzy najkrótszą odległość do poligonu państwowego, co wiąże się z kilkoma istotnymi zastrzeżeniami:
- Geografia terenu: W Kazunguli granice przebiegają przez koryta rzek Chobe i Zambezi. W Alpach promień 35 km w linii prostej oznacza konieczność pokonania masywów górskich, dolin i przełęczy.
- Infrastruktura i prawo: Model nie uwzględnia przejść granicznych, procedur celnych, wiz, stref zdemilitaryzowanych ani faktycznej dostępności terenu. „Państwo oddalone o 35 km” w sensie geometrii komputerowej nie oznacza, że można je legalnie i bez przeszkód osiągnąć w linii prostej.
- Dokładność danych: Granice w plikach wektorowych Natural Earth są opracowaniami kartograficznymi o określonej rozdzielczości. W miejscach spornych lub o skomplikowanym przebiegu linia demarkacyjna może się nieznacznie różnić od oficjalnych traktatów państwowych.
- Polityka: Inne dane można mieć jeśli uwzględni się Kosovo i inne miejsca o nieuregulowanym statusie politycznym.
Uzyskane wartości reprezentują najmniejsze promienie znalezione w toku zoptymalizowanego przeszukiwania przestrzennego, a nie bezwzględny dowód matematyczny. Zwiększenie rozdzielczości siatki obliczeniowej mogłoby w nielicznych przypadkach nieznacznie skorygować ostateczne odległości.
Kod:
#!/usr/bin/env python3
# search-meeting-borders-v14.py
#
# Szybkie wyszukiwanie punktu, z którego w promieniu R km znajduje się
# maksymalna liczba suwerennych państw.
#
# v14:
# - interaktywne zapytanie o zakres promieni na starcie (jeśli brak argumentów CLI)
# - śledzenie minimalnego promienia dla każdej unikalnej liczby państw
# - generowanie podsumowania minimalnych promieni na końcu pliku raportu .md
# - dokładny seed w lokalnym AEQD zamiast odległości w EPSG:6933
# - seed jest przekazywany dalej do branch-and-bound, więc rekord nie znika
# przy braku komórek frontier
# - raport filtruje listę państw po rzeczywistym dystansie <= R
from __future__ import annotations
import argparse
import math
import pickle
import sys
from datetime import datetime
from pathlib import Path
import geopandas as gpd
import numpy as np
import requests
from pyproj import CRS, Transformer
from shapely.geometry import box, Point
from shapely.ops import transform
# ---------------------------------------------------------------------------
# KONFIGURACJA
# ---------------------------------------------------------------------------
URL = (
"https://naturalearth.s3.amazonaws.com/10m_cultural/"
"ne_10m_admin_0_sovereignty.zip"
)
DATA_DIR = Path(__file__).resolve().parent
ZIP_PATH = DATA_DIR / "ne_10m_admin_0_sovereignty.zip"
CACHE_PATH = DATA_DIR / "country-radius-v14-cache.pkl"
GRID_CRS = CRS.from_epsg(6933)
START_CELL_KM = 1536.0
MAX_SUBDIVISIONS = 7
FINAL_GRID_N = 11
LOCAL_REFINE_N = 17
SHIFTS = (
(0.00, 0.00),
(0.25, 0.25),
)
DROP_NAMES = {
"antarctica",
"brazilian island",
"cyprus no mans area",
"kosovo",
}
# ---------------------------------------------------------------------------
# DANE
# ---------------------------------------------------------------------------
def download_data():
if ZIP_PATH.exists():
return
print("Pobieram warstwę Natural Earth Sovereignty...")
r = requests.get(URL, timeout=120)
r.raise_for_status()
ZIP_PATH.write_bytes(r.content)
print(f"Zapisano: {ZIP_PATH}")
def normal_name(x):
return " ".join(str(x).strip().lower().split())
def load_countries():
download_data()
print("Wczytuję granice państw...")
gdf = gpd.read_file(f"zip://{ZIP_PATH}")
field = "SOVEREIGNT"
if field not in gdf.columns:
raise RuntimeError(
f"Brak pola {field}. Dostępne pola: {list(gdf.columns)}"
)
print(f"Obiektów w warstwie Sovereignty: {len(gdf)}")
gdf = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty].copy()
gdf["_name_norm"] = gdf[field].map(normal_name)
gdf = gdf[~gdf["_name_norm"].isin(DROP_NAMES)].copy()
gdf = (
gdf[[field, "geometry"]]
.dissolve(by=field, as_index=False)
.rename(columns={field: "country"})
)
gdf["country_norm"] = gdf["country"].map(normal_name)
gdf = gdf[~gdf["country_norm"].isin(DROP_NAMES)].copy()
gdf = gdf.sort_values("country").reset_index(drop=True)
gdf["cid"] = np.arange(len(gdf), dtype=np.int32)
print(f"Państw do liczenia: {len(gdf)}")
projected = gdf.to_crs(GRID_CRS)
geoms = projected.geometry.to_numpy()
return gdf, geoms
# ---------------------------------------------------------------------------
# CACHE SIATKI
# ---------------------------------------------------------------------------
def load_or_create_grid_cache():
if CACHE_PATH.exists():
try:
with CACHE_PATH.open("rb") as f:
cache = pickle.load(f)
if (
cache.get("version") == 14
and cache.get("start_cell_km") == START_CELL_KM
):
print(f"Używam cache siatki: {CACHE_PATH}")
return cache
except Exception as e:
print(f"Nie można użyć cache ({e}), tworzę od nowa.")
cache = {
"version": 14,
"start_cell_km": START_CELL_KM,
"shifts": SHIFTS,
"cells": {},
}
with CACHE_PATH.open("wb") as f:
pickle.dump(cache, f, protocol=pickle.HIGHEST_PROTOCOL)
print(f"Utworzono cache siatki: {CACHE_PATH}")
return cache
# ---------------------------------------------------------------------------
# GLOBALNY ZASIĘG EPSG:6933
# ---------------------------------------------------------------------------
XMIN = -17367530.45
XMAX = 17367530.45
YMIN = -7324184.56
YMAX = 7324184.56
WORLD_W = XMAX - XMIN
WORLD_H = YMAX - YMIN
def cell_key(level, shift_id, ix, iy):
return (int(level), int(shift_id), int(ix), int(iy))
def make_level_cells(level: int, shift_id: int):
shift_x, shift_y = SHIFTS[shift_id]
size = START_CELL_KM * 1000.0 / (2 ** level)
dx = size * shift_x
dy = size * shift_y
x0 = XMIN - dx
y0 = YMIN - dy
nx = int(math.ceil(WORLD_W / size)) + 2
ny = int(math.ceil(WORLD_H / size)) + 2
cells = []
for iy in range(ny):
ymin = y0 + iy * size
ymax = ymin + size
if ymax <= YMIN or ymin >= YMAX:
continue
for ix in range(nx):
xmin = x0 + ix * size
xmax = xmin + size
parts = []
if xmin < XMIN:
parts.append(box(xmin + WORLD_W, ymin, XMAX, ymax))
parts.append(box(XMIN, ymin, xmax, ymax))
elif xmax > XMAX:
parts.append(box(xmin, ymin, XMAX, ymax))
parts.append(box(XMIN, ymin, xmax - WORLD_W, ymax))
else:
parts.append(box(xmin, ymin, xmax, ymax))
clipped = []
for p in parts:
p = p.intersection(box(XMIN, YMIN, XMAX, YMAX))
if not p.is_empty:
clipped.append(p)
if not clipped:
continue
cells.append((ix, iy, clipped))
return size, cells
def get_level_cells(level, shift_id, cache):
key = ("level", int(level), int(shift_id))
if key in cache["cells"]:
return cache["cells"][key]
size, cells = make_level_cells(level, shift_id)
serializable = [
(int(ix), int(iy), parts)
for ix, iy, parts in cells
]
cache["cells"][key] = {
"size_m": float(size),
"cells": serializable,
}
return cache["cells"][key]
# ---------------------------------------------------------------------------
# GEOMETRYCZNE SPRAWDZANIE KOMÓREK
# ---------------------------------------------------------------------------
def query_expanded_cell_country_ids(parts, radius_m, tree):
ids = set()
for p in parts:
q = p.buffer(radius_m)
for idx in tree.query(q, predicate="intersects"):
ids.add(int(idx))
q_left = transform(
lambda x, y, z=None: (x - WORLD_W, y), q
)
q_right = transform(
lambda x, y, z=None: (x + WORLD_W, y), q
)
for qq in (q_left, q_right):
for idx in tree.query(qq, predicate="intersects"):
ids.add(int(idx))
return tuple(sorted(ids))
# ---------------------------------------------------------------------------
# STARTOWY REKORD
# ---------------------------------------------------------------------------
def build_lonlat_seed_tree(gdf):
"""R-tree w EPSG:4326 używany wyłącznie do bezpiecznego wyboru kandydatów."""
return gpd.GeoSeries(gdf.geometry, crs="EPSG:4326").sindex
def exact_seed_candidate(lon, lat, radius_m, gdf, seed_tree):
"""Dokładnie policz państwa w promieniu od pojedynczego punktu.
Najpierw wybieramy kandydatów przez konserwatywny prostokąt geograficzny,
a następnie liczymy rzeczywistą odległość w lokalnym odwzorowaniu AEQD.
Dzięki temu seed jest prawdziwym dolnym ograniczeniem i nie może zawyżać
rekordu przez dystans mierzony w EPSG:6933.
"""
lat_delta = radius_m / 110000.0
min_lat = max(-90.0, lat - lat_delta)
max_lat = min(90.0, lat + lat_delta)
# Dla długości geograficznej uwzględniamy cały zakres szerokości
# prostokąta. Użycie najmniejszego cos(lat) sprawia, że filtr nie może
# odrzucić kraju tylko dlatego, że na brzegu prostokąta południk jest
# "tańszy" lub "droższy" niż w punkcie centralnym.
lat_abs_max = max(abs(min_lat), abs(max_lat))
cos_min = abs(math.cos(math.radians(min(89.999999, lat_abs_max))))
if lat_abs_max >= 89.999999 or cos_min < 1e-6:
lon_delta = 180.0
else:
lon_delta = min(180.0, radius_m / (110000.0 * cos_min))
min_lon = lon - lon_delta
max_lon = lon + lon_delta
query_boxes = []
if min_lon < -180.0:
query_boxes.append(box(min_lon + 360.0, min_lat, 180.0, max_lat))
query_boxes.append(box(-180.0, min_lat, max_lon, max_lat))
elif max_lon > 180.0:
query_boxes.append(box(min_lon, min_lat, 180.0, max_lat))
query_boxes.append(box(-180.0, min_lat, max_lon - 360.0, max_lat))
else:
query_boxes.append(box(min_lon, min_lat, max_lon, max_lat))
candidate_ids = set()
for q in query_boxes:
candidate_ids.update(int(i) for i in seed_tree.query(q, predicate="intersects"))
if not candidate_ids:
return ()
aeqd = CRS.from_proj4(
f"+proj=aeqd +lat_0={lat} +lon_0={lon} "
"+datum=WGS84 +units=m +no_defs"
)
to_local = Transformer.from_crs(4326, aeqd, always_xy=True).transform
p_local = Point(0.0, 0.0)
ids = []
for cid in candidate_ids:
local_geom = transform(to_local, gdf.geometry.iloc[cid])
if local_geom.distance(p_local) <= radius_m:
ids.append(cid)
return tuple(sorted(ids))
def cheap_seed(radius_km, gdf, geoms, tree, seed_tree):
"""Znajdź dokładny startowy rekord dla danego promienia.
Seed testuje punkty reprezentatywne i centroidy wszystkich państw.
Ponieważ każdy kandydat jest liczony dokładnie w lokalnym AEQD, wynik
może być użyty jako bezpieczne dolne ograniczenie w branch-and-bound.
"""
radius_m = radius_km * 1000.0
to_ll = Transformer.from_crs(GRID_CRS, 4326, always_xy=True).transform
candidates = []
seen = set()
for geom in geoms:
for p in (geom.representative_point(), geom.centroid):
lon, lat = to_ll(p.x, p.y)
key = (round(lon, 8), round(lat, 8))
if key not in seen:
seen.add(key)
candidates.append((lon, lat))
best = 0
best_point = candidates[0]
best_ids = ()
for lon, lat in candidates:
ids = exact_seed_candidate(lon, lat, radius_m, gdf, seed_tree)
if len(ids) > best:
best = len(ids)
best_point = (lon, lat)
best_ids = ids
return best, best_point, best_ids
# ---------------------------------------------------------------------------
# REKURENCJA
# ---------------------------------------------------------------------------
def search_radius(radius_km, gdf, geoms, tree, cache, initial_best, seed_xy, seed_ids):
radius_m = radius_km * 1000.0
best_count = initial_best
best_xy = seed_xy
best_ids = tuple(seed_ids)
print("\n" + "=" * 70)
print(f"RADIUS: {radius_km:g} KM")
print("=" * 70)
print(f"Startowy rekord: {best_count}")
frontier = []
for shift_id in range(2):
data = get_level_cells(0, shift_id, cache)
cells = data["cells"]
promising = []
for ix, iy, parts in cells:
ids = query_expanded_cell_country_ids(parts, radius_m, tree)
if len(ids) > best_count:
promising.append((ix, iy, parts, ids))
print(
f"Siatka {shift_id + 1}: "
f"{len(cells)} początkowych komórek, "
f"{len(promising)} wymagających dalszego sprawdzania"
)
for item in promising:
frontier.append((shift_id, 0, *item))
for level in range(1, MAX_SUBDIVISIONS + 1):
if not frontier:
break
next_frontier = []
for shift_id, parent_level, p_ix, p_iy, parent_parts, parent_ids in frontier:
for parent in parent_parts:
minx, miny, maxx, maxy = parent.bounds
mx = (minx + maxx) / 2.0
my = (miny + maxy) / 2.0
children = (
box(minx, miny, mx, my),
box(mx, miny, maxx, my),
box(minx, my, mx, maxy),
box(mx, my, maxx, maxy),
)
for child in children:
ids = query_expanded_cell_country_ids(
[child], radius_m, tree
)
if len(ids) > best_count:
next_frontier.append(
(shift_id, level, p_ix, p_iy, [child], ids)
)
frontier = next_frontier
cell_km = START_CELL_KM / (2 ** level)
print(
f"Poziom {level}: ~{cell_km:.2f} km, "
f"komórek do dalszego sprawdzania: {len(frontier)}"
)
if not frontier:
return best_count, best_xy, best_ids
print(f"Końcowe obszary do dokładnego sprawdzenia: {len(frontier)}")
seen = set()
for shift_id, level, ix, iy, parts, possible_ids in frontier:
for part in parts:
key = tuple(round(v, 1) for v in part.bounds)
if key in seen:
continue
seen.add(key)
cx = (part.bounds[0] + part.bounds[2]) / 2.0
cy = (part.bounds[1] + part.bounds[3]) / 2.0
to_ll = Transformer.from_crs(
GRID_CRS, 4326, always_xy=True
).transform
lon0, lat0 = to_ll(cx, cy)
if not (-90 <= lat0 <= 90):
continue
aeqd = CRS.from_proj4(
f"+proj=aeqd +lat_0={lat0} +lon_0={lon0} "
"+datum=WGS84 +units=m +no_defs"
)
to_local = Transformer.from_crs(
4326, aeqd, always_xy=True
).transform
to_local_grid = Transformer.from_crs(
GRID_CRS, aeqd, always_xy=True
).transform
from_local = Transformer.from_crs(
aeqd, 4326, always_xy=True
).transform
local_geoms = []
local_names = []
for cid in possible_ids:
lg = transform(to_local, gdf.geometry.iloc[cid])
local_geoms.append(lg)
local_names.append(gdf.country.iloc[cid])
if not local_geoms:
continue
local_part = transform(to_local_grid, part)
search_area = local_part.buffer(radius_m)
minx, miny, maxx, maxy = search_area.bounds
xs = np.linspace(minx, maxx, FINAL_GRID_N)
ys = np.linspace(miny, maxy, FINAL_GRID_N)
local_best = -1
local_best_xy = None
def count_at(x, y):
p = Point(float(x), float(y))
n = 0
for g in local_geoms:
if g.distance(p) <= radius_m:
n += 1
return n
for x in xs:
for y in ys:
p = Point(float(x), float(y))
if not search_area.contains(p):
continue
n = count_at(x, y)
if n > local_best:
local_best = n
local_best_xy = (x, y)
for g in local_geoms:
rp = g.representative_point()
x, y = rp.x, rp.y
if search_area.distance(Point(x, y)) == 0:
n = count_at(x, y)
if n > local_best:
local_best = n
local_best_xy = (x, y)
if local_best_xy is not None:
bx, by = local_best_xy
span = max(maxx - minx, maxy - miny) / 8.0
for x in np.linspace(bx - span, bx + span, LOCAL_REFINE_N):
for y in np.linspace(by - span, by + span, LOCAL_REFINE_N):
p = Point(float(x), float(y))
if not search_area.contains(p):
continue
n = count_at(x, y)
if n > local_best:
local_best = n
local_best_xy = (x, y)
if local_best > best_count and local_best_xy is not None:
best_count = local_best
lon, lat = from_local(*local_best_xy)
best_xy = (lon, lat)
ids = []
p = Point(*local_best_xy)
for cid, g in zip(possible_ids, local_geoms):
if g.distance(p) <= radius_m:
ids.append(cid)
best_ids = tuple(ids)
print(
f" nowy rekord {best_count}: "
f"{lat:.6f}, {lon:.6f}"
)
return best_count, best_xy, best_ids
# ---------------------------------------------------------------------------
# WYNIK I RAPORT .md
# ---------------------------------------------------------------------------
def process_result_data(radius, count, xy, ids, gdf):
if xy is None:
return None
lon, lat = xy
aeqd = CRS.from_proj4(
f"+proj=aeqd +lat_0={lat} +lon_0={lon} "
"+datum=WGS84 +units=m +no_defs"
)
to_local = Transformer.from_crs(4326, aeqd, always_xy=True).transform
p = transform(to_local, Point(lon, lat))
distances = []
radius_km = float(radius)
for cid in ids:
geom = transform(to_local, gdf.geometry.iloc[cid])
d = geom.distance(p) / 1000.0
if d <= radius_km + 1e-9:
distances.append((d, gdf.country.iloc[cid]))
sorted_distances = sorted(distances)
return {
"radius": radius,
"count": len(sorted_distances),
"lat": lat,
"lon": lon,
"sorted_distances": sorted_distances
}
def save_iter_result(res, md_report_path):
radius = res["radius"]
count = res["count"]
lat = res["lat"]
lon = res["lon"]
sorted_distances = res["sorted_distances"]
print("-" * 70)
print(f"WYNIK {radius:g} KM: {count} państw")
print(f"Punkt: {lat:.8f}, {lon:.8f}")
for d, name in sorted_distances:
print(f" {name:<35} {d:8.3f} km")
wynik_tekst = f"""### Promień: {radius:g} km
- **Liczba państw:** {count}
- **Współrzędne punktu:** `{lat:.8f}, {lon:.8f}`
- **Zestawienie odległości:**
"""
for d, name in sorted_distances:
wynik_tekst += f" - {name} – {d:.3f} km\n"
wynik_tekst += "\n---\n"
with md_report_path.open("a", encoding="utf-8") as f:
f.write(wynik_tekst)
print("-" * 70)
def save_summary_report(min_radius_results, md_report_path):
summary_text = "\n\n## Podsumowanie: Minimalne promienie dla liczby państw\n\n"
for count in sorted(min_radius_results.keys()):
res = min_radius_results[count]
summary_text += f"### Najmniejszy testowany promień dla {count} państw: {res['radius']:g} km\n"
summary_text += f"- **Współrzędne punktu:** `{res['lat']:.8f}, {res['lon']:.8f}`\n"
summary_text += "- **Zestawienie odległości:**\n"
for d, name in res['sorted_distances']:
summary_text += f" - {name} – {d:.3f} km\n"
summary_text += "\n"
with md_report_path.open("a", encoding="utf-8") as f:
f.write(summary_text)
print(f"\nDodano podsumowanie minimalnych promieni do raportu: {md_report_path}")
# ---------------------------------------------------------------------------
# MAIN
# ---------------------------------------------------------------------------
def main():
parser = argparse.ArgumentParser()
parser.add_argument(
"--radius-range",
nargs=2,
type=int,
default=None,
help="Zakres promieni od do, np. --radius-range 10 20",
)
parser.add_argument(
"--radii",
nargs="+",
type=float,
default=None,
help="Lista konkretnych promieni, np. --radii 10 15 20",
)
parser.add_argument(
"--rebuild-cache",
action="store_true",
help="usuń cache definicji siatki i zbuduj ponownie",
)
args = parser.parse_args()
# Określanie listy promieni (argumenty CLI lub interaktywne zapytanie)
if args.radius_range:
start_r, end_r = args.radius_range
radii_list = list(range(start_r, end_r + 1))
elif args.radii:
radii_list = args.radii
else:
print("\n--- KONFIGURACJA ZAKRESU PROMIENI ---")
try:
wejscie = input(
"Podaj zakres promieni (np. '1 50') "
"lub konkretne wartości po spacji: "
).strip()
parts = wejscie.split()
if len(parts) == 2:
start_r, end_r = int(parts[0]), int(parts[1])
radii_list = list(range(start_r, end_r + 1))
elif len(parts) > 2:
radii_list = [float(x) for x in parts]
elif len(parts) == 1:
radii_list = [float(parts[0])]
else:
print("Brak poprawnego wejścia, domyślnie zakres 1..100.")
radii_list = list(range(1, 101))
except (KeyboardInterrupt, EOFError):
print("\nPrzerwano.")
sys.exit(0)
except Exception:
print("Błąd parsowania, domyślnie zakres 1..100.")
radii_list = list(range(1, 101))
# Sortujemy promienie rosnąco, aby poprawnie wychwycić minimalny promień dla danej liczby państw
radii_list = sorted(list(set(radii_list)))
timestamp = datetime.now().strftime("%Y-%m-%d-%H-%M")
report_filename = f"raport-wynikow-v14-{timestamp}.md"
md_report_path = DATA_DIR / report_filename
if args.rebuild_cache and CACHE_PATH.exists():
CACHE_PATH.unlink()
gdf, geoms = load_countries()
print("Buduję indeks przestrzenny państw...")
tree = gpd.GeoSeries(geoms, crs=GRID_CRS).sindex
seed_tree = build_lonlat_seed_tree(gdf)
cache = load_or_create_grid_cache()
min_radius_results = {}
for radius in radii_list:
initial_best, seed_xy, seed_ids = cheap_seed(
radius,
gdf,
geoms,
tree,
seed_tree,
)
print(f"Startowy seed dla {radius:g} km: {initial_best}")
count, xy, ids = search_radius(
radius,
gdf,
geoms,
tree,
cache,
initial_best,
seed_xy,
seed_ids,
)
res = process_result_data(radius, count, xy, ids, gdf)
if res is not None:
save_iter_result(res, md_report_path)
# Zapisujemy jako minimalny testowany promień dla danej liczby państw,
# o ile jeszcze go nie ma (ponieważ radii_list jest posortowane rosnąco,
# pierwsze natrafienie na dany 'count' oznacza najmniejszy testowany promień).
c_val = res["count"]
if c_val not in min_radius_results:
min_radius_results[c_val] = res
with CACHE_PATH.open("wb") as f:
pickle.dump(cache, f, protocol=pickle.HIGHEST_PROTOCOL)
# Generowanie podsumowania minimalnych promieni na końcu pliku raportu
if min_radius_results:
save_summary_report(min_radius_results, md_report_path)
print(f"\nPełny raport zapisany w pliku: {md_report_path}")
if __name__ == "__main__":
main()
Dodaj komentarz