Granice państw offline: Shapely zamiast reverse geocodingu

Żeby rozliczyć diety zagraniczne kierowcy, BusiKM musi wiedzieć, ile godzin spędził w każdym państwie. Dane to punkty GPS z trasy, więc potrzebna jest jedna informacja: w którym państwie leży dana współrzędna. Zamiast pytać o to API do reverse geocodingu, trzymam granice państw w pliku w repozytorium i sprawdzam je lokalnie przez Shapely. W tym poście opisuję, dlaczego tak, jak działa wyszukiwanie w dwóch krokach i ile kosztuje uproszczona mapa.

Dlaczego nie Mapbox#

Mapbox ma reverse geocoding i BusiKM i tak korzysta z Mapboxa przy mapach. Problem jest w skali i w tym, czego naprawdę potrzebuję. Reverse geocoding to jedno zapytanie na punkt. Nocne przeliczenie dni za granicą dla całej floty to dziesiątki tysięcy punktów, czyli tyle samo zapytań, do tego limity API, koszt i wynik zależny od sieci w chwili przeliczania.

A z całej odpowiedzi API potrzebuję dwóch liter: kodu kraju. Nie ulicy, nie miasta, nie kodu pocztowego. Do tego wystarczy wiedzieć, wewnątrz którego wielokąta leży punkt, a granice państw zmieniają się rzadko. To zadanie, które da się policzyć lokalnie, bez sieci i za każdym razem tak samo.

Dane: Natural Earth przycięte do kodu i kształtu#

Granice pochodzą z Natural Earth w skali 1:50m. To dane w domenie publicznej, więc mogą leżeć w repozytorium bez problemów z licencją. Oryginalny plik ma około 90 kolumn atrybutów na państwo: nazwy w kilkunastu językach, statystyki, kody pomocnicze. Osobny skrypt pobiera go i zostawia tylko kod ISO-2 i geometrię. Plik ma wtedy 2,2 MB zamiast 3,1 MB i zawiera 239 państw.

Trzy terytoria nie mają kodu ISO-2: Cypr Północny, lodowiec Siachen i Somaliland. Skrypt je pomija, bo i tak nie dałoby się ich zapisać w polu z kodem kraju. Skrypt uruchamiam ręcznie, gdy trzeba odświeżyć dane. Nie jest krokiem wdrożenia ani CI, więc build nie zależy od tego, czy serwer z danymi akurat działa.

Wyszukiwanie w dwóch krokach#

Sprawdzanie punktu po kolei we wszystkich 239 wielokątach działałoby, ale robiłoby dużo niepotrzebnej pracy. Shapely ma STRtree, czyli drzewo R zbudowane z prostokątów otaczających każdą geometrię. Zapytanie o punkt zwraca tylko te państwa, których prostokąt go zawiera, i to jest szybkie.

Prostokąt to jednak nie granica. Punkt przy granicy polsko-niemieckiej leży w prostokątach obu państw, a punkt na Bałtyku w prostokącie Polski. Dlatego drugi krok sprawdza dokładny kształt, ale już tylko dla kandydatów:

class CountryBorders:
    def __init__(self, path: Path = BORDERS_PATH) -> None:
        payload = json.loads(path.read_text(encoding="utf-8"))
        self._codes = [f["properties"]["iso_a2"] for f in payload["features"]]
        self._geometries = [shape(f["geometry"]) for f in payload["features"]]
        self._tree = STRtree(self._geometries)
 
    def country_for(self, latitude: float, longitude: float) -> str | None:
        """ISO-2 code of the country the point lies in. None at sea."""
        point = Point(longitude, latitude)
        # The R-tree returns candidates by bounding box only;
        # covers() decides whether the point is really inside the border.
        for index in self._tree.query(point):
            if self._geometries[index].covers(point):
                return self._codes[index]
        return None
 
 
@lru_cache(maxsize=1)
def borders() -> CountryBorders:
    """Border index, built once per process."""
    return CountryBorders()

Bez covers() punkt na Bałtyku trafiałby do pierwszego państwa, którego prostokąt go obejmuje. Dla diet to nie jest drobiazg: godziny na morzu zapisane jako pobyt w danym kraju zawyżyłyby należność. Punkt poza lądem zwraca None i nie jest nigdzie przypisywany.

