Vertical level support for API

api/formatters/covjson.py

@@ -27,7 +27,7 @@
 
 # mime_type = "application/prs.coverage+json"
 
-Dom = namedtuple("Dom", ["lat", "lon", "times"])
+Dom = namedtuple("Dom", ["lat", "lon", "z", "times"])
 Data = namedtuple("Data", ["dom", "values", "ts_mdata"])
 
 
@@ -73,12 +73,16 @@ def convert_to_covjson(observations):
 
     # Need to sort before using groupBy. Also sort on parameter_name to get consistently sorted output
     data.sort(key=lambda x: (x.dom, x.ts_mdata.parameter_name))
-    for (lat, lon, times), group in groupby(data, lambda x: x.dom):
+    for (lat, lon, z, times), group in groupby(data, lambda x: x.dom):
         referencing = [
             ReferenceSystemConnectionObject(
                 coordinates=["x", "y"],
                 system=ReferenceSystem(type="GeographicCRS", id="http://www.opengis.net/def/crs/OGC/1.3/CRS84"),
             ),
+            ReferenceSystemConnectionObject(
+                coordinates=["z"],
+                system=ReferenceSystem(type="VerticalCRS", id="http://www.opengis.net/def/crs/EPSG/0/5715"),
+            ),
             ReferenceSystemConnectionObject(
                 coordinates=["t"],
                 system=ReferenceSystem(type="TemporalRS", calendar="Gregorian"),
@@ -89,6 +93,7 @@ def convert_to_covjson(observations):
             axes=Axes(
                 x=ValuesAxis[float](values=[lon]),
                 y=ValuesAxis[float](values=[lat]),
+                z=ValuesAxis[float](values=[z]),
                 t=ValuesAxis[AwareDatetime](values=times),
             ),
             referencing=referencing,
@@ -106,7 +111,7 @@ def convert_to_covjson(observations):
             parameters[parameter_id] = make_parameter(data.ts_mdata)
 
             ranges[parameter_id] = NdArrayFloat(
-                values=values_no_nan, axisNames=["t", "x", "y"], shape=[len(values_no_nan), 1, 1]
+                values=values_no_nan, axisNames=["t", "x", "y", "z"], shape=[len(values_no_nan), 1, 1, 1]
             )
 
         custom_fields = {"metocean:wigosId": data.ts_mdata.platform}
@@ -124,9 +129,12 @@ def convert_to_covjson(observations):
 def _collect_data(ts_mdata, obs_mdata):
     lat = obs_mdata[0].geo_point.lat  # HACK: For now assume they all have the same position
     lon = obs_mdata[0].geo_point.lon
+    # HACK: For now assume they all have the same elevation
+    # Unset int64 defaults to 0, which is a valid elevation for pointseries
+    z = convert_cm_to_m(obs_mdata[0].camsl)
     tuples = (
         (o.obstime_instant.ToDatetime(tzinfo=timezone.utc), float(o.value)) for o in obs_mdata
-    )  # HACK: str -> float
+    )  # HACK: str -> float))
     (times, values) = zip(*tuples)
 
-    return Data(Dom(lat, lon, times), values, ts_mdata)
+    return Data(Dom(lat, lon, z, times), values, ts_mdata)

api/formatters/geojson.py

@@ -51,6 +51,7 @@ def convert_to_geojson(observations, base_url):
                 coordinates=[
                     ts.obs_mdata[0].geo_point.lon,
                     ts.obs_mdata[0].geo_point.lat,
+                    convert_cm_to_m(ts.obs_mdata[0].camsl),
                 ],
             ),
         )

api/metadata_endpoints.py

@@ -16,6 +16,8 @@
 from edr_pydantic.extent import Custom
 from edr_pydantic.extent import Extent
 from edr_pydantic.extent import Spatial
+
+# from edr_pydantic.extent import Vertical
 from edr_pydantic.extent import Temporal
 from edr_pydantic.link import EDRQueryLink
 from edr_pydantic.link import Link
