feat: Accept POST request (#105)

* Add post support for position queries

* Improve perfoormance of multipoint reprojection

* Add Area query POST support, update README

* Fix up tests

* Some cleanup

Author: mpiannucci

README.md

@@ -52,6 +52,8 @@ rest = xpublish.Rest(
 
 This package attempts to follow [the spec](https://docs.ogc.org/is/19-086r6/19-086r6.html) where reasonable, adding functionality where the value is demonstrable.
 
+> **Note:** `POST` is supported on `/position` and `/area` as a non-spec extension so that requests with large geometries (many points, complex polygons) can submit them in the request body instead of being limited by URL length. All selection parameters (`datetime`, `z`, `parameter-name`, `crs`, `f`, `method`) are still passed as query string parameters. See the per-query tables below for supported body content types.
+
 ### [collections](https://docs.ogc.org/is/19-086r6/19-086r6.html#_e55ba0f5-8f24-4f1b-a7e3-45775e39ef2e) and Resource Paths Support
 
 `xpublish-edr` does not currently support the `/collections/{collectionId}/query` path template described in the spec. Instead the path resource appears as `/{dataset_id}/edr/{query}`. This is because of the path structure of xpublish. In the future, if `xpublish` supports [`DataTree`](https://docs.xarray.dev/en/stable/generated/xarray.DataTree.html) it could provide a path to supporting the spec compliant `collections` resource path.
@@ -64,27 +66,29 @@ This package attempts to follow [the spec](https://docs.ogc.org/is/19-086r6/19-0
 
 | Query  | Compliant | Comments
 | ------------- | ------------- | ------------- |
-| `coords`  | ✅ | |
+| `coords`  | ✅ | Required for `GET`; for `POST` the points are read from the request body |
 | `z`  | ✅ | |
 | `datetime`  | ✅ | |
 | `parameter-name`  | ✅   | |
 | `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` |
 | `f`  | ✅ | Supports `cf_covjson`, `csv`, `geojson` `netcdf`, `parquet` |
 | `method`  | ➕ | Optional: controls data selection. Use "nearest" for nearest neighbor selection, or "linear" for interpolated selection. Uses `nearest` if not specified |
+| `POST` body  | ➕ | Non-spec extension. Supported content types: `text/csv` (columns `x`/`y`, `lon`/`lat`, or `longitude`/`latitude`); `application/geo+json` (Point, MultiPoint, Feature, FeatureCollection, or GeometryCollection) |
 
 > Any additional query parameters are assumed to be additional selections to make on the dimensions/coordinates. These queries will use the specified selections `method`.
 
 [8.2.3 Area query](https://docs.ogc.org/is/19-086r6/19-086r6.html#_c92d1888-dc80-454f-8452-e2f070b90dcd)
 
 | Query  | Compliant | Comments
 | ------------- | ------------- | ------------- |
-| `coords`  | ✅ | Only `POLYGON` supported currently |
+| `coords`  | ✅ | `POLYGON` and `MULTIPOLYGON` supported. Required for `GET`; for `POST` the polygon is read from the request body |
 | `z`  | ✅   | |
 | `datetime`  | ✅ | |
 | `parameter-name`  | ✅   | |
 | `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` |
 | `f`  | ✅   | Supports `cf_covjson`, `csv`, `geojson` `netcdf`, `parquet` |
 | `method`  | ➕ | Optional: controls data selection. Use "nearest" for nearest neighbor selection, or "linear" for interpolated selection. Uses `nearest` if not specified |
+| `POST` body  | ➕ | Non-spec extension. Supported content types: `application/geo+json` (Polygon, MultiPolygon, Feature, FeatureCollection, or GeometryCollection); `application/wkt` / `text/plain` (raw WKT Polygon or MultiPolygon) |
 
 > `method` is not applicable for the coordinates of area queries, only for selecting datetime, z, or additional dimensions.
 

tests/test_cf_router.py

@@ -1,3 +1,4 @@
+import json
 from io import BytesIO
 
 import cf_xarray  # noqa: F401
@@ -615,6 +616,45 @@ def test_cf_multiple_position_csv(cf_client):
         assert key in csv_data[0], f"column {key} should be in the header"
 
 
+def test_cf_post_position_csv_body(cf_client):
+    response = cf_client.post(
+        "/datasets/air/edr/position",
+        content="lon,lat\n202,43\n205,45\n",
+        headers={"content-type": "text/csv"},
+    )
+    assert response.status_code == 200, response.text
+    axes = response.json()["domain"]["axes"]
+    assert axes["x"] == {"values": [202.5, 205.0]}
+    assert axes["y"] == {"values": [42.5, 45.0]}
+
+
+def test_cf_post_position_geojson_body(cf_client):
+    body = {
+        "type": "FeatureCollection",
+        "features": [
+            {
+                "type": "Feature",
+                "geometry": {"type": "Point", "coordinates": [202, 43]},
+                "properties": {},
+            },
+            {
+                "type": "Feature",
+                "geometry": {"type": "Point", "coordinates": [205, 45]},
+                "properties": {},
+            },
+        ],
+    }
+    response = cf_client.post(
+        "/datasets/air/edr/position",
+        json=body,
+        headers={"content-type": "application/geo+json"},
+    )
+    assert response.status_code == 200, response.text
+    axes = response.json()["domain"]["axes"]
+    assert axes["x"] == {"values": [202.5, 205.0]}
+    assert axes["y"] == {"values": [42.5, 45.0]}
+
+
 def test_cf_area_query(cf_client, cf_air_dataset):
     coords = "POLYGON((201 41, 201 49, 209 49, 209 41, 201 41))"
     response = cf_client.get(f"/datasets/air/edr/area?coords={coords}&f=cf_covjson")
@@ -716,6 +756,56 @@ def test_cf_area_geojson_query(cf_client, cf_air_dataset):
     assert len(features) == 36, "There should be 36 data points"
 
 
+@pytest.mark.parametrize(
+    "content_type,body,expected_status",
+    [
+        (
+            "application/geo+json",
+            json.dumps(
+                {
+                    "type": "Feature",
+                    "geometry": {
+                        "type": "Polygon",
+                        "coordinates": [
+                            [[201, 41], [201, 49], [209, 49], [209, 41], [201, 41]],
+                        ],
+                    },
+                    "properties": {},
+                },
+            ),
+            200,
+        ),
+        (
+            "application/wkt",
+            "POLYGON((201 41, 201 49, 209 49, 209 41, 201 41))",
+            200,
+        ),
+        ("application/geo+json", "", 422),
+        (
+            "application/geo+json",
+            json.dumps({"type": "Point", "coordinates": [202, 43]}),
+            422,
+        ),
+    ],
+    ids=["geojson_feature", "wkt", "empty_body", "point_rejected"],
+)
+def test_cf_post_area_body(cf_client, content_type, body, expected_status):
+    response = cf_client.post(
+        "/datasets/air/edr/area",
+        content=body,
+        headers={"content-type": content_type},
+    )
+    assert response.status_code == expected_status, response.text
+    if expected_status == 200:
+        axes = response.json()["domain"]["axes"]
+        assert axes["x"] == {
+            "values": [202.5, 205.0, 207.5, 202.5, 205.0, 207.5, 202.5, 205.0, 207.5],
+        }
+        assert axes["y"] == {
+            "values": [47.5, 47.5, 47.5, 45.0, 45.0, 45.0, 42.5, 42.5, 42.5],
+        }
+
+
 def test_cf_area_nc_query(cf_client, cf_air_dataset):
     coords = "POLYGON((201 41, 201 49, 209 49, 209 41, 201 41))"
     response = cf_client.get(f"/datasets/air/edr/area?coords={coords}&f=nc")

tests/test_select.py

@@ -208,7 +208,7 @@ def test_select_query_error(regular_xy_dataset):
         parameters="air",
     )
 
-    with pytest.raises(TypeError):
+    with pytest.raises(ValueError, match="Invalid datetime"):
         query.select(regular_xy_dataset, {})
 
     query = EDRPositionQuery(

xpublish_edr/geometry/common.py

@@ -6,11 +6,12 @@
 from functools import lru_cache, partial
 from typing import Union
 
+import numpy as np
 import pyproj
 import rioxarray  # noqa
+import shapely
 import xarray as xr
 from shapely import Geometry
-from shapely.ops import transform
 
 VECTORIZED_DIM = "pts"
 
@@ -93,7 +94,13 @@ def project_geometry(ds: xr.Dataset, geometry_crs: str, geometry: Geometry) -> G
         crs_from=geometry_crs,
         crs_to=data_crs,
     )
-    return transform(transformer.transform, geometry)
+
+    def _transform(coords: np.ndarray) -> np.ndarray:
+        """Vectorized callback for shapely.transform: project all coords in one call."""
+        x, y = transformer.transform(coords[:, 0], coords[:, 1])
+        return np.column_stack([x, y])
+
+    return shapely.transform(geometry, _transform)
 
 
 def project_bbox(