Feat/wheel proj

python/eckit/src/_eckit/_eckit_geo.pyx

@@ -10,6 +10,8 @@
 
 cimport eckit_geo
 from cython.operator cimport dereference
+from libcpp.memory cimport unique_ptr
+from libcpp.string cimport string
 from libcpp.utility cimport pair
 from libcpp.vector cimport vector
 
@@ -30,6 +32,22 @@ def cache_dir_purge() -> None:
     eckit_geo.LibEcKitGeo.purgeCacheDir()
 
 
+def ensure_proj_database(fallback_db, fallback_search_paths=()) -> bool:
+    """Ensure eckit's PROJ subsystem can find a usable proj.db.
+
+    Respects PROJ_DATA/PROJ_LIB and any database PROJ can already resolve on
+    its own; only if neither yields a database does it fall back to the given
+    proj.db + search paths. No-op (returns False) when eckit was built without
+    PROJ. See LibEcKitGeo::ensureProjDatabase for the full policy.
+    """
+    cdef vector[string] paths
+    cdef string db
+    for p in fallback_search_paths:
+        paths.push_back(p.encode("utf-8") if isinstance(p, str) else p)
+    db = fallback_db.encode("utf-8") if isinstance(fallback_db, str) else fallback_db
+    return eckit_geo.LibEcKitGeo.ensureProjDatabase(db, paths)
+
+
 cdef class Area:
     cdef const eckit_geo.Area* _area
 
@@ -72,6 +90,170 @@ cdef class Area:
         return self._area.type()
 
 
+cdef class BoundingBox:
+    cdef eckit_geo.BoundingBox* _bbox
+
+    def __dealloc__(self):
+        if self._bbox != NULL:
+            del self._bbox
+
+    def __cinit__(self, north=None, west=None, south=None, east=None):
+        self._bbox = NULL
+
+        cdef unique_ptr[eckit_geo.BoundingBox] bbox
+        bbox = eckit_geo.BoundingBox.make_from_area(
+            float(north), float(west), float(south), float(east),
+        )
+        self._bbox = bbox.release()
+
+    def intersects(self, other) -> bool:
+        cdef BoundingBox other_bbox
+        cdef bint intersects
+        if not isinstance(other, BoundingBox):
+            raise TypeError("other must be a BoundingBox")
+        other_bbox = <BoundingBox>other
+        intersects = eckit_geo.bbox_intersects(
+            dereference(self._bbox), dereference(other_bbox._bbox)
+        )
+        return intersects
+
+    def contains(self, other) -> bool:
+        cdef BoundingBox other_bbox
+        if not isinstance(other, BoundingBox):
+            raise TypeError("other must be a BoundingBox")
+        other_bbox = <BoundingBox>other
+        return self._bbox.contains(dereference(other_bbox._bbox))
+
+    def contains_point(self, lon, lat) -> bool:
+        return eckit_geo.bbox_contains_lonlat(
+            dereference(self._bbox), float(lon), float(lat)
+        )
+
+    def as_list(self) -> list:
+        return [self.north, self.west, self.south, self.east]
+
+    @property
+    def spec_str(self) -> str:
+        return self._bbox.spec_str()
+
+    @property
+    def spec(self) -> dict:
+        from yaml import safe_load
+        return safe_load(self.spec_str)
+
+    @property
+    def north(self) -> float:
+        return self._bbox.north()
+
+    @property
+    def west(self) -> float:
+        return self._bbox.west()
+
+    @property
+    def south(self) -> float:
+        return self._bbox.south()
+
+    @property
+    def east(self) -> float:
+        return self._bbox.east()
+
+    @property
+    def global_(self) -> bool:
+        return self._bbox.is_global()
+
+    @property
+    def periodic(self) -> bool:
+        return self._bbox.periodic()
+
+    @property
+    def empty(self) -> bool:
+        return self._bbox.empty()
+
+    @property
+    def area(self) -> float:
+        return self._bbox.area()
+
+    def __repr__(self) -> str:
+        return str(self.as_list())
+
+    __str__ = __repr__
+
+
+cdef class Figure:
+    cdef const eckit_geo.Figure* _figure
+
+    def __dealloc__(self):
+        if self._figure != NULL:
+            del self._figure
+
+    def __cinit__(self, spec = None, **kwargs):
+        self._figure = NULL
+        assert bool(spec) != bool(kwargs)
+
+        if kwargs or isinstance(spec, dict):
+            from yaml import dump
+            spec = dump(kwargs if kwargs else spec, default_flow_style=True).strip()
+
+        try:
+            assert isinstance(spec, str)
+            self._figure = eckit_geo.FigureFactory.make_from_string(spec)
+
+        except RuntimeError as e:
+            # opportunity to do something interesting
+            raise
+
+    def __eq__(self, other) -> bool:
+        if not isinstance(other, Figure):
+            return NotImplemented
+        return self.spec_str == other.spec_str
+
+    def area(self, bbox = None) -> float:
+        cdef BoundingBox _bbox
+        if bbox is None:
+            return self._figure.area()
+        if not isinstance(bbox, BoundingBox):
+            raise TypeError("bbox must be a BoundingBox")
+        _bbox = <BoundingBox>bbox
+        return self._figure.area(dereference(_bbox._bbox))
+
+    @property
+    def R(self) -> float:
+        return eckit_geo.figure_R(dereference(self._figure))
+
+    @property
+    def a(self) -> float:
+        return self._figure.a()
+
+    @property
+    def b(self) -> float:
+        return self._figure.b()
+
+    @property
+    def spec_str(self) -> str:
+        return self._figure.spec_str()
+
+    @property
+    def spec(self) -> dict:
+        from yaml import safe_load
+        return safe_load(self.spec_str)
+
+    @property
+    def proj_str(self) -> str:
+        return self._figure.proj_str()
+
+    @property
+    def spherical(self) -> bool:
+        return self._figure.spherical()
+
+    @property
+    def eccentricity(self) -> float:
+        return self._figure.eccentricity()
+
+    @property
+    def flattening(self) -> float:
+        return self._figure.flattening()
+
+
 cdef class Grid:
     cdef const eckit_geo.Grid* _grid
 