Drugi szczegół to lru_cache. Wczytanie 2,2 MB GeoJSON-a i zbudowanie drzewa trwa ułamek sekundy, ale robione dla każdego punktu byłoby nie do przyjęcia. Indeks powstaje raz na proces i potem obsługuje wszystkie zapytania. Jedna uwaga przy okazji: Point przyjmuje najpierw długość geograficzną, potem szerokość, czyli odwrotnie niż zwykle zapisuje się współrzędne GPS. Funkcja przyjmuje je w kolejności latitude, longitude, a zamienia dopiero przy tworzeniu punktu, żeby ten błąd nie przeciekał do reszty kodu.

Test, który pilnuje decyzji#

Wybór offline ma sens tylko wtedy, gdy nikt kiedyś nie dopisze tam zapytania do sieci. Dlatego jeden z testów podmienia urllib.request.urlopen na wyjątek i sprawdza, że Berlin dalej wychodzi jako DE:

def test_detection_needs_no_network():
    """The point of the decision: offline, no API limits, deterministic in CI."""
    with mock.patch("urllib.request.urlopen", side_effect=AssertionError("no network")):
        assert country_for_point(*BERLIN) == "DE"

Pozostałe testy sprawdzają miasta w kilku państwach, punkt na Bałtyku i trasy przez dwa i trzy państwa tego samego dnia. Ponieważ dane leżą w repozytorium, CI liczy dokładnie to samo co środowisko lokalne.

Z punktów na godziny#

Kod kraju dla punktu to dopiero połowa pracy. Punkty GPS są próbkami: mówią, gdzie kierowca był w danej chwili, a nie ile czasu spędził w którym państwie. Godziny powstają więc z odstępów między kolejnymi próbkami, a każdy odstęp należy do państwa wcześniejszego punktu, bo na jego początku kierowca był właśnie tam.

MAX_GAP = timedelta(minutes=30)
 
for current, following in zip(ordered, ordered[1:]):
    country = country_for_point(float(current["latitude"]), float(current["longitude"]))
    if country is None:
        continue  # at sea: no country to give the time to
    gap = following["timestamp"] - current["timestamp"]
    if gap <= timedelta(0):
        continue
    # A long break is a gap in the data, not a stay. Count it only up to the limit.
    counted_end = current["timestamp"] + min(gap, MAX_GAP)
    for day, part_seconds in _split_across_days(current["timestamp"], counted_end):
        seconds[(day, country)] += part_seconds

Najważniejsza decyzja w tym fragmencie to to, czego kod nie robi. Jeśli między próbkami minęły dwie godziny, kierowca mógł w tym czasie przejechać przez dwa państwa. Wymyślenie, które to były, byłoby nieprawdą zapisaną w rozliczeniu płacowym. Dlatego przerwa dłuższa niż 30 minut liczy się tylko do tego progu, reszta zostaje nieprzypisana, a o luce decyduje człowiek. Wpisy poprawione ręcznie są chronione: nocny przebieg ich nie nadpisuje, bo na nich opiera się wypłata.

Dzień liczę po dacie lokalnej, nie w UTC, bo dieta należy się za dzień kalendarzowy. Odstęp przez północ dzielę na dwie części, po jednej dla każdego dnia. O tym, jak łatwo pomylić datę lokalną z datą serwera, pisałem w poście o date.today() i datach faktur.

Ile kosztuje uproszczona mapa#

Skala 1:50m to uproszczone linie granic. Punkt kilkaset metrów od granicy może trafić do sąsiedniego państwa. Natural Earth ma też dokładniejszą skalę 1:10m, ale jej plik jest wielokrotnie większy, a wczytanie i budowa indeksu trwałyby dłużej. Dla diet liczy się godzina przekroczenia granicy, nie metr, a próbki GPS i tak padają co kilkadziesiąt sekund. Błąd przy samej granicy jest rzędu jednej próbki i na takie przypadki jest korekta ręczna w rejestrze dni za granicą.

Drugi koszt to aktualność danych. Granice nie zmieniają się często, ale skrypt odświeżający muszę uruchomić sam, gdy będzie trzeba. Nie jest częścią wdrożenia, więc nowe dane trafiają do repozytorium jak każda inna zmiana: przez PR i testy. Przy danych, od których zależą rozliczenia, to raczej zaleta niż koszt.

Ten wybór dobrze pokazuje, że pierwsza myśl przy mapach, czyli „użyj API”, nie zawsze jest najlepsza. Jeśli z całej odpowiedzi potrzebujesz jednej wartości, a dane zmieniają się rzadko, lokalny plik i dobry indeks są prostsze, tańsze i łatwiejsze do przetestowania. Jeśli masz podobny problem z danymi geograficznymi w swoim systemie, napisz do mnie.