Issue259 allow epsg code as string

openeo/rest/datacube.py

@@ -37,7 +37,7 @@
 from openeo.rest.service import Service
 from openeo.rest.udp import RESTUserDefinedProcess
 from openeo.rest.vectorcube import VectorCube
-from openeo.util import get_temporal_extent, dict_no_none, rfc3339, guess_format
+from openeo.util import get_temporal_extent, dict_no_none, rfc3339, guess_format, crs_to_epsg_code
 
 if typing.TYPE_CHECKING:
     # Imports for type checking only (circular import issue at runtime).
@@ -333,7 +333,8 @@ def filter_bbox(
                               " Use keyword arguments or tuple/list argument instead.")
                 west, east, north, south = args[:4]
                 if len(args) > 4:
-                    crs = args[4]
+                    # TODO #259 convert CRS strings to integer EPSG code
+                    crs = crs_to_epsg_code(args[4])
             elif len(args) == 1 and (isinstance(args[0], (list, tuple)) and len(args[0]) == 4
                                      or isinstance(args[0], (dict, shapely.geometry.base.BaseGeometry, Parameter))):
                 bbox = args[0]
@@ -833,9 +834,20 @@ def _get_geometry_argument(
             ))
         if crs:
             # TODO: don't warn when the crs is Lon-Lat like EPSG:4326?
-            warnings.warn("Geometry with non-Lon-Lat CRS {c!r} is only supported by specific back-ends.".format(c=crs))
+            warnings.warn(f"Geometry with non-Lon-Lat CRS {crs!r} is only supported by specific back-ends.")
             # TODO #204 alternative for non-standard CRS in GeoJSON object?
-            geometry["crs"] = {"type": "name", "properties": {"name": crs}}
+            # TODO: #259 Should we also normalize EPSG codes passed as a string here?
+            #   It seems this method already supports epsg codes of the format "EPSG:<int>";
+
+            epsg_code = crs_to_epsg_code(crs)
+            if epsg_code is not None:
+                # proj did recognize the CRS
+                crs_name = f"EPSG:{epsg_code}"
+            else:
+                # proj did not recognise this CRS
+                warnings.warn(f"non-Lon-Lat CRS {crs!r} is not known to the proj library and might not be supported.")
+                crs_name = crs
+            geometry["crs"] = {"type": "name", "properties": {"name": crs_name}}
         return geometry
 
     @openeo_process

openeo/util.py

@@ -15,6 +15,7 @@
 from urllib.parse import urljoin
 
 import requests
+import pyproj.crs
 import shapely.geometry.base
 from deprecated import deprecated
 
@@ -531,6 +532,8 @@ class BBoxDict(dict):
     def __init__(self, *, west: float, south: float, east: float, north: float, crs: Optional[str] = None):
         super().__init__(west=west, south=south, east=east, north=north)
         if crs is not None:
+            # TODO: #259, should we covert EPSG strings to int here with crs_to_epsg_code?
+            # self.update(crs=crs_to_epsg_code(crs))
             self.update(crs=crs)
 
     # TODO: provide west, south, east, north, crs as @properties? Read-only or read-write?
@@ -627,3 +630,74 @@ def get(self, fraction: float) -> str:
         width = self.width - len(self.left) - len(self.right)
         bar = self.bar * int(round(width * clip(fraction, min=0, max=1)))
         return f"{self.left}{bar:{self.fill}<{width}s}{self.right}"
