Add CRS support for queries

README.md

@@ -69,7 +69,7 @@ In the future, when `xpublish` supports [`DataTree`](https://docs.xarray.dev/en/
 | `z`  | ✅ | |
 | `datetime`  | ✅ | |
 | `parameter-name`  | ✅   | |
-| `crs`  | ❌ | Not currently supported, all coordinates should be in the reference system of the queried dataset |
+| `crs`  | ✅  | Requires a CF compliant [grid mapping](https://cf-xarray.readthedocs.io/en/latest/grid_mappings.html) on the target dataset. Default is `EPSG:4326` |
 | `parameter-name`  | ✅ | |
 | `f`  | ✅ | |
 | `method`  | ➕ | Optional: controls data selection. Use "nearest" for nearest neighbor selection, or "linear" for interpolated selection. Uses `nearest` if not specified |
@@ -84,7 +84,7 @@ In the future, when `xpublish` supports [`DataTree`](https://docs.xarray.dev/en/
 | `z`  | ✅   | |
 | `datetime`  | ✅ | |
 | `parameter-name`  | ✅   | |
-| `crs`  | ❌ | Not currently supported, all coordinates should be in the reference system of the queried dataset |
+| `crs`  | ✅  | Requires a CF compliant [grid mapping](https://cf-xarray.readthedocs.io/en/latest/grid_mappings.html) on the target dataset. Default is `EPSG:4326` |
 | `parameter-name`  | ✅   | |
 | `f`  | ✅   | |
 | `method`  | ➕ | Optional: controls data selection. Use "nearest" for nearest neighbor selection, or "linear" for interpolated selection. Uses `nearest` if not specified |

tests/test_select.py

@@ -17,6 +17,14 @@
     return xr.tutorial.load_dataset("air_temperature")
 
 
+@pytest.fixture(scope="function")
+def projected_xy_dataset():
+    """Loads a sample dataset with projected X and Y coordinates"""
+    from cf_xarray.datasets import rotds
+
+    return rotds
+
+
 def test_select_query(regular_xy_dataset):
     query = EDRQuery(
         coords="POINT(200 45)",
@@ -127,13 +135,27 @@
     assert "lat" in ds, "Dataset does not contain the lat variable"
     assert "lon" in ds, "Dataset does not contain the lon variable"
 
     assert ds["air"].shape == ds["time"].shape, "Dataset shape is incorrect"
     npt.assert_array_equal(ds["lat"], 45.0), "Latitude is incorrect"
     npt.assert_array_equal(ds["lon"], 205.0), "Longitude is incorrect"
     npt.assert_approx_equal(ds["air"][0], 280.2), "Temperature is incorrect"
     npt.assert_approx_equal(ds["air"][-1], 279.19), "Temperature is incorrect"
 
 
+def test_select_position_projected_xy(projected_xy_dataset):
+    from xpublish_edr.geometry.common import project_geometry
+
+    point = Point((64.59063409, 66.66454929))
+    projected_point = project_geometry(projected_xy_dataset, "EPSG:4326", point)
+    npt.assert_approx_equal(projected_point.x, 18.045), "Longitude is incorrect"
+    npt.assert_approx_equal(projected_point.y, 21.725), "Latitude is incorrect"
+
+    ds = select_by_position(projected_xy_dataset, projected_point)
+    npt.assert_approx_equal(ds.rlon.values, 18.045), "Longitude is incorrect"
+    npt.assert_approx_equal(ds.rlat.values, 21.725), "Latitude is incorrect"
+    npt.assert_approx_equal(ds.temp.values, 0.89959461), "Temperature is incorrect"
+
+
 def test_select_position_regular_xy_interpolate(regular_xy_dataset):
     point = Point((204, 44))
     ds = select_by_position(regular_xy_dataset, point, method="linear")

xpublish_edr/geometry/common.py

@@ -2,10 +2,19 @@
 Common geometry handling functions
 """
 
+import itertools
+from functools import lru_cache
+
+import pyproj
 import xarray as xr
+from shapely import Geometry
+from shapely.ops import transform
 
 VECTORIZED_DIM = "pts"
 
+# https://pyproj4.github.io/pyproj/stable/advanced_examples.html#caching-pyproj-objectshttps://pyproj4.github.io/pyproj/stable/advanced_examples.html#caching-pyproj-objects
+transformer_from_crs = lru_cache(pyproj.Transformer.from_crs)
+
 
 def coord_is_regular(da: xr.DataArray) -> bool:
     """
@@ -19,3 +28,26 @@ def is_regular_xy_coords(ds: xr.Dataset) -> bool:
     Check if the dataset has 2D coordinates
     """
     return coord_is_regular(ds.cf["X"]) and coord_is_regular(ds.cf["Y"])
+
+
+def project_geometry(ds: xr.Dataset, geometry_crs: str, geometry: Geometry) -> Geometry:
+    """
+    Get the projection from the dataset
+    """
+    grid_mapping_names = ds.cf.grid_mapping_names
+    if len(grid_mapping_names) == 0:
+        # TODO: Should we require a grid mapping? For now return as is
+        return geometry
+    if len(grid_mapping_names) > 1:
+        raise ValueError(f"Multiple grid mappings found: {grid_mapping_names!r}!")
+    (grid_mapping_var,) = tuple(itertools.chain(*ds.cf.grid_mapping_names.values()))
+
+    grid_mapping = ds[grid_mapping_var]
+    data_crs = pyproj.crs.CRS.from_cf(grid_mapping.attrs)
+
+    transformer = transformer_from_crs(
+        crs_from=geometry_crs,
+        crs_to=data_crs,
+        always_xy=True,
+    )
+    return transform(transformer.transform, geometry)

xpublish_edr/geometry/position.py

@@ -45,7 +45,7 @@ def _select_by_position_regular_xy_grid(
     """
     # Find the nearest X and Y coordinates to the point
     if method == "nearest":
-        return ds.cf.sel(X=point.x, Y=point.y, method=method)
+        return ds.cf.sel(X=[point.x], Y=[point.y], method=method)
     else:
         return ds.cf.interp(X=point.x, Y=point.y, method=method)
 

xpublish_edr/plugin.py

@@ -101,7 +101,7 @@ def get_position(
             logger.debug(f"Dataset filtered by query params {ds}")
 
             try:
-                ds = select_by_position(ds, query.geometry, query.method)
+                ds = select_by_position(ds, query.project_geometry(ds), query.method)
             except KeyError:
                 raise HTTPException(
                     status_code=404,
@@ -148,7 +148,7 @@ def get_area(
             logger.debug(f"Dataset filtered by query params {ds}")
 
             try:
-                ds = select_by_area(ds, query.geometry)
+                ds = select_by_area(ds, query.project_geometry(ds))
             except KeyError:
                 raise HTTPException(
                     status_code=404,

xpublish_edr/query.py

@@ -7,8 +7,9 @@
 import xarray as xr
 from fastapi import Query
 from pydantic import BaseModel, Field
-from shapely import wkt
+from shapely import Geometry, wkt
 
+from xpublish_edr.geometry.common import project_geometry
 from xpublish_edr.logger import logger
 
 
@@ -25,15 +26,24 @@ class EDRQuery(BaseModel):
     z: Optional[str] = None
     datetime: Optional[str] = None
     parameters: Optional[str] = None
-    crs: Optional[str] = None
+    crs: str = Field(
+        "EPSG:4326",
+        title="Coordinate Reference System",
+        description="Coordinate Reference System for the query. Default is EPSG:4326",
+    )
     format: Optional[str] = None
     method: Literal["nearest", "linear"] = "nearest"
 
     @property
-    def geometry(self):
+    def geometry(self) -> Geometry:
         """Shapely point from WKT query params"""
         return wkt.loads(self.coords)
 
+    def project_geometry(self, ds: xr.Dataset) -> Geometry:
+        """Project the geometry to the dataset's CRS"""
+        geometry = self.geometry
+        return project_geometry(ds, self.crs, geometry)
+
     def select(self, ds: xr.Dataset, query_params: dict) -> xr.Dataset:
         """Select data from a dataset based on the query"""
         if self.z:
@@ -117,8 +127,8 @@ def edr_query(
         alias="parameter-name",
         description="xarray variables to query",
     ),
-    crs: Optional[str] = Query(
-        None,
+    crs: str = Query(
+        "EPSG:4326",
         deprecated=True,
         description="CRS is not yet implemented",
     ),