@@ -176,11 +178,13 @@ async def get_collection_metadata(base_url: str, collection_id: str, is_self) ->
 
         all_parameters[ts.parameter_name] = parameter
 
+    # TODO: Get the vertical extent values with GetExtentsRequest
     extent_request = dstore.GetExtentsRequest()
     extent_response = await get_extents_request(extent_request)
     spatial_extent = extent_response.spatial_extent
     interval_start = extent_response.temporal_extent.start.ToDatetime(tzinfo=timezone.utc)
     interval_end = extent_response.temporal_extent.end.ToDatetime(tzinfo=timezone.utc)
+    # vertical_extent = extent_response.vertical_extent
 
     # TODO: Check if these make /collections significantly slower. If yes, do we need DB indices on these? And parallel
     levels = [convert_cm_to_m(level) for level in await get_unique_values_for_metadata("level")]
@@ -212,6 +216,12 @@ async def get_collection_metadata(base_url: str, collection_id: str, is_self) ->
                 values=[f"{datetime_to_iso_string(interval_start)}/{datetime_to_iso_string(interval_end)}"],
                 trs="Gregorian",
             ),
+            # TODO: Add vertical extent to the response when we get the values
+            # vertical=Vertical(
+            #    interval=[[vertical_extent[0], vertical_extent[-1]]],
+            #    values=[vertical_extent],
+            #    vrs="EPSG:5714",
+            # ),
             custom=[
                 Custom(
                     id="standard_name",

api/openapi/edr_query_parameter_descriptions.py

@@ -1,5 +1,10 @@
 bbox = "Bounding box to query data from. The bounding box is an area defined by two longitudes and two latitudes."
 datetime = "Time range to query data from."
+z = (
+    "Vertical level(s) to query data from in meters. Can be a comma separated list, a range or a repeating interval. "
+    "Repeating intervals are defined in the format of "
+    "'__R__ *number of intervals / min-height / height to increment by*'."
+)
 parameter_name = (
     "Comma separated list of parameter names. Each consists of four components separated by colons."
     " The components are standard name, level in meters, aggregation method, and period. "
@@ -8,7 +13,7 @@
 )
 standard_name = "Comma separated list of parameter standard_name(s) to query."
 level = (
-    "Define the vertical level(s) to return data from using either a comma separated list, "
+    "Define the standard measurement level(s) to return data from using either a comma separated list, "
     "a range or a repeating interval. <br /> Repeating intervals are defined in the format of "
     "'__R__ *number of intervals / min-level / height to increment by*'."
 )

api/openapi/openapi_examples.py

@@ -109,9 +109,20 @@
     "10.0": {"summary": "10 meters above ground", "value": "10.0"},
 }
 
-period = {
+duration = {
     "PT1M": {"summary": "1 minute", "value": "PT1M"},
     "PT10M": {"summary": "10 minutes", "value": "PT10M"},
 }
 
 method = {"mean": {"summary": "Mean", "value": "mean"}, "maximum": {"summary": "Maximum", "value": "maximum"}}
+
+z = {
+    "Default": {"value": ""},
+    "10.0": {"summary": "10 meters above mean sea level", "value": "10.0"},
+    "5/200": {"summary": "Range from 5 to 200 meters above mean sea level", "value": "5/200"},
+    "R10/5/10": {
+        "summary": "Repeating interval of 10 levels, starting at 5 meters above mean sea level, "
+        "with an increment of 10 meters",
+        "value": "R10/5/10",
+    },
+}

api/routers/edr.py

@@ -44,6 +44,7 @@
     "unit",
     "obstime_instant",
     "value",
+    "camsl",
 ]
 
 
@@ -72,6 +73,10 @@ async def get_locations(
             openapi_examples=openapi_examples.datetime,
         ),
     ] = None,
