add reverse geocoder based on natural earth datasets

pyproject.toml

@@ -24,10 +24,19 @@ requires-python = ">=3.9"
 dependencies = [
     "pyaro@git+https://github.com/metno/pyaro",
     "requests",
-    "geocoder",
+    "fiona",
+    "shapely",
+    "rtree",
     "tqdm",
 ]
 
+[tool.setuptools]
+packages = ["pyaro_readers", "geocoder_reverse_natural_earth"]
+package-dir = {"" = "src"}
+
+[tool.setuptools.package-data]
+geocoder_reverse_natural_earth = ["*.zip"]
+
 [project.urls]
 "Homepage" = "https://github.com/metno/pyaro-readers"
 "Bug Tracker" = "https://github.com/metno/pyaro-readers/issues"
@@ -54,6 +63,7 @@ skip_missing_interpreters = True
 isolated_build = True
 
 envlist =
+    py310
     py311
     format
 

src/geocoder_reverse_natural_earth/__init__.py

@@ -0,0 +1 @@
+from .geocoder_reverse_ne import Geocoder_Reverse_NE, Geocoder_Reverse_Exception

src/geocoder_reverse_natural_earth/geocoder_reverse_ne.py

@@ -0,0 +1,87 @@
+import os
+import fiona
+import rtree
+import shapely.geometry
+import shapely.speedups
+
+
+class Geocoder_Reverse_Exception(Exception):
+    pass
+
+
+class Geocoder_Reverse_NE:
+    _index = None
+
+    @staticmethod
+    def __generator_function():
+        ne = os.path.join(os.path.dirname(__file__), "ne_10m_admin_0_countries_deu.zip")
+        with fiona.open(f"zip://{ne}", "r") as fi:
+            for fid, feat in enumerate(fi):
+                geom = shapely.geometry.shape(feat["geometry"])
+                yield (fid, geom.bounds, (feat["properties"], geom))
+
+    """Get geocode from the Natural Earth Admin country shape-files"""
+
+    def __init__(self) -> None:
+        if Geocoder_Reverse_NE._index is None:
+            Geocoder_Reverse_NE._index = rtree.index.Index(
+                Geocoder_Reverse_NE.__generator_function()
+            )
+
+    def intersect(self, geometry) -> [dict()]:
+        """Find intersecting area with geometry
+
+        :param geometry: a shapely geometry
+        :return: a list of dicts of natural-earth properties
+        """
+        ret_list = []
+        for obj in self._index.intersection(geometry.bounds, objects=True):
+            props, geom = obj.object
+            if geometry.intersects(geom):
+                ret_list.append(props)
+        return ret_list
+
+    def lookup(self, latitude: float, longitude: float) -> dict():
+        """Lookup coordinates for a country
+
+        :param self: _description_
+        :raises Geocoder_Reverse_Exception: _description_
+        :return: _description_
+        """
+        point = shapely.geometry.Point((longitude, latitude))
+        objs = self.intersect(point)
+        if len(objs) > 0:
+            return objs[0]
+        else:
+            raise Geocoder_Reverse_Exception(
+                f"coordinates (latitude, longitude) not found in NE-admin"
+            )
+
+    def lookup_nearest(self, latitude: float, longitude: float) -> dict():
+        point = shapely.geometry.Point((longitude, latitude))
+        objs = self.intersect(point)
+        if len(objs) == 0:
+            nearest = list(self._index.nearest(point.bounds, 1, objects=True))
+            props, geom = nearest[0].object
+            return props
+        return objs[0]
+
+
+if __name__ == "__main__":
+    from tqdm import tqdm
+
+    geo = Geocoder_Reverse_NE()
+    print(geo.lookup(60, 10)["ISO_A2_EH"])
+    print(geo.lookup(78, 15)["ISO_A2_EH"])
+    print(geo.lookup_nearest(78.2361926, 15.3692614)["ISO_A2_EH"])
+    count = 0
+    for j in tqdm(range(-60, 85)):
+        for i in range(-180, 180):
+            count += 1
+            try:
+                geo.lookup(i, j)
+            except Geocoder_Reverse_Exception as grex:
+                pass
+            if False:
+                if geo.lookup_nearest(i, j) is None:
+                    print(f"error: {i},{j}")

src/pyaro_readers/aeronetsunreader/AeronetSunTimeseriesReader.py

@@ -4,7 +4,7 @@
 from urllib.request import urlopen
 from zipfile import BadZipFile, ZipFile
 
-import geocoder
+from geocoder_reverse_natural_earth import Geocoder_Reverse_NE
 import numpy as np
 import requests
 from pyaro.timeseries import (
@@ -107,6 +107,7 @@ def __init__(
         # get fields from header line although csv can do that as well since we might want to adjust these names
         self._fields = lines.pop(0).strip().split(",")
 
+        gcd = Geocoder_Reverse_NE()
         crd = csv.DictReader(lines, fieldnames=self._fields, **csvreader_kwargs)
         bar = tqdm(desc=tqdm_desc, total=len(lines))
         for _ridx, row in enumerate(crd):
@@ -120,10 +121,7 @@ def __init__(
                 alt = float(row["Site_Elevation(m)"])
                 if fill_country_flag:
                     try:
-                        country = geocoder.osm([lat, lon], method="reverse").json[
-                            "country_code"
-                        ]
-                        country = country.upper()
+                        country = gcd.lookup(lat, lon)["ISO_A2_EH"]
                     except:
                         country = "NN"
                 else:

tests/testdata/test_geocoder_reverse_natural_earth.py

@@ -0,0 +1,17 @@
+import unittest
+from geocoder_reverse_natural_earth import Geocoder_Reverse_NE, Geocoder_Reverse_Exception
+
+
+class TestGeocoderRevereseNE(unittest.TestCase):
+    def test_lookup(self):
+        geo = Geocoder_Reverse_NE()
+        self.assertEqual(geo.lookup(60,10)["ISO_A2_EH"], "NO")
+        self.assertEqual(geo.lookup(78, 15)["ISO_A2_EH"], "NO")
+        with self.assertRaises(Geocoder_Reverse_Exception):
+            geo.lookup(78.2361926, 15.3692614)
+
+    def test_lookup(self):
+        geo = Geocoder_Reverse_NE()
+        self.assertEqual(geo.lookup_nearest(78.2361926, 15.3692614)["ISO_A2_EH"], "NO")
+
+