Offline Country Borders: Shapely vs Reverse Geocoding

To settle a driver's foreign per diems, BusiKM has to know how many hours he spent in each country. The input is GPS points from the route, so I need one piece of information: which country a given coordinate falls in. Instead of asking a reverse geocoding API, I keep country borders in a file in the repository and check them locally with Shapely. In this post I describe why, how the two-step lookup works, and what the simplified map costs.

Why not Mapbox#

Mapbox offers reverse geocoding, and BusiKM already uses Mapbox for maps. The problem is scale and what I actually need. Reverse geocoding is one request per point. A nightly recalculation of days abroad for the whole fleet means tens of thousands of points, so just as many requests, plus API limits, cost, and a result that depends on the network at the moment of recalculation.

And from the whole API response I need two letters: the country code. Not the street, not the city, not the postal code. All that's required is knowing which polygon contains the point, and country borders rarely change. This is a task you can compute locally, without a network, and get the same result every time.

Data: Natural Earth trimmed to code and shape#

The borders come from Natural Earth at 1:50m scale. It's public domain data, so it can live in the repository without licensing issues. The original file has about 90 attribute columns per country: names in a dozen or so languages, statistics, helper codes. A separate script downloads it and keeps only the ISO-2 code and the geometry. The file then weighs 2.2 MB instead of 3.1 MB and contains 239 countries.

Three territories have no ISO-2 code: Northern Cyprus, the Siachen Glacier and Somaliland. The script skips them, since they couldn't be stored in a country code field anyway. I run the script manually when the data needs refreshing. It is not a deployment or CI step, so the build doesn't depend on whether the data server happens to be up.

A two-step lookup#

Checking a point against all 239 polygons one by one would work, but it would do a lot of unnecessary work. Shapely has STRtree, an R-tree built from the bounding boxes of each geometry. A point query returns only the countries whose bounding box contains it, and that is fast.

A bounding box is not a border, though. A point near the Polish-German border lies in the boxes of both countries, and a point in the Baltic Sea lies in Poland's box. That's why the second step checks the exact shape, but only for the candidates:

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()

Without covers(), a point in the Baltic would be assigned to the first country whose box contains it. For per diems that's not a minor detail: hours at sea recorded as time spent in a given country would inflate the entitlement. A point outside land returns None and isn't assigned anywhere.

The second detail is lru_cache. Loading 2.2 MB of GeoJSON and building the tree takes a fraction of a second, but doing it for every point would be unacceptable. The index is built once per process and then serves all queries. One side note: Point takes longitude first, then latitude, the reverse of how GPS coordinates are usually written. The function accepts them in the order latitude, longitude and swaps them only when creating the point, so this mistake doesn't leak into the rest of the code.

A test that guards the decision#

Going offline only makes sense if nobody later adds a network call there. So one of the tests replaces urllib.request.urlopen with an exception and checks that Berlin still comes out as 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"

The remaining tests cover cities in several countries, a point in the Baltic Sea, and routes through two and three countries on the same day. Because the data lives in the repository, CI computes exactly the same thing as the local environment.

From points to hours#

The country code for a point is only half the work. GPS points are samples: they say where the driver was at a given moment, not how much time he spent in which country. So hours come from the gaps between consecutive samples, and each gap belongs to the country of the earlier point, because at its start the driver was there.

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

The most important decision in this fragment is what the code doesn't do. If two hours passed between samples, the driver could have crossed two countries in that time. Guessing which ones would be a falsehood recorded in payroll settlement. So a gap longer than 30 minutes counts only up to that threshold, the rest stays unassigned, and a human decides about the gap. Manually corrected entries are protected: the nightly run doesn't overwrite them, because pay is based on them.

I count the day by local date, not UTC, because a per diem is due for a calendar day. A gap spanning midnight is split into two parts, one for each day. I wrote about how easy it is to confuse local date with server date in the post on date.today() and invoice dates.

What the simplified map costs#

The 1:50m scale means simplified border lines. A point a few hundred meters from the border may land in the neighboring country. Natural Earth also has a more detailed 1:10m scale, but its file is many times larger, and loading it and building the index would take longer. For per diems, what matters is the hour of border crossing, not the meter, and GPS samples arrive every few dozen seconds anyway. The error right at the border is on the order of one sample, and for such cases there is a manual correction in the register of days abroad.

The second cost is data freshness. Borders don't change often, but I have to run the refresh script myself when needed. It's not part of the deployment, so new data goes into the repository like any other change: through a PR and tests. For data that settlements depend on, that's more of an advantage than a cost.

This choice shows well that the first thought when working with maps, "use an API", isn't always the best. If you need one value from the whole response and the data rarely changes, a local file and a good index are simpler, cheaper and easier to test. If you have a similar problem with geographic data in your system, get in touch.