+    z: Annotated[
+        str | None,
+        Query(description=edr_query_parameter_descriptions.z, openapi_examples=openapi_examples.z),
+    ] = None,
     parameter_name: Annotated[
         str | None,
         Query(
@@ -119,7 +124,7 @@ async def get_locations(
             [dstore.Point(lat=coord[1], lon=coord[0]) for coord in poly.exterior.coords],
         )
 
-    await add_request_parameters(loc_request, parameter_name, datetime, standard_name, level, method, duration)
+    await add_request_parameters(loc_request, parameter_name, datetime, z, standard_name, level, method, duration)
 
     grpc_response = await get_locations_request(loc_request)
     locations = grpc_response.locations
@@ -145,6 +150,8 @@ async def get_locations(
                     + "collections/observations/items?platform="
                     + loc.platform,
                 },
+                # TODO: loc_request does not return elevation for locations currently
+                # so for now we just leave it out entirely untill #256 is done
                 geometry=Point(
                     type="Point",
                     coordinates=(loc.geo_point.lon, loc.geo_point.lat),
@@ -188,6 +195,10 @@ async def get_data_location_id(
         str | None,
         Query(description=edr_query_parameter_descriptions.datetime, openapi_examples=openapi_examples.datetime),
     ] = None,
+    z: Annotated[
+        str | None,
+        Query(description=edr_query_parameter_descriptions.z, openapi_examples=openapi_examples.z),
+    ] = None,
     f: Annotated[
         formatters.Formats, Query(description=edr_query_parameter_descriptions.format)
     ] = formatters.Formats.covjson,
@@ -227,7 +238,7 @@ async def get_data_location_id(
         included_response_fields=response_fields_needed_for_data_api,
     )
 
-    await add_request_parameters(request, parameter_name, datetime, standard_name, level, method, duration)
+    await add_request_parameters(request, parameter_name, datetime, z, standard_name, level, method, duration)
 
     grpc_response = await get_obs_request(request)
     observations = grpc_response.observations
@@ -262,6 +273,10 @@ async def get_data_position(
             openapi_examples=openapi_examples.datetime,
         ),
     ] = None,
+    z: Annotated[
+        str | None,
+        Query(description=edr_query_parameter_descriptions.z, openapi_examples=openapi_examples.z),
+    ] = None,
     f: Annotated[
         formatters.Formats, Query(description=edr_query_parameter_descriptions.format)
     ] = formatters.Formats.covjson,
@@ -314,13 +329,23 @@ async def get_data_position(
             detail={"coords": f"Unexpected error occurred during wkt parsing: {coords}"},
         )
 
+    # TODO: Allow elevation in Point and update backend to use 3D spatial indexing?
     request = dstore.GetObsRequest(
         # 10 meters around the point
         spatial_circle=dstore.Circle(center=dstore.Point(lat=point.y, lon=point.x), radius=0.01),
         included_response_fields=response_fields_needed_for_data_api,
     )
 
-    await add_request_parameters(request, parameter_name, datetime, standard_name, level, method, duration)
+    # OGC EDR SPEC states:
+    # If a client request has a coords value which includes a height value and defines a z query parameter,
+    # the z query parameter will be the requested height value.
+    elevation = None
+    if point.has_z:
+        elevation = str(point.z)
+    if z:
+        elevation = z
+
+    await add_request_parameters(request, parameter_name, datetime, elevation, standard_name, level, method, duration)
 
     grpc_response = await get_obs_request(request)
     observations = grpc_response.observations
@@ -352,6 +377,10 @@ async def get_data_area(
         str | None,
         Query(description=edr_query_parameter_descriptions.datetime, openapi_examples=openapi_examples.datetime),
     ] = None,
+    z: Annotated[
+        str | None,
+        Query(description=edr_query_parameter_descriptions.z, openapi_examples=openapi_examples.z),
+    ] = None,
     f: Annotated[
         formatters.Formats, Query(description=edr_query_parameter_descriptions.format)
     ] = formatters.Formats.covjson,
@@ -404,14 +433,27 @@ async def get_data_area(
             detail={"coords": f"Unexpected error occurred during wkt parsing: {coords}"},
         )
 
+    # TODO: Allow elevation in Polygon and update backend to use 3D spatial indexing?
     request = dstore.GetObsRequest(
         spatial_polygon=dstore.Polygon(
             points=[dstore.Point(lat=coord[1], lon=coord[0]) for coord in poly.exterior.coords]
         ),
         included_response_fields=response_fields_needed_for_data_api,
     )
 
-    await add_request_parameters(request, parameter_name, datetime, standard_name, level, method, duration)
+    # TODO: This does not explicitly handle the case of a 3D polygon, but instead uses the ranges
+    # from poly.exterior.coords and creates a 3D flat polygon from the lowest to the highest z value.
+    elevation = None
+    if len(poly.exterior.coords[0]) == 3:
+        z_values = [coord[2] for coord in poly.exterior.coords]
+        elevation = f"{min(z_values)}/{max(z_values)}"
+    # OGC EDR SPEC states:
+    # If a client request has a coords value which includes a height value and defines a z query parameter,
+    # the z query parameter will be the requested height value.
+    if z:
+        elevation = z
+
+    await add_request_parameters(request, parameter_name, datetime, elevation, standard_name, level, method, duration)
 
     grpc_response = await get_obs_request(request)
     observations = grpc_response.observations

api/routers/feature.py

@@ -21,6 +21,8 @@
 from shapely import geometry
 from utilities import get_datetime_range
 from utilities import split_and_strip
+from utilities import get_levels_values
+from utilities import get_durations_or_range
 
 router = APIRouter(prefix="/collections/observations")
 
@@ -49,6 +51,13 @@ async def search_timeseries(
             description="E-SOH database only contains data from the last 24 hours",
         ),
     ] = None,
+    z: Annotated[
+        str | None,
+        Query(
+            description="Vertical level in meters above or below mean sea level",
+            openapi_examples=openapi_examples.z,
+        ),
+    ] = None,
     id: Annotated[str | None, Query(description="E-SOH time series id")] = None,
     parameter_name: Annotated[str | None, Query(alias="parameter-name", description="E-SOH parameter name")] = None,
     naming_authority: Annotated[
@@ -85,9 +94,9 @@ async def search_timeseries(
             openapi_examples=openapi_examples.level,
         ),
     ] = None,
-    period: Annotated[
+    duration: Annotated[
         str | None,
-        Query(description="Duration of collection period in ISO8601", openapi_examples=openapi_examples.period),
+        Query(description="Duration of collection period in ISO8601", openapi_examples=openapi_examples.duration),
     ] = None,
     method: Annotated[
         str | None,
@@ -117,9 +126,10 @@ async def search_timeseries(
             standard_name=dstore.Strings(values=split_and_strip(standard_name) if standard_name else None),
             unit=dstore.Strings(values=split_and_strip(unit) if unit else None),
             instrument=dstore.Strings(values=split_and_strip(instrument) if instrument else None),
-            level=dstore.Strings(values=split_and_strip(level) if level else None),
-            period=dstore.Strings(values=split_and_strip(period) if period else None),
+            level=dstore.Strings(values=get_levels_values(level) if level else None),
+            period=dstore.Strings(values=get_durations_or_range(duration) if duration else None),
             function=dstore.Strings(values=split_and_strip(method) if method else None),
+            camsl=dstore.Strings(values=get_levels_values(z) if z else None),
         ),
         spatial_polygon=(
             dstore.Polygon(points=[dstore.Point(lat=coord[1], lon=coord[0]) for coord in poly.exterior.coords])
@@ -163,6 +173,7 @@ async def get_dataset_metadata(request: Request):
 
     # need to get spatial extent.
     spatial_request = dstore.GetExtentsRequest()
+    # TODO: Add vertical extent to the response in datastore and add it to the jinja tempalate
     extent = await get_extents_request(spatial_request)
     dynamic_fields = {
         "spatial_extents": [