+
+
+class EPSGCodeNotFound(Exception):
+    """Could not find matching EPSG code for this CRS"""
+
+    def __init__(self, crs: Union[str, int, None], *args: object) -> None:
+        super().__init__(*args)
+        self.crs = crs
+
+
+def crs_to_epsg_code(crs: Union[str, int, None]) -> Optional[int]:
+    """Convert a CRS string or int to an integer EPGS code, where CRS usually comes from user input.
+
+    Three cases:
+
+    - If it is already an integer we just keep it.
+    - If it is None it stays None, and empty strings become None as well.
+    - If it is a string we try to parse it with the pyproj library.
+        - Strings of the form "EPSG:<int>" will be converted to teh value <int>
+        - For any other strings formats, it will work if pyproj supports is,
+          otherwise it won't.
+
+    The result is **always** an EPSG code, so the CRS should be one that is
+    defined in EPSG. For any other definitions pyproj will only give you the
+    closest EPSG match and that result is possibly inaccurate.
+
+    For a list of CRS input formats that proj supports
+    see: https://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.from_user_input
+
+    :param crs:
+        Input from user for the Coordinate Reference System to convert to an
+        EPSG code.
+
+    :raises ValueError:
+        When the crs is a negative value (as int or as a str representing an int)
+    :raises TypeError:
+        When crs is none of the supported types: str, int, None
+
+    :return: An EPGS code if it could be found, otherwise None
+
+    """
+
+    if crs is None or crs == "":
+        return None
+
+    # TODO: decide: are more fine-grained checks more helpful than always raising EPSGCodeNotFound?
+    if isinstance(crs, int):
+        if crs <= 0:
+            raise ValueError("If crs is an integer then it must be > 0.")
+        return crs
+
+    if not isinstance(crs, str):
+        raise TypeError("The allowed type for the parameter 'crs' are: str, int and None")
+
+    try:
+        crs_int = int(crs)
+    except:
+        # Need to process it with pyproj, below.
+        pass
+    else:
+        if crs_int <= 0:
+            raise ValueError("If crs is a string that represents an integer then the number must be > 0.")
+        return crs_int
+
+    try:
+        converted_crs = pyproj.crs.CRS.from_user_input(crs)
+    except pyproj.exceptions.CRSError as exc:
+        logger.error("Could not convert CRS string to EPSG code: {crs=}, exception: {exc}", exc_info=True)
+        raise EPSGCodeNotFound(crs) from exc
+    else:
+        return converted_crs.to_epsg()

tests/rest/datacube/test_datacube100.py

@@ -61,6 +61,64 @@
 ]
 
 
+WKT2_FOR_EPSG23631 = """
+PROJCRS["WGS 84 / UTM zone 31N",
+    BASEGEOGCRS["WGS 84",
+        ENSEMBLE["World Geodetic System 1984 ensemble",
+            MEMBER["World Geodetic System 1984 (Transit)"],
+            MEMBER["World Geodetic System 1984 (G730)"],
+            MEMBER["World Geodetic System 1984 (G873)"],
+            MEMBER["World Geodetic System 1984 (G1150)"],
+            MEMBER["World Geodetic System 1984 (G1674)"],
+            MEMBER["World Geodetic System 1984 (G1762)"],
+            MEMBER["World Geodetic System 1984 (G2139)"],
+            ELLIPSOID["WGS 84",6378137,298.257223563,
+                LENGTHUNIT["metre",1]],
+            ENSEMBLEACCURACY[2.0]],
+        PRIMEM["Greenwich",0,
+            ANGLEUNIT["degree",0.0174532925199433]],
+        ID["EPSG",4326]],
+    CONVERSION["UTM zone 31N",
+        METHOD["Transverse Mercator",
+            ID["EPSG",9807]],
+        PARAMETER["Latitude of natural origin",0,
+            ANGLEUNIT["degree",0.0174532925199433],
+            ID["EPSG",8801]],
+        PARAMETER["Longitude of natural origin",3,
+            ANGLEUNIT["degree",0.0174532925199433],
+            ID["EPSG",8802]],
+        PARAMETER["Scale factor at natural origin",0.9996,
+            SCALEUNIT["unity",1],
+            ID["EPSG",8805]],
+        PARAMETER["False easting",500000,
+            LENGTHUNIT["metre",1],
+            ID["EPSG",8806]],
+        PARAMETER["False northing",0,
+            LENGTHUNIT["metre",1],
+            ID["EPSG",8807]]],
+    CS[Cartesian,2],
+        AXIS["(E)",east,
+            ORDER[1],
+            LENGTHUNIT["metre",1]],
+        AXIS["(N)",north,
+            ORDER[2],
+            LENGTHUNIT["metre",1]],
+    USAGE[
+        SCOPE["Engineering survey, topographic mapping."],
+        AREA["Between 0°E and 6°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Andorra. Belgium. Benin. Burkina Faso. Denmark - North Sea. France. Germany - North Sea. Ghana. Luxembourg. Mali. Netherlands. Niger. Nigeria. Norway. Spain. Togo. United Kingdom (UK) - North Sea."],
+        BBOX[0,0,84,6]],
+    ID["EPSG",32631]]
+"""
+
+UTM31_CRS_STRINGS = [
+    "EPSG:32631",
+    "32631",
+    32631,
+    "+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs",  # is also EPSG:32631, in proj format
+    WKT2_FOR_EPSG23631,
+]
+
+
 def _get_leaf_node(cube: DataCube) -> dict:
     """Get leaf node (node with result=True), supporting old and new style of graph building."""
     flat_graph = cube.flat_graph()
@@ -276,12 +334,13 @@ def test_aggregate_spatial_types(con100: Connection, polygon, expected_geometrie
     }
 
 