python/eckit/src/_eckit/eckit_geo.pxd

@@ -8,6 +8,8 @@
 # does it submit to any jurisdiction.
 
 
+from libcpp cimport bool
+from libcpp.memory cimport unique_ptr
 from libcpp.string cimport string
 from libcpp.utility cimport pair
 from libcpp.vector cimport vector
@@ -21,6 +23,9 @@ cdef extern from "eckit/geo/LibEcKitGeo.h" namespace "eckit":
         @staticmethod
         void purgeCacheDir()
 
+        @staticmethod
+        bool ensureProjDatabase(const string& fallback_db, const vector[string]& fallback_search_paths) except +
+
         string version()
         string gitsha1(unsigned int n)  # n=40 for full sha1
 
@@ -37,11 +42,42 @@ cdef extern from "eckit/geo/Area.h" namespace "eckit::geo":
 
 cdef extern from "eckit/geo/area/BoundingBox.h" namespace "eckit::geo::area":
     cdef cppclass BoundingBox(Area):
+        BoundingBox(double north, double west, double south, double east) except +
+
         double north() const
         double west() const
         double south() const
         double east() const
 
+        bool intersects(BoundingBox&) const
+        bool contains(const BoundingBox&) const
+        bool is_global "global"() const
+        bool periodic() const
+        bool empty() const
+        double area() const
+
+        @staticmethod
+        unique_ptr[BoundingBox] make_from_area(
+            double n, double w, double s, double e
+        ) except +
+
+
+cdef extern from "eckit/geo/Figure.h" namespace "eckit::geo":
+    cdef cppclass Figure:
+        double a() const
+        double b() const
+        double area() const
+        double area(const BoundingBox&) const
+        string spec_str() const
+        string proj_str() const
+        bool spherical() const
+        double eccentricity() const
+        double flattening() const
+
+    cdef cppclass FigureFactory:
+        @staticmethod
+        Figure* make_from_string(const string) except +
+
 
 cdef extern from "eckit/geo/Range.h" namespace "eckit::geo":
     cdef cppclass Range:
@@ -75,7 +111,10 @@ cdef extern from * namespace "eckit::geo::python":
     """
     #include <vector>
 
+    #include "eckit/geo/area/BoundingBox.h"
+    #include "eckit/geo/Figure.h"
     #include "eckit/geo/Grid.h"
+    #include "eckit/geo/PointLonLat.h"
 
     namespace eckit::geo::python {
 
@@ -86,9 +125,25 @@ cdef extern from * namespace "eckit::geo::python":
     inline v grid_lon_values(const Grid& grid) { return grid.lon().values(); }
     inline v grid_lat_values(const Grid& grid) { return grid.lat().values(); }
 
+    inline bool bbox_intersects(const area::BoundingBox& lhs, area::BoundingBox& rhs) {
+        return lhs.intersects(rhs);
+    }
+
+    inline bool bbox_contains_lonlat(
+        const area::BoundingBox& bbox, double lon, double lat
+    ) {
+        return bbox.contains(PointLonLat{lon, lat});
+    }
+
+    inline double figure_R(const Figure& figure) { return figure.R(); }
+
     }
     """
     vector[double] grid_x_values(const Grid& grid) except +
     vector[double] grid_y_values(const Grid& grid) except +
     vector[double] grid_lon_values(const Grid& grid) except +
     vector[double] grid_lat_values(const Grid& grid) except +
+
+    bint bbox_intersects(const BoundingBox& lhs, BoundingBox& rhs) except +
+    bint bbox_contains_lonlat(const BoundingBox& bbox, double lon, double lat) except +
+    double figure_R(const Figure& figure) except +

python/eckit/src/eckit/geo/__init__.py

@@ -7,12 +7,43 @@
 # granted to it by virtue of its status as an intergovernmental organisation nor
 # does it submit to any jurisdiction.
 
+from pathlib import Path
+
 import findlibs
 
 findlibs.load("eckit", "eckitlib")
 findlibs.load("eckit_maths", "eckitlib")
 findlibs.load("eckit_geo", "eckitlib")
 
+
+def _configure_eckit_proj() -> None:
+    """
+    Point eckit's PROJ to the proj.db in the eckitlib wheel, so projections work out of the box.
+
+    All the policy is in LibEcKitGeo::ensureProjDatabase. Here we only discover the candidate
+    bundled location (it knows where pip installed the eckitlib package).
+    """
+    try:
+        import eckitlib
+    except ImportError:
+        return
+    eckitlib_file = getattr(eckitlib, "__file__", None)
+    if not eckitlib_file:
+        return
+
+    proj_dir = Path(eckitlib_file).parent / "share" / "proj"
+    from eckit.geo._eckit_geo import ensure_proj_database
+
+    if not ensure_proj_database(str(proj_dir / "proj.db"), [str(proj_dir)]):
+        import warnings
+
+        warnings.warn(
+            f"Could not find bundled proj.db at {proj_dir / 'proj.db'}, projection support is limited. Please install PROJ and/or set the PROJ_DATA environment variable."
+        )
+
+
+_configure_eckit_proj()
+
 from eckit._certifi import configure_ca_bundle_from_certifi
 from eckit.geo._eckit_geo import *
 

python/eckit/tests/test_area.py → python/eckit/tests/geo/test_area.py

@@ -8,20 +8,8 @@
 # does it submit to any jurisdiction.
 
 
-import pytest
-
-SPECS = [
-    (dict(north=90, west=0, south=-90, east=360), dict(area=[90, 0, -90, 360])),
-    (
-        dict(north=10, west=360 * 99 + 1, south=0, east=360 * 42 + 10),
-        dict(area=[10, 1, 0, 10]),
-    ),
-]
-
-
-@pytest.mark.parametrize("_spec, _expected", SPECS)
-def test_area(_spec, _expected):
+def test_area():
     from eckit.geo import Area
 
-    area = Area(_spec)
-    assert area.spec == _expected
+    area = Area(dict(north=90, west=0, south=-90, east=360))
+    assert area.spec == dict(area=[90, 0, -90, 360])