diff --git a/python/eckit/src/_eckit/_eckit_geo.pyx b/python/eckit/src/_eckit/_eckit_geo.pyx index 7637a6f5d..4d6749f2b 100644 --- a/python/eckit/src/_eckit/_eckit_geo.pyx +++ b/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 = 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 = 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 = 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 diff --git a/python/eckit/src/_eckit/eckit_geo.pxd b/python/eckit/src/_eckit/eckit_geo.pxd index 42c57706c..0028148d8 100644 --- a/python/eckit/src/_eckit/eckit_geo.pxd +++ b/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 + #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 + diff --git a/python/eckit/src/eckit/geo/__init__.py b/python/eckit/src/eckit/geo/__init__.py index 067191e0d..16302cd14 100644 --- a/python/eckit/src/eckit/geo/__init__.py +++ b/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 * diff --git a/python/eckit/tests/test_area.py b/python/eckit/tests/geo/test_area.py similarity index 52% rename from python/eckit/tests/test_area.py rename to python/eckit/tests/geo/test_area.py index f71d00866..a5a6e37a0 100644 --- a/python/eckit/tests/test_area.py +++ b/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]) diff --git a/python/eckit/tests/geo/test_bbox.py b/python/eckit/tests/geo/test_bbox.py new file mode 100644 index 000000000..409431aa6 --- /dev/null +++ b/python/eckit/tests/geo/test_bbox.py @@ -0,0 +1,69 @@ +# (C) Copyright 1996- ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation nor +# does it submit to any jurisdiction. + + +import pytest + +from eckit.geo import BoundingBox + +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_bbox_spec(_spec, _expected): + bbox = BoundingBox(**_spec) + assert bbox.spec == _expected + + +def test_bbox_global(): + bbox = BoundingBox(north=90, west=0, south=-90, east=360) + + assert bbox.as_list() == [90, 0, -90, 360] + assert bbox.global_ + assert not bbox.empty + assert bbox.periodic + + +def test_bbox_not_global(): + bbox = BoundingBox(north=10, west=0, south=0, east=10) + + assert not bbox.global_ + assert not bbox.empty + assert not bbox.periodic + + +def test_bbox_intersects_and_contains(): + outer = BoundingBox(north=10, west=0, south=-10, east=10) + inner = BoundingBox(north=5, west=2, south=-5, east=8) + disjoint = BoundingBox(north=20, west=20, south=15, east=30) + + assert outer.intersects(inner) + assert outer.contains(inner) + assert not inner.contains(outer) + assert not outer.intersects(disjoint) + assert not outer.contains(disjoint) + + +def test_bbox_contains_point(): + bbox = BoundingBox(north=10, west=0, south=-10, east=10) + + assert bbox.contains_point(5, 5) + assert not bbox.contains_point(50, 50) + + +def test_bbox_area(): + bbox = BoundingBox(north=90, west=0, south=-90, east=360) + + assert bbox.area > 0 diff --git a/python/eckit/tests/geo/test_figure.py b/python/eckit/tests/geo/test_figure.py new file mode 100644 index 000000000..033e66615 --- /dev/null +++ b/python/eckit/tests/geo/test_figure.py @@ -0,0 +1,53 @@ +# (C) Copyright 1996- ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation nor +# does it submit to any jurisdiction. + + +import pytest + +from eckit.geo import BoundingBox +from eckit.geo import Figure + + +def test_figure_earth(): + R = 6371229.0 + A = 510101140207791.56 + + figure = Figure(R=R) + + assert figure.R == R + assert figure.a == figure.R + assert figure.b == figure.R + assert figure.spherical + assert figure.eccentricity == 0 + assert figure.flattening == 0 + assert figure.spec == dict(r=R) + + assert A == pytest.approx(figure.area()) + assert A == pytest.approx(figure.area(BoundingBox(north=90, west=0, south=-90, east=360))) + assert A / 2 == pytest.approx(figure.area(BoundingBox(north=90, west=0, south=0, east=360))) + + +def test_figure_wgs84(): + A = 510065621724079.0 + + figure = Figure(figure="wgs84") + + with pytest.raises(Exception): + figure.R # not defined for a != b + + assert figure.a == pytest.approx(6378137.0) + assert figure.b == pytest.approx(6356752.314245) + assert not figure.spherical + assert 0 < figure.eccentricity < 1 + assert figure.flattening == pytest.approx(1 - figure.b / figure.a) + assert figure.spec == dict(figure="wgs84") + + assert A == pytest.approx(figure.area()) + assert A == pytest.approx(figure.area(BoundingBox(north=90, west=0, south=-90, east=360))) + assert A / 2 == pytest.approx(figure.area(BoundingBox(north=90, west=0, south=0, east=360))) diff --git a/python/eckit/tests/test_grid.py b/python/eckit/tests/geo/test_grid.py similarity index 100% rename from python/eckit/tests/test_grid.py rename to python/eckit/tests/geo/test_grid.py diff --git a/python/eckitlib/buildconfig b/python/eckitlib/buildconfig index 5fa392b54..99ec516e2 100644 --- a/python/eckitlib/buildconfig +++ b/python/eckitlib/buildconfig @@ -13,6 +13,6 @@ NAME="eckit" _HERE="$(dirname -- "$(readlink -f -- "${BASH_SOURCE[0]}")")" PYTHON_INSTALL_DIR=$(cd -- "$_HERE/../eckit/src" && pwd) # dont use readlink -m here, not macos portable -CMAKE_PARAMS="-DENABLE_MPI=0 -DENABLE_ECKIT_GEO=1 -DENABLE_BUILD_TOOLS=OFF -DENABLE_AEC=0 -DENABLE_EIGEN=0 -DENABLE_LZ4=1 -DENABLE_PYTHON=1 -DPYTHONEXT_INSTALL_DIR=$PYTHON_INSTALL_DIR" +CMAKE_PARAMS="-DENABLE_MPI=0 -DENABLE_ECKIT_GEO=1 -DENABLE_BUILD_TOOLS=OFF -DENABLE_AEC=0 -DENABLE_EIGEN=0 -DENABLE_LZ4=1 -DENABLE_PYTHON=1 -DPYTHONEXT_INSTALL_DIR=$PYTHON_INSTALL_DIR -DENABLE_PROJ=1" PYPROJECT_DIR="python/eckitlib" DEPENDENCIES='[]' diff --git a/python/eckitlib/post-build.sh b/python/eckitlib/post-build.sh index c66cccca8..4e001fc16 100755 --- a/python/eckitlib/post-build.sh +++ b/python/eckitlib/post-build.sh @@ -2,18 +2,196 @@ set -euo pipefail +SCRIPT_DIR=$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")" && pwd) +ECKIT_SRC_DIR=$(cd -- "$SCRIPT_DIR/../eckit/src" && pwd) + +compute_relative_up() { + local rel_dir="$1" + if [ -z "$rel_dir" ] || [ "$rel_dir" = "." ]; then + echo "" + return + fi + local count + count=$(echo "$rel_dir" | awk -F'/' '{print NF}') + local up="" + for ((i=0; i/dev/null || true + + echo "running install_name_tool -add_rpath $rpath_eckit $bin_file" + install_name_tool -add_rpath "$rpath_eckit" "$bin_file" 2>/dev/null || true +} + +fix_another_wheel_linux() { + local vendor_dir="$1" + + echo "will crawl $ECKIT_SRC_DIR" + if [ -d "$ECKIT_SRC_DIR" ]; then + while read -r bin_file; do + [ -z "$bin_file" ] && continue + echo "found $bin_file, will process it" + process_binary_linux "$bin_file" "$ECKIT_SRC_DIR" "$vendor_dir" + done < <(find "$ECKIT_SRC_DIR" -type f -name "*.so") + fi +} + +fix_another_wheel_macos() { + local wheel_file="$1" + + echo "inspecting delocated wheel $wheel_file" + local inspect_dir="/tmp/eckit/delocate_inspect" + rm -rf "$inspect_dir" + mkdir -p "$inspect_dir" + unzip -q "$wheel_file" -d "$inspect_dir" + + local vendor_dir + vendor_dir=$(find "$inspect_dir" -maxdepth 2 -type d \( -name "*.dylibs" -o -name ".dylibs" -o -name "*.libs" \) | head -n 1) + + if [ -n "$vendor_dir" ] && [ -d "$vendor_dir" ]; then + echo "found vendor dylibs at $vendor_dir" + else + echo "warning: vendor dylibs directory not found in $wheel_file" + vendor_dir="$inspect_dir" + fi + + echo "will crawl $ECKIT_SRC_DIR" + if [ -d "$ECKIT_SRC_DIR" ]; then + while read -r bin_file; do + [ -z "$bin_file" ] && continue + echo "found $bin_file, will process it" + process_binary_macos "$bin_file" "$ECKIT_SRC_DIR" "$vendor_dir" + done < <(find "$ECKIT_SRC_DIR" -type f \( -name "*.so" -o -name "*.dylib" \)) + fi + + rm -rf "$inspect_dir" +} + # NOTE in case of problems like we had with eccodes, replace with noop here if [ "$(uname)" != "Darwin" ] ; then rm -rf /tmp/eckit/auditwheel auditwheel repair -w /tmp/eckit/auditwheel /tmp/eckit/build/wheel/*whl cd /tmp/eckit/auditwheel - F=$(ls *whl) - unzip $F + F=$(basename "$(ls ./*.whl)") + unzip "$F" + # shellcheck disable=SC2016 patchelf --add-rpath '$ORIGIN' eckitlib.libs/* - rm $F - zip -r $F ./* + + fix_another_wheel_linux "/tmp/eckit/auditwheel/eckitlib.libs" + + rm "$F" + zip -r "$F" ./* rm /tmp/eckit/build/wheel/* - mv /tmp/eckit/auditwheel/$F /tmp/eckit/build/wheel + mv "/tmp/eckit/auditwheel/$F" /tmp/eckit/build/wheel cd - fi @@ -21,4 +199,8 @@ fi # is disabled anyway if [ "$(uname)" = "Darwin" ] ; then delocate-wheel /tmp/eckit/build/wheel/*whl + WHEEL_FILE=$(find /tmp/eckit/build/wheel -name "*.whl" | head -n 1) + if [ -n "$WHEEL_FILE" ]; then + fix_another_wheel_macos "$WHEEL_FILE" + fi fi diff --git a/python/eckitlib/pre-compile.sh b/python/eckitlib/pre-compile.sh index 08bc49497..417494455 100755 --- a/python/eckitlib/pre-compile.sh +++ b/python/eckitlib/pre-compile.sh @@ -14,6 +14,11 @@ set -euo pipefail +# Source the buildconfig so we can inspect $CMAKE_PARAMS below (e.g. to decide +# whether to bundle PROJ's data files). +# shellcheck source=/dev/null +source python/eckitlib/buildconfig + mkdir -p python/eckitlib/src/copying mkdir -p /tmp/eckit/target/eckit/lib64/ @@ -24,12 +29,37 @@ if [ "$(uname)" != "Darwin" ] ; then # git clone https://github.com/lz4/lz4 /src/lz4 && cd /src/lz4 # make -j10 && make install DESTDIR=/tmp/lz4 # cd - + PROJ_ROOT="${PROJ_ROOT:-/cxx-deps}" else echo "no deps installation for platform $(uname)" + PROJ_ROOT="${PROJ_ROOT:-/tmp/cxx-deps}" fi wget https://raw.githubusercontent.com/lz4/lz4/dev/LICENSE -O python/eckitlib/src/copying/liblz4.txt -echo '{"liblz4": {"path": "copying/liblz4.txt", "home": "https://github.com/lz4/lz4"}}' > python/eckitlib/src/copying/list.json + +# Bundle PROJ's data files (proj.db, proj.ini, grids) inside the wheel if we +# built eckit with PROJ support. The runtime code in eckit/geo/__init__.py +# points eckit's bundled libproj at this dir via the PROJ context API, so it +# works offline without any system PROJ install and without leaking PROJ_DATA +# into the wider process (which would interfere with pyproj/fiona/GDAL/etc.). +if echo " $CMAKE_PARAMS " | grep -qE '[[:space:]]-DENABLE_PROJ=(1|ON)[[:space:]]' ; then + if [ ! -d "$PROJ_ROOT/share/proj" ] ; then + echo "ERROR: ENABLE_PROJ=1 but no PROJ data found at $PROJ_ROOT/share/proj" >&2 + echo " (override with PROJ_ROOT=/path/to/proj/prefix)" >&2 + exit 1 + fi + mkdir -p /tmp/eckit/target/eckit/share/proj + cp -r "$PROJ_ROOT/share/proj/." /tmp/eckit/target/eckit/share/proj/ + wget https://raw.githubusercontent.com/OSGeo/PROJ/master/COPYING -O python/eckitlib/src/copying/libproj.txt + cat > python/eckitlib/src/copying/list.json <<'JSON' +{ + "liblz4": {"path": "copying/liblz4.txt", "home": "https://github.com/lz4/lz4"}, + "libproj": {"path": "copying/libproj.txt", "home": "https://proj.org"} +} +JSON +else + echo '{"liblz4": {"path": "copying/liblz4.txt", "home": "https://github.com/lz4/lz4"}}' > python/eckitlib/src/copying/list.json +fi uv pip install cython diff --git a/share/eckit/geo/ORCA.yaml b/share/eckit/geo/ORCA.yaml index 4eef547d4..0c2699629 100644 --- a/share/eckit/geo/ORCA.yaml +++ b/share/eckit/geo/ORCA.yaml @@ -6,7 +6,7 @@ grid_names: name: ORCA2 intgrid: O42 arrangement: F - shape: [182, 149] + shape: [149, 182] uid: "174487fbace54b00d959d971e88b71e7" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA2_F.atlas @@ -15,7 +15,7 @@ grid_names: name: ORCA2 intgrid: O42 arrangement: T - shape: [182, 149] + shape: [149, 182] uid: "d5bde4f52ff3a9bea5629cd9ac514410" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA2_T.atlas @@ -24,7 +24,7 @@ grid_names: name: ORCA2 intgrid: O42 arrangement: U - shape: [182, 149] + shape: [149, 182] uid: "857f7affa3a381e3882d38d321384e49" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA2_U.atlas @@ -33,7 +33,7 @@ grid_names: name: ORCA2 intgrid: O42 arrangement: V - shape: [182, 149] + shape: [149, 182] uid: "ca637bc5dc9a54e2ea4b9750e1b79e6e" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA2_V.atlas @@ -42,7 +42,7 @@ grid_names: name: ORCA2 intgrid: O42 arrangement: W - shape: [182, 149] + shape: [149, 182] uid: "edea6f71eb558dc056b5f576d5b904f7" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA2_T.atlas @@ -51,7 +51,7 @@ grid_names: name: ORCA1 intgrid: O96 arrangement: F - shape: [362, 292] + shape: [292, 362] uid: "a832a12030c73928133553ec3a8d2a7e" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA1_F.atlas @@ -60,7 +60,7 @@ grid_names: name: ORCA1 intgrid: O96 arrangement: T - shape: [362, 292] + shape: [292, 362] uid: "f4c91b6233fe55dec992160ec12b38df" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA1_T.atlas @@ -69,7 +69,7 @@ grid_names: name: ORCA1 intgrid: O96 arrangement: U - shape: [362, 292] + shape: [292, 362] uid: "1b0f8d234753f910197c975c906b4da5" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA1_U.atlas @@ -78,7 +78,7 @@ grid_names: name: ORCA1 intgrid: O96 arrangement: V - shape: [362, 292] + shape: [292, 362] uid: "c637340454795b395f982851b840943d" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA1_V.atlas @@ -87,7 +87,7 @@ grid_names: name: ORCA1 intgrid: O96 arrangement: W - shape: [362, 292] + shape: [292, 362] uid: "d50061c43e83c46c3810002591ea21e1" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA1_T.atlas @@ -96,7 +96,7 @@ grid_names: name: eORCA1 intgrid: O96 arrangement: F - shape: [362, 332] + shape: [332, 362] uid: "3c6d95561710c6f39b394809ff6c588c" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_F.atlas @@ -105,7 +105,7 @@ grid_names: name: eORCA1 intgrid: O96 arrangement: T - shape: [362, 332] + shape: [332, 362] uid: "ba65665a9e68d1a8fa0352ecfcf8e496" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_T.atlas @@ -114,7 +114,7 @@ grid_names: name: eORCA1 intgrid: O96 arrangement: U - shape: [362, 332] + shape: [332, 362] uid: "4eb1054957dcae914e219faf9a4068e3" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_U.atlas @@ -123,7 +123,7 @@ grid_names: name: eORCA1 intgrid: O96 arrangement: V - shape: [362, 332] + shape: [332, 362] uid: "09131429766e7737c087d3a8d7073dc9" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_V.atlas @@ -132,7 +132,7 @@ grid_names: name: eORCA1 intgrid: O96 arrangement: W - shape: [362, 332] + shape: [332, 362] uid: "5c678d8f9aa2edfbf57246d11d9c1278" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_T.atlas @@ -141,7 +141,7 @@ grid_names: name: ORCA025 intgrid: O400 arrangement: F - shape: [1442, 1021] + shape: [1021, 1442] uid: "efbc280d8d4b6048797880da2605bacb" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA025_F.atlas @@ -150,7 +150,7 @@ grid_names: name: ORCA025 intgrid: O400 arrangement: T - shape: [1442, 1021] + shape: [1021, 1442] uid: "15c961c269ac182ca226d7195f3921ba" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA025_T.atlas @@ -159,7 +159,7 @@ grid_names: name: ORCA025 intgrid: O400 arrangement: U - shape: [1442, 1021] + shape: [1021, 1442] uid: "3f4a68bc5b54c9f867fbcc12aacc723d" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA025_U.atlas @@ -168,7 +168,7 @@ grid_names: name: ORCA025 intgrid: O400 arrangement: V - shape: [1442, 1021] + shape: [1021, 1442] uid: "9c87699ee2026c0feee07d2a972eaccd" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA025_V.atlas @@ -177,7 +177,7 @@ grid_names: name: ORCA025 intgrid: O400 arrangement: W - shape: [1442, 1021] + shape: [1021, 1442] uid: "74ca68f1c8524811f3d3aad99536adc2" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA025_T.atlas @@ -186,7 +186,7 @@ grid_names: name: eORCA025 intgrid: O400 arrangement: F - shape: [1442, 1207] + shape: [1207, 1442] uid: "770e5bbb667a253d55db8a98a3b2d3a9" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA025_F.atlas @@ -195,7 +195,7 @@ grid_names: name: eORCA025 intgrid: O400 arrangement: T - shape: [1442, 1207] + shape: [1207, 1442] uid: "983412216c9768bc794c18dc92082895" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA025_T.atlas @@ -204,7 +204,7 @@ grid_names: name: eORCA025 intgrid: O400 arrangement: U - shape: [1442, 1207] + shape: [1207, 1442] uid: "b1b2922e9b57ee9c6eeddad218b6e4f3" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA025_U.atlas @@ -213,7 +213,7 @@ grid_names: name: eORCA025 intgrid: O400 arrangement: V - shape: [1442, 1207] + shape: [1207, 1442] uid: "9b06bf73a8f14e927bd9b0f1f0c04f74" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA025_V.atlas @@ -222,7 +222,7 @@ grid_names: name: eORCA025 intgrid: O400 arrangement: W - shape: [1442, 1207] + shape: [1207, 1442] uid: "4a1ba3b11b8888aefc96992b6b1cab62" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA025_T.atlas @@ -231,7 +231,7 @@ grid_names: name: ORCA12 intgrid: O1600 arrangement: F - shape: [4322, 3059] + shape: [3059, 4322] uid: "29693ad8a7af3ae3ee0f02d090f0ec7b" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA12_F.atlas @@ -240,7 +240,7 @@ grid_names: name: ORCA12 intgrid: O1600 arrangement: T - shape: [4322, 3059] + shape: [3059, 4322] uid: "b117d01170ac77bca68560ab10e559de" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA12_T.atlas @@ -249,7 +249,7 @@ grid_names: name: ORCA12 intgrid: O1600 arrangement: U - shape: [4322, 3059] + shape: [3059, 4322] uid: "fff193b92d94d03e847ff2fa62b493f4" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA12_U.atlas @@ -258,7 +258,7 @@ grid_names: name: ORCA12 intgrid: O1600 arrangement: V - shape: [4322, 3059] + shape: [3059, 4322] uid: "986e3450774b716f6e75c1987e370b10" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA12_V.atlas @@ -267,7 +267,7 @@ grid_names: name: ORCA12 intgrid: O1600 arrangement: W - shape: [4322, 3059] + shape: [3059, 4322] uid: "ccfe953619a8dd49a7f765923882a274" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/ORCA12_T.atlas @@ -276,7 +276,7 @@ grid_names: name: eORCA12 intgrid: O1600 arrangement: F - shape: [4322, 3606] + shape: [3606, 4322] uid: "25da53ed581b3931fa310840fa9aefd9" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA12_F.atlas @@ -285,7 +285,7 @@ grid_names: name: eORCA12 intgrid: O1600 arrangement: T - shape: [4322, 3606] + shape: [3606, 4322] uid: "1553b66f5885cf5f83ad4b4fdf25f460" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA12_T.atlas @@ -294,7 +294,7 @@ grid_names: name: eORCA12 intgrid: O1600 arrangement: U - shape: [4322, 3606] + shape: [3606, 4322] uid: "3e87c826643da440b4e9d9f67588a576" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA12_U.atlas @@ -303,7 +303,7 @@ grid_names: name: eORCA12 intgrid: O1600 arrangement: V - shape: [4322, 3606] + shape: [3606, 4322] uid: "cc1e3fc06a2cd18c0653e557510b8a71" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA12_V.atlas @@ -312,7 +312,7 @@ grid_names: name: eORCA12 intgrid: O1600 arrangement: W - shape: [4322, 3606] + shape: [3606, 4322] uid: "462469edbd0e0586a0cf17424cc58c89" url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA12_T.atlas @@ -358,7 +358,7 @@ grid_uids: name: eORCA1 intgrid: O96 arrangement: T - shape: [362, 332] + shape: [332, 362] uid: "16076978a048410747dd7c9876677b28" # (uid older version) url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_T.atlas @@ -367,7 +367,7 @@ grid_uids: name: eORCA1 intgrid: O96 arrangement: U - shape: [362, 332] + shape: [332, 362] uid: "7378487847e050559b82d0792374a705" # (uid older version) url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_U.atlas @@ -376,6 +376,6 @@ grid_uids: name: eORCA1 intgrid: O96 arrangement: V - shape: [362, 332] + shape: [332, 362] uid: "d9622b55f3120eafb3dbaf5c742bc56c" # (uid older version) url: https://sites.ecmwf.int/repository/atlas/grids/orca/v0/eORCA1_V.atlas diff --git a/src/eckit/geo/CMakeLists.txt b/src/eckit/geo/CMakeLists.txt index 1982aa424..a0bc06060 100644 --- a/src/eckit/geo/CMakeLists.txt +++ b/src/eckit/geo/CMakeLists.txt @@ -174,6 +174,7 @@ set(eckit_geo_include_dirs ) set(eckit_geo_libs eckit_codec eckit_maths eckit_spec) +set(eckit_geo_private_libs) set(eckit_GEO_SHARE_AREA "~eckit/share/eckit/geo/area.yaml") set(eckit_GEO_SHARE_GRID "~eckit/share/eckit/geo/grid.yaml" @@ -193,7 +194,7 @@ if(eckit_HAVE_PROJ) projection/PROJ.cc projection/PROJ.h ) - list(APPEND eckit_geo_libs PROJ::proj) + list(APPEND eckit_geo_private_libs PROJ::proj) endif() if(eckit_HAVE_GEO_AREA_SHAPEFILE) @@ -203,7 +204,7 @@ if(eckit_HAVE_GEO_AREA_SHAPEFILE) cache/Unzip.cc cache/Unzip.h ) - list(APPEND eckit_geo_libs shapelib::shp libzip::zip) + list(APPEND eckit_geo_private_libs shapelib::shp libzip::zip) endif() string(REPLACE ";" ":" eckit_GEO_SHARE_GRID "${eckit_GEO_SHARE_GRID}") @@ -224,6 +225,7 @@ ecbuild_add_library( INSTALL_HEADERS ALL HEADER_DESTINATION ${INSTALL_INCLUDE_DIR}/eckit/geo PUBLIC_LIBS ${eckit_geo_libs} + PRIVATE_LIBS ${eckit_geo_private_libs} PUBLIC_INCLUDES ${eckit_geo_include_dirs} SOURCES ${eckit_geo_srcs} ) diff --git a/src/eckit/geo/Figure.cc b/src/eckit/geo/Figure.cc index 35659476e..ac7271c7b 100644 --- a/src/eckit/geo/Figure.cc +++ b/src/eckit/geo/Figure.cc @@ -15,6 +15,7 @@ #include #include "eckit/geo/Exceptions.h" +#include "eckit/geo/eckit_geo_config.h" #include "eckit/geo/figure/Earth.h" #include "eckit/geo/figure/OblateSpheroid.h" #include "eckit/geo/figure/Sphere.h" @@ -23,6 +24,10 @@ #include "eckit/spec/Custom.h" #include "eckit/types/FloatCompare.h" +#if eckit_HAVE_PROJ +#include "eckit/geo/projection/PROJ.h" +#endif + namespace eckit::geo { @@ -75,6 +80,16 @@ std::string Figure::spec_str() const { } +std::string Figure::proj_str() const { +#if eckit_HAVE_PROJ + std::unique_ptr custom(spec()); + return projection::PROJ::proj_str(*custom); +#else + NOTIMP; +#endif +} + + bool Figure::spherical() const { return types::is_approximately_equal(a(), b()); } diff --git a/src/eckit/geo/LibEcKitGeo.cc b/src/eckit/geo/LibEcKitGeo.cc index 6a4752817..911cf366f 100644 --- a/src/eckit/geo/LibEcKitGeo.cc +++ b/src/eckit/geo/LibEcKitGeo.cc @@ -21,6 +21,10 @@ #include "eckit/geo/eckit_geo_config.h" #include "eckit/utils/StringTools.h" +#if eckit_HAVE_PROJ +#include "eckit/geo/projection/PROJ.h" +#endif + namespace eckit { @@ -172,6 +176,18 @@ bool LibEcKitGeo::proj() { } +bool LibEcKitGeo::ensureProjDatabase(const std::string& fallback_db, + const std::vector& fallback_search_paths) { +#if eckit_HAVE_PROJ + return geo::projection::PROJ::proj_database_available(fallback_db, fallback_search_paths); +#else + static_cast(fallback_db); + static_cast(fallback_search_paths); + return false; +#endif +} + + const void* LibEcKitGeo::addr() const { return this; } diff --git a/src/eckit/geo/LibEcKitGeo.h b/src/eckit/geo/LibEcKitGeo.h index 1d5cc7fdc..798e7ecdd 100644 --- a/src/eckit/geo/LibEcKitGeo.h +++ b/src/eckit/geo/LibEcKitGeo.h @@ -12,6 +12,7 @@ #pragma once +#include #include #include "eckit/system/Library.h" @@ -50,6 +51,13 @@ class LibEcKitGeo final : public system::Library { static bool proj(); + /** + * If PROJ has a usable database. + * @return true if PROJ has a usable database, false otherwise + */ + static bool ensureProjDatabase(const std::string& fallback_db, + const std::vector& fallback_search_paths = {}); + std::string version() const override; std::string gitsha1(unsigned int count) const override; diff --git a/src/eckit/geo/area/library/Shapefile.cc b/src/eckit/geo/area/library/Shapefile.cc index d06f1b203..119e856ff 100644 --- a/src/eckit/geo/area/library/Shapefile.cc +++ b/src/eckit/geo/area/library/Shapefile.cc @@ -17,6 +17,8 @@ #include #include +#include "shapefil.h" + #include "eckit/geo/Exceptions.h" #include "eckit/geo/LibEcKitGeo.h" #include "eckit/geo/area/Polygon.h" @@ -31,6 +33,15 @@ namespace eckit::geo::area::library { +// PIMPL body: holds the shapelib-typed state so stays out of the +// installed public header. Kept in the .cc so that the unique_ptr in the +// Shapefile class only sees a complete Implementation type where its +// destructor is instantiated (i.e. in ~Shapefile() below). +struct Shapefile::Implementation { + SHPInfo* shp = nullptr; +}; + + namespace { @@ -91,16 +102,19 @@ Shapefile::Shapefile(const PathName& file) : Shapefile(file, "") {} Shapefile::Shapefile(const PathName& shp, const PathName& dbf, const std::string& name) : - shpPath_(path_shp(shp)), dbfPath_(path_dbf(dbf, shpPath_)), nEntities_(0) { + shpPath_(path_shp(shp)), + dbfPath_(path_dbf(dbf, shpPath_)), + impl_(std::make_unique()), + nEntities_(0) { Log::debug() << "eckit::geo::area::library::Shapefile(shp='" << shpPath_.realName() << "',dbf='" << dbfPath_.realName() << "',name='" << name << "')" << std::endl; - if ((shp_ = SHPOpen(shpPath_.localPath(), "rb")) == nullptr) { + if ((impl_->shp = SHPOpen(shpPath_.localPath(), "rb")) == nullptr) { throw CantOpenFile(shpPath_ + " (as .shp)", Here()); } int type = 0; - SHPGetInfo(shp_, &nEntities_, &type, nullptr, nullptr); + SHPGetInfo(impl_->shp, &nEntities_, &type, nullptr, nullptr); if (type != SHPT_ARC && type != SHPT_POLYGON) { throw ReadError("Shapefile: unsupported shape type", Here()); @@ -142,7 +156,7 @@ Shapefile::Shapefile(const PathName& shp, const PathName& dbf, const std::string Shapefile::~Shapefile() { - SHPClose(shp_); + SHPClose(impl_->shp); } @@ -200,7 +214,7 @@ Area* Shapefile::make_area(size_t entity) const { explicit Object(SHPObject* ptr) : unique_ptr{ptr, SHPDestroyObject} { ASSERT(operator bool()); } }; - Object obj(SHPReadObject(shp_, static_cast(entity))); + Object obj(SHPReadObject(impl_->shp, static_cast(entity))); ASSERT(obj); ASSERT(obj->nSHPType == SHPT_ARC || obj->nSHPType == SHPT_POLYGON); diff --git a/src/eckit/geo/area/library/Shapefile.h b/src/eckit/geo/area/library/Shapefile.h index 3922e0301..fd5bada78 100644 --- a/src/eckit/geo/area/library/Shapefile.h +++ b/src/eckit/geo/area/library/Shapefile.h @@ -14,13 +14,12 @@ #include #include +#include #include #include "eckit/filesystem/PathName.h" #include "eckit/geo/area/Library.h" -#include "shapefil.h" - namespace eckit::geo { class Area; @@ -65,12 +64,16 @@ class Shapefile : public Library { private: + // -- Types + + struct Implementation; + // -- Members const PathName shpPath_; const PathName dbfPath_; - SHPInfo* shp_; + std::unique_ptr impl_; std::string name_; int nEntities_; diff --git a/src/eckit/geo/figure/OblateSpheroid.cc b/src/eckit/geo/figure/OblateSpheroid.cc index 0958ab478..be8741738 100644 --- a/src/eckit/geo/figure/OblateSpheroid.cc +++ b/src/eckit/geo/figure/OblateSpheroid.cc @@ -86,7 +86,7 @@ double OblateSpheroid::_area(double a, double b, const area::BoundingBox& bbox) const auto dlam = util::DEGREE_TO_RADIAN * (bbox.east() - bbox.west()); const auto e = eccentricity(a, b); - const auto A = dlam * a * b * (f(phi2, e) - f(phi1, e)); + const auto A = dlam * b * b / 2. * (f(phi2, e) - f(phi1, e)); return A; } diff --git a/src/eckit/geo/projection/PROJ.cc b/src/eckit/geo/projection/PROJ.cc index 5d6ebc5eb..38297c5f8 100644 --- a/src/eckit/geo/projection/PROJ.cc +++ b/src/eckit/geo/projection/PROJ.cc @@ -14,10 +14,13 @@ #include +#include #include #include #include +#include +#include "eckit/filesystem/PathName.h" #include "eckit/geo/Exceptions.h" #include "eckit/geo/Figure.h" #include "eckit/spec/Custom.h" @@ -32,11 +35,23 @@ static ProjectionRegisterType PROJECTION("proj"); namespace { -constexpr auto CTX = PJ_DEFAULT_CTX; +constexpr auto CTX = PJ_DEFAULT_CTX; +constexpr PJ_AREA* DEFAULT_AREA = nullptr; struct pj_t : std::unique_ptr { - explicit pj_t(element_type* ptr) : unique_ptr(ptr, &proj_destroy) {} + explicit pj_t(element_type* ptr) : unique_ptr(ptr, &proj_destroy) { + if (!operator bool()) { + // common errors are "proj.db not found" or "invalid CRS string" + const auto err = proj_context_errno(CTX); + throw exception::ProjectionError( + "PROJ: failed to create object (err=" + std::to_string(err) + ", description='" + + proj_errno_string(err) + + "'). Ensure proj.db is available (https://proj.org/en/stable/resource_files.html) or install " + "the eckitlib wheel which bundles its own).", + Here()); + } + } }; @@ -93,7 +108,7 @@ struct XYZ final : Convert { Figure* make_figure(const std::string& proj_str) { - pj_t identity(proj_create_crs_to_crs(CTX, proj_str.c_str(), proj_str.c_str(), nullptr)); + pj_t identity(proj_create_crs_to_crs(CTX, proj_str.c_str(), proj_str.c_str(), DEFAULT_AREA)); pj_t crs(proj_get_target_crs(CTX, identity.get())); pj_t ellipsoid(proj_get_ellipsoid(CTX, crs.get())); @@ -142,7 +157,7 @@ PROJ::PROJ(const std::string& source, const std::string& target, double lon_mini ASSERT(!target_.empty()); auto make_convert = [lon_minimum](const std::string& string) -> Convert* { - pj_t identity(proj_create_crs_to_crs(CTX, string.c_str(), string.c_str(), nullptr)); + pj_t identity(proj_create_crs_to_crs(CTX, string.c_str(), string.c_str(), DEFAULT_AREA)); pj_t crs(proj_get_target_crs(CTX, identity.get())); pj_t cs(proj_crs_get_coordinate_system(CTX, crs.get())); ASSERT(cs); @@ -158,11 +173,10 @@ PROJ::PROJ(const std::string& source, const std::string& target, double lon_mini }; // projection, normalised - auto ctx = PJ_DEFAULT_CTX; + pj_t p(proj_create_crs_to_crs(CTX, source_.c_str(), target_.c_str(), DEFAULT_AREA)); + p.reset(proj_normalize_for_visualization(CTX, p.release())); - implementation_ = std::make_unique( - proj_normalize_for_visualization(ctx, proj_create_crs_to_crs(ctx, source_.c_str(), target_.c_str(), nullptr)), - ctx, make_convert(source_), make_convert(target_)); + implementation_ = std::make_unique(p.release(), CTX, make_convert(source_), make_convert(target_)); ASSERT(implementation_); } @@ -260,6 +274,67 @@ const std::string& PROJ::proj_default() { } +bool PROJ::proj_database_available(const std::string& fallback_db, + const std::vector& fallback_search_paths) { + // Resolution order: + // 1. If PROJ can resolve a database on its own (via its compiled-in default search paths, e.g. a system install) + // 2. If @p fallback_db points at a @c proj.db file, with @p fallback_search_paths used to locate @c proj.ini and + // grid files. + // + // The fallback is applied via the PROJ per-context API, so it affects only eckit's libproj instance and never leaks + // into other PROJ users in the process (pyproj, fiona, GDAL, ...). It is a no-op when eckit was built without PROJ + // support. Intended to be called once, at initialisation, before any PROJ-backed projection is created. + + struct Database { + Database(const std::string& fallback_db, const std::vector& fallback_search_paths) { + auto database_available = [&]() -> bool { + const auto previous_level = proj_log_level(CTX, PJ_LOG_NONE); + + auto avail = false; + try { + pj_t crs(proj_create_from_database(CTX, "EPSG", "4326", PJ_CATEGORY_CRS, false, nullptr)); + avail = static_cast(crs); + } + catch (...) { + } + + proj_log_level(CTX, previous_level); + return avail; + }; + + // (1) If PROJ already finds a usable database on its own, leave it alone. + if (database_available()) { + available = true; + return; + } + + // (2) Only now fall back to the provided (bundled) database, if present. + if (!fallback_db.empty() && PathName{fallback_db}.exists()) { + proj_context_set_database_path(CTX, fallback_db.c_str(), nullptr, nullptr); + + std::vector paths; + paths.reserve(fallback_search_paths.size()); + for (const auto& p : fallback_search_paths) { + if (!p.empty() && PathName{p}.exists()) { + paths.push_back(p.c_str()); + } + } + + if (!fallback_search_paths.empty()) { + proj_context_set_search_paths(CTX, static_cast(paths.size()), paths.data()); + } + + available = database_available(); + } + } + + bool available = false; + } static const DATABASE(fallback_db, fallback_search_paths); + + return DATABASE.available; +} + + void PROJ::fill_spec(spec::Custom& custom) const { custom.set("type", "proj"); if (source_ != proj_default()) { diff --git a/src/eckit/geo/projection/PROJ.h b/src/eckit/geo/projection/PROJ.h index 8ec1a7190..eff500f63 100644 --- a/src/eckit/geo/projection/PROJ.h +++ b/src/eckit/geo/projection/PROJ.h @@ -51,6 +51,13 @@ class PROJ : public Projection { static std::string proj_str(const spec::Custom&); static const std::string& proj_default(); + /** + * If PROJ has a usable database. + * @return true if PROJ has a usable database, false otherwise + */ + static bool proj_database_available(const std::string& fallback_db, + const std::vector& fallback_search_paths = {}); + private: // -- Types