-def test_aggregate_spatial_with_crs(con100: Connection, recwarn):
+@pytest.mark.parametrize("crs", UTM31_CRS_STRINGS)
+def test_aggregate_spatial_with_crs(con100: Connection, recwarn, crs: str):
     img = con100.load_collection("S2")
     polygon = shapely.geometry.box(0, 0, 1, 1)
-    masked = img.aggregate_spatial(geometries=polygon, reducer="mean", crs="EPSG:32631")
+    masked = img.aggregate_spatial(geometries=polygon, reducer="mean", crs=crs)
     warnings = [str(w.message) for w in recwarn]
-    assert "Geometry with non-Lon-Lat CRS 'EPSG:32631' is only supported by specific back-ends." in warnings
+    assert f"Geometry with non-Lon-Lat CRS {crs!r} is only supported by specific back-ends." in warnings
     assert sorted(masked.flat_graph().keys()) == ["aggregatespatial1", "loadcollection1"]
     assert masked.flat_graph()["aggregatespatial1"] == {
         "process_id": "aggregate_spatial",
@@ -300,6 +359,51 @@ def test_aggregate_spatial_with_crs(con100: Connection, recwarn):
     }
 
 
+@pytest.mark.parametrize(
+    "crs",
+    [
+        "WGS",  # WGS is not really specific enough, though WGS84 would have been fine
+        "does-not-exist-crs",
+        "EEPSG:32165",  # Simulate a user typo
+    ],
+)
+def test_aggregate_spatial_with_unusual_crs(con100: Connection, recwarn, crs: str):
+    img = con100.load_collection("S2")
+    polygon = shapely.geometry.box(0, 0, 1, 1)
+
+    import openeo.util
+
+    with pytest.raises(openeo.util.EPSGCodeNotFound):
+        masked = img.aggregate_spatial(geometries=polygon, reducer="mean", crs=crs)
+
+
+@pytest.mark.parametrize(
+    "crs",
+    [
+        "does-not-exist-crs",
+        "EEPSG:32165",  # Simulate a user typo
+        -1,  # integer non-sense, negative value can not be valid EPSG code
+        "-1",  # integer non-sense, negative value can not be valid EPSG code
+        1.0,  # floating point non-sense: type is not supported
+        "1.0",  # floating point non-sense: type is not supported
+        "0.0",  # floating point non-sense: type is not supported
+        {1: 1},  # type is not supported
+        [1],  # type is not supported
+    ],
+)
+def test_aggregate_spatial_with_invalid_crs(con100: Connection, recwarn, crs: str):
+    """Test that it refuses invalid input for the CRS: incorrect types and negative integers,
+    i.e. things that can not be a CRS at all, soo it is not just a CRS proj does not know about.
+    """
+    img = con100.load_collection("S2")
+    polygon = shapely.geometry.box(0, 0, 1, 1)
+
+    import openeo.util
+
+    with pytest.raises((ValueError, TypeError, openeo.util.EPSGCodeNotFound)):
+        img.aggregate_spatial(geometries=polygon, reducer="mean", crs=crs)
+
+
 def test_aggregate_spatial_target_dimension(con100: Connection):
     img = con100.load_collection("S2")
     polygon = shapely.geometry.box(0, 0, 1, 1)
@@ -476,19 +580,21 @@ def test_mask_polygon_types(con100: Connection, polygon, expected_mask):
     }
 
 
-def test_mask_polygon_with_crs(con100: Connection, recwarn):
+@pytest.mark.parametrize("crs", UTM31_CRS_STRINGS)
+def test_mask_polygon_with_crs(con100: Connection, recwarn, crs: str):
     img = con100.load_collection("S2")
     polygon = shapely.geometry.box(0, 0, 1, 1)
-    masked = img.mask_polygon(mask=polygon, srs="EPSG:32631")
+    masked = img.mask_polygon(mask=polygon, srs=crs)
     warnings = [str(w.message) for w in recwarn]
-    assert "Geometry with non-Lon-Lat CRS 'EPSG:32631' is only supported by specific back-ends." in warnings
+    assert f"Geometry with non-Lon-Lat CRS {crs!r} is only supported by specific back-ends." in warnings
     assert sorted(masked.flat_graph().keys()) == ["loadcollection1", "maskpolygon1"]
     assert masked.flat_graph()["maskpolygon1"] == {
         "process_id": "mask_polygon",
         "arguments": {
             "data": {"from_node": "loadcollection1"},
             "mask": {
                 "type": "Polygon", "coordinates": (((1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0), (1.0, 0.0)),),
+                # All listed test inputs for crs should be converted to "EPSG:32631"
                 "crs": {"type": "name", "properties": {"name": "EPSG:32631"}},
             },
         },