Fix GFI for fvcom

xpublish_wms/grids/fvcom.py

@@ -1,4 +1,4 @@
-from typing import Optional, Sequence
+from typing import Optional, Sequence, Union
 
 import dask.array as dask_array
 import matplotlib.tri as tri
@@ -9,6 +9,7 @@
 from xpublish_wms.utils import (
     barycentric_weights,
     lat_lng_find_tri,
+    lat_lng_find_tri_ind,
     strip_float,
     to_mercator,
 )
@@ -98,6 +99,18 @@ def sel_lat_lng(
         elif "siglev" in subset.dims:
             subset = subset.rename_dims({"siglev": "sigmaa"})
 
+        if "nele" in subset.dims:
+            return self.sel_lat_lng_nele(subset, lng, lat, parameters)
+        else:
+            return self.sel_lat_lng_node(subset, lng, lat, parameters)
+
+    def sel_lat_lng_node(
+        self,
+        subset: xr.Dataset,
+        lng,
+        lat,
+        parameters,
+    ) -> tuple[xr.Dataset, list, list]:
         temp_arrays = dict()
         # create new dataarrays so that nele variables can be adjusted appropriately
         for parameter in parameters:
@@ -185,9 +198,14 @@ def sel_lat_lng(
         lng_values = subset.cf["longitude"].values
         lat_values = subset.cf["latitude"].values
 
+        if lng < 0 and np.min(lng_values) > 0:
+            lng += 360
+        elif lng > 180 and np.min(lng_values) < 0:
+            lng -= 360
+
         # find if the selected lng/lat is within a triangle
         valid_tri = lat_lng_find_tri(
-            lng + 360,
+            lng,
             lat,
             lng_values,
             lat_values,
@@ -209,7 +227,7 @@ def sel_lat_lng(
             p1 = [lng_values[valid_tri[0]], lat_values[valid_tri[0]]]
             p2 = [lng_values[valid_tri[1]], lat_values[valid_tri[1]]]
             p3 = [lng_values[valid_tri[2]], lat_values[valid_tri[2]]]
-            w1, w2, w3 = barycentric_weights([lng + 360, lat], p1, p2, p3)
+            w1, w2, w3 = barycentric_weights([lng, lat], p1, p2, p3)
 
             for parameter in parameters:
                 values = subset[parameter].values
@@ -232,6 +250,50 @@ def sel_lat_lng(
         y_axis = [strip_float(ret_subset.cf["latitude"])]
         return ret_subset, x_axis, y_axis
 
+    def sel_lat_lng_nele(
+        self,
+        subset: xr.Dataset,
+        lng,
+        lat,
+        parameters,
+    ) -> tuple[xr.Dataset, list, list]:
+        lng_values = self.ds.lon.values
+        lat_values = self.ds.lat.values
+
+        if lng < 0 and np.max(lng_values) > 180:
+            lng += 360
+        elif lng > 180 and np.min(lng_values) < 0:
+            lng -= 360
+
+        # find if the selected lng/lat is within a triangle
+        valid_tri = lat_lng_find_tri_ind(
+            lng,
+            lat,
+            lng_values,
+            lat_values,
+            self.tessellate(subset),
+        )
+
+        ret_subset = subset.isel(nele=0)
+
+        # if no -> set all values to nan
+        if valid_tri is None:
+            for parameter in parameters:
+                ret_subset.__setitem__(
+                    parameter,
+                    (
+                        ret_subset[parameter].dims,
+                        np.full(ret_subset[parameter].shape, np.nan),
+                        ret_subset[parameter].attrs,
+                    ),
+                )
+        else:
+            ret_subset = subset.isel(nele=valid_tri[0])
+
+        x_axis = [strip_float(ret_subset.cf["longitude"])]
+        y_axis = [strip_float(ret_subset.cf["latitude"])]
+        return ret_subset, x_axis, y_axis
+
     def select_by_elevation(
         self,
         da: xr.DataArray,
@@ -371,7 +433,7 @@ def project(
 
         return da, None, None
 
-    def tessellate(self, da: xr.DataArray) -> np.ndarray:
+    def tessellate(self, da: Union[xr.DataArray, xr.Dataset]) -> np.ndarray:
         nv = self.ds.nv
         if len(nv.shape) > 2:
             for i in range(len(nv.shape) - 2):

xpublish_wms/utils.py

@@ -226,9 +226,9 @@ def barycentric_weights(point, v1, v2, v3):
     return w1, w2, w3
 
 
-def lat_lng_find_tri(lng, lat, lng_values, lat_values, triangles):
+def lat_lng_find_tri_ind(lng, lat, lng_values, lat_values, triangles):
     """
-    Find the triangle that the inputted lng/lat is within
+    Find the triangle index that the inputted lng/lat is within
 
     Inputs
     ------
@@ -248,7 +248,6 @@ def lat_lng_find_tri(lng, lat, lng_values, lat_values, triangles):
     to interpolate the value accurate to the lng/lat requested, or None if the point is outside
     the triangular mesh
     """
-
     lnglat_data = np.stack((lng_values[triangles], lat_values[triangles]), axis=2)
 
     d1 = (
@@ -270,7 +269,36 @@ def lat_lng_find_tri(lng, lat, lng_values, lat_values, triangles):
     if len(tri_index) == 0 or len(tri_index[0]) == 0:
         return None
     else:
-        return triangles[tri_index[0]][0]
+        return tri_index[0]
+
+
+def lat_lng_find_tri(lng, lat, lng_values, lat_values, triangles):
+    """
+    Find the triangle that the inputted lng/lat is within
+
+    Inputs
+    ------
+    lng: float, int
+        Longitude of comparison point.
+    lat: float, int
+        Latitude of comparison point.
+    lng_values: xr.DataArray
+        Longitudes of points corresponding to the indices in triangles
+    lat_values: xr.DataArray
+        Latitudes of points corresponding to the indices in triangles
+    triangles: ndarray of shape (X, 3)
+        Triangle mesh of indices as generated by tri.Triangulation
+    Returns
+    -------
+    Triangle of indices that the lng/lat is within, which can be used with barycentric_weights
+    to interpolate the value accurate to the lng/lat requested, or None if the point is outside
+    the triangular mesh
+    """
+    tri_index = lat_lng_find_tri_ind(lng, lat, lng_values, lat_values, triangles)
+    if tri_index is None:
+        return None
+    else:
+        return triangles[tri_index][0]
 
 
 def bilinear_interp(percent_point, percent_quad, value_quad):