From cf00eeb51864353f8e4b8396a0eedb7d9f1a024c Mon Sep 17 00:00:00 2001
From: Matthew Iannucci <matt.iannucci@rpsgroup.com>
Date: Wed, 16 Oct 2024 15:12:17 -0400
Subject: [PATCH] Normalize 4326 coordinates (#84)

* Normalize regular grid projected coordinates

* Normalize coordinates

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Remove unused method

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
---
 xpublish_wms/grids/grid.py           |  7 +++++--
 xpublish_wms/grids/regular.py        | 22 ++++++++++++++++++++++
 xpublish_wms/wms/__init__.py         |  1 +
 xpublish_wms/wms/get_feature_info.py |  3 +++
 4 files changed, 31 insertions(+), 2 deletions(-)

diff --git a/xpublish_wms/grids/grid.py b/xpublish_wms/grids/grid.py
index 5fbb587..3c09728 100644
--- a/xpublish_wms/grids/grid.py
+++ b/xpublish_wms/grids/grid.py
@@ -133,9 +133,12 @@ def sel_lat_lng(
         parameters,
     ) -> tuple[xr.Dataset, list, list]:
         """Select the given dataset by the given lon/lat and optional elevation"""
+        subset = self.mask(subset[parameters])
 
-        subset = self.mask(subset)
-        subset = subset.cf.interp(longitude=lng, latitude=lat)
+        subset = subset.cf.interp(
+            longitude=lng,
+            latitude=lat,
+        )
 
         x_axis = [strip_float(subset.cf["longitude"])]
         y_axis = [strip_float(subset.cf["latitude"])]
diff --git a/xpublish_wms/grids/regular.py b/xpublish_wms/grids/regular.py
index 9cf873b..7e0304f 100644
--- a/xpublish_wms/grids/regular.py
+++ b/xpublish_wms/grids/regular.py
@@ -40,6 +40,11 @@ def project(
             if coord != da.cf["longitude"].name and coord != da.cf["latitude"].name:
                 coords[coord] = da.coords[coord]
 
+        # normalize longitude to be between -180 and 180
+        da = da.cf.assign_coords(
+            longitude=(((da.cf["longitude"] + 180) % 360) - 180),
+        ).sortby(da.cf["longitude"].name)
+
         # build new x coordinate
         coords["x"] = ("x", da.cf["longitude"].values, da.cf["longitude"].attrs)
         # build new y coordinate
@@ -61,3 +66,20 @@ def project(
             da = da.unify_chunks()
 
         return da, None, None
+
+    def sel_lat_lng(
+        self,
+        subset: xr.Dataset,
+        lng,
+        lat,
+        parameters,
+    ) -> tuple[xr.Dataset, list, list]:
+        # normalize longitude to be between -180 and 180
+        subset = subset.cf.assign_coords(
+            longitude=(((subset.cf["longitude"] + 180) % 360) - 180),
+        ).sortby(subset.cf["longitude"].name)
+
+        if lng > 180:
+            lng = (lng + 180) % 360 - 180
+
+        return super().sel_lat_lng(subset, lng, lat, parameters)
diff --git a/xpublish_wms/wms/__init__.py b/xpublish_wms/wms/__init__.py
index fbcf68b..f2def79 100644
--- a/xpublish_wms/wms/__init__.py
+++ b/xpublish_wms/wms/__init__.py
@@ -30,6 +30,7 @@ def wms_handler(
     method = query_params.get("request", "").lower()
     logger.debug(f"Received wms request: {request.url}")
     logger.info(f"WMS: {method}")
+
     if method == "getcapabilities":
         return get_capabilities(dataset, request, query_params)
     elif method == "getmap":
diff --git a/xpublish_wms/wms/get_feature_info.py b/xpublish_wms/wms/get_feature_info.py
index 66c9dae..3b4bbf3 100644
--- a/xpublish_wms/wms/get_feature_info.py
+++ b/xpublish_wms/wms/get_feature_info.py
@@ -134,6 +134,9 @@ def get_feature_info(ds: xr.Dataset, query: dict) -> Response:
     any_has_vertical_axis = True in has_vertical_axis
 
     crs = query.get("crs", None) or query.get("srs")
+    if crs != "EPSG:4326":
+        raise HTTPException(501, "Only EPSG:4326 is supported")
+
     bbox = [float(x) for x in query["bbox"].split(",")]
     width = int(query["width"])
     height = int(query["height"])