Fix coordinate discontinuities on x- and y-coordinates for 2D grids

examples/maplibre/tiles.html

@@ -61,7 +61,7 @@
                     type: "raster",
                     // use the tiles option to specify a raster tile source URL
                     // https://docs.mapbox.com/style-spec/reference/sources/
-                    url: "http://localhost:8080/tiles/WebMercatorQuad/tilejson.json?variables=gust&width=512&height=512&style=raster/viridis&colorscalerange=0,30",
+                    url: "http://localhost:8080/tiles/WebMercatorQuad/tilejson.json?variables=foo&width=512&height=512&style=raster/viridis&colorscalerange=-1,1",
                     tileSize: 512,
                 });
                 map.addLayer(

pyproject.toml

@@ -88,6 +88,7 @@ dev = [
     "xpublish-tiles[testing]",
     "coiled>=1.118.3",
     "memray>=1.18.0",
+    "pre-commit>=4.5.1",
 ]
 
 [project.entry-points."xpublish.plugin"]

src/xpublish_tiles/pipeline.py

@@ -12,6 +12,7 @@
 import xarray as xr
 from xpublish_tiles.config import config
 from xpublish_tiles.grids import (
+    Curvilinear,
     GridMetadata,
     GridSystem,
     GridSystem2D,
@@ -426,13 +427,22 @@ def coarsen(
         # a discontinuity at the anti-meridian; which we end up averaging over below.
         # So fix that here.
         if grid.lon_spans_globe:
-            newX = fix_coordinate_discontinuities(
-                da[grid.X].data,
-                # FIXME: test 0->360 also!
-                transformer_from_crs(grid.crs, grid.crs),
-                axis=da[grid.X].get_axis_num(grid.Xdim),
-                bbox=grid.bbox,
-            )
+            if grid.Xdim in coarsen_factors:
+                newX = fix_coordinate_discontinuities(
+                    da[grid.X].data,
+                    # FIXME: test 0->360 also!
+                    transformer_from_crs(grid.crs, grid.crs),
+                    axis=da[grid.X].get_axis_num(grid.Xdim),
+                    bbox=grid.bbox,
+                )
+            if grid.Ydim in coarsen_factors:
+                newX = fix_coordinate_discontinuities(
+                    da[grid.X].data,
+                    # FIXME: test 0->360 also!
+                    transformer_from_crs(grid.crs, grid.crs),
+                    axis=da[grid.X].get_axis_num(grid.Ydim),
+                    bbox=grid.bbox,
+                )
             da = da.assign_coords({grid.X: da[grid.X].copy(data=newX)})
         with NUMBA_THREADING_LOCK:
             coarsened = da.coarsen(coarsen_factors, boundary="pad").mean()  # type: ignore[unresolved-attribute]
@@ -854,6 +864,7 @@ def subset_to_bbox(
             left=bbox.west, right=bbox.east, top=bbox.north, bottom=bbox.south
         )
         if grid.crs.is_geographic:
+            # Handle antimeridian crossing: west > east means bbox crosses -180/180
             west = west - 360 if west > east else west
 
         input_bbox = BBox(west=west, south=south, east=east, north=north)
@@ -889,11 +900,21 @@ def subset_to_bbox(
             )
         )
 
+        has_discontinuity_x = False
+        has_discontinuity_y = False
         if grid.crs.is_geographic:
             if isinstance(grid, GridSystem2D):
-                has_discontinuity = has_coordinate_discontinuity(
+                # Check for discontinuities along X dimension
+                has_discontinuity_x = has_coordinate_discontinuity(
                     subset[grid.X].data, axis=subset[grid.X].get_axis_num(grid.Xdim)
                 )
+                # For curvilinear grids (2D coordinates), also check along Y dimension
+                # (e.g., HYCOM grids can have discontinuities along Y at certain X indices)
+                if isinstance(grid, Curvilinear):
+                    has_discontinuity_y = has_coordinate_discontinuity(
+                        subset[grid.X].data, axis=subset[grid.X].get_axis_num(grid.Ydim)
+                    )
+                has_discontinuity = has_discontinuity_x or has_discontinuity_y
             elif isinstance(grid, Triangular):
                 anti = next(iter(slicers[grid.Xdim])).antimeridian_vertices
                 has_discontinuity = anti["pos"].size > 0 or anti["neg"].size > 0
@@ -913,13 +934,25 @@ def subset_to_bbox(
         # Fix coordinate discontinuities in transformed coordinates if detected
         if has_discontinuity:
             if isinstance(grid, GridSystem2D):
-                fixed = fix_coordinate_discontinuities(
-                    newX.data,
-                    input_to_output,
-                    axis=newX.get_axis_num(grid.Xdim),
-                    bbox=bbox,
-                )
-                newX = newX.copy(data=fixed)
+                # Fix discontinuities along X dimension
+                if has_discontinuity_x:
+                    fixed = fix_coordinate_discontinuities(
+                        newX.data,
+                        input_to_output,
+                        axis=newX.get_axis_num(grid.Xdim),
+                        bbox=bbox,
+                    )
+                    newX = newX.copy(data=fixed)
+
+                # For curvilinear grids, also fix discontinuities along Y dimension
+                if has_discontinuity_y and isinstance(grid, Curvilinear):
+                    fixed = fix_coordinate_discontinuities(
+                        newX.data,
+                        input_to_output,
+                        axis=newX.get_axis_num(grid.Ydim),
+                        bbox=bbox,
+                    )
+                    newX = newX.copy(data=fixed)
             elif isinstance(grid, Triangular):
                 anti = next(iter(slicers[grid.dim])).antimeridian_vertices
                 for verts in [anti["pos"], anti["neg"]]:

src/xpublish_tiles/testing/datasets.py

@@ -1059,6 +1059,114 @@ def global_nans_grid(
 )
 
 
+def _create_curvilinear_grid_like_hycom(
+    *,
+    regional_subset: bool,
+    dims: tuple[Dim, ...],
+    dtype: npt.DTypeLike,
+    attrs: dict[str, Any],
+) -> xr.Dataset:
+    """Build a global HYCOM-like curvilinear grid matching actual HYCOM/RTOFS dimensions.
+
+    Creates a simplified curvilinear grid with:
+    - Full HYCOM dimensions: 4500 (lon) x 3298 (lat)
+    - Latitude range: -80° to 90°
+    - Longitude: -180° to 180° with wraparound at antimeridian
+    - 2D coordinate arrays (curvilinear)
+    """
+    ds = uniform_grid(dims=dims, dtype=dtype, attrs=attrs)
+
+    ny, nx = ds.sizes["Y"], ds.sizes["X"]
+
+    lat_1d = np.linspace(-80.0, 90.0, ny, dtype=np.float32)
+    lon_1d = np.linspace(-180.0, 180.0, nx, dtype=np.float32)
+
+    lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)
+
+    lat_variation = 0.1 * np.sin(2 * np.pi * np.arange(nx) / nx)
+    lat = lat_2d + lat_variation[np.newaxis, :]
+
+    lon_variation = 0.5 * np.sin(2 * np.pi * np.arange(ny) / ny)
+    lon = lon_2d + lon_variation[:, np.newaxis]
+
+    lon = ((lon + 180.0) % 360.0) - 180.0
+
+    ds["foo"] = ds["foo"].chunk(X=1000, Y=1000)
+    ds.coords["latitude"] = (
+        ["Y", "X"],
+        lat.astype(np.float32),
+        {"standard_name": "latitude", "units": "degrees_north"},
+    )
+    ds.coords["longitude"] = (
+        ["Y", "X"],
+        lon.astype(np.float32),
+        {
+            "standard_name": "longitude",
+            "units": "degrees_east",
+            "modulo": "360 degrees",
+        },
+    )
+    if regional_subset:
+        mask = (
+            (ds.longitude > -180)
+            & (ds.longitude < -120)
+            & (ds.latitude > 0)
+            & (ds.latitude < 80)
+        ).compute()
+        ds.attrs["bbox"] = BBox(west=-180.0, south=0.0, east=-120.0, north=80.0)
+        return ds.where(mask, drop=True)
+    else:
+        ds.attrs["bbox"] = BBox(west=-180.0, south=-80.0, east=180.0, north=90.0)
+        return ds
+
+
+GLOBAL_HYCOM = Dataset(
+    name="global_hycom",
+    setup=partial(_create_curvilinear_grid_like_hycom, regional_subset=False),
+    dims=(
+        Dim(
+            name="X",
+            size=4500,
+            chunk_size=4500,
+            data=np.arange(4500),
+            attrs={
+                "standard_name": "projection_x_coordinate",
+                "axis": "X",
+                "point_spacing": "even",
+            },
+        ),
+        Dim(
+            name="Y",
+            size=3298,
+            chunk_size=3298,
+            data=np.arange(3298),
+            attrs={
+                "standard_name": "projection_y_coordinate",
+                "axis": "Y",
+                "point_spacing": "even",
+            },
+        ),
+    ),
+    dtype=np.float64,
+    attrs={
+        "valid_min": 5.0,
+        "valid_max": 15.0,
+        "coordinates": "latitude longitude",
+    },
+)
+
+REGIONAL_HYCOM = Dataset(
+    name="regional_hycom",
+    setup=partial(_create_curvilinear_grid_like_hycom, regional_subset=True),
+    dims=GLOBAL_HYCOM.dims,  # gets subset later
+    dtype=np.float64,
+    attrs={
+        "valid_min": 5.0,
+        "valid_max": 15.0,
+        "coordinates": "latitude longitude",
+    },
+)
+
 POPDS = xr.Dataset(
     {
         "TEMP": (
@@ -1295,6 +1403,8 @@ def create_n320(
     "utm33s_hires": UTM33S_HIRES,
     "utm50s_hires": UTM50S_HIRES,
     "curvilinear": CURVILINEAR,
+    "regional_hycom": REGIONAL_HYCOM,
+    "global_hycom": GLOBAL_HYCOM,
     "hrrr_multiple": HRRR_MULTIPLE,
     "global_nans": GLOBAL_NANS,
     "redgauss_n320": REDGAUSS_N320,

src/xpublish_tiles/testing/tiles.py

@@ -304,6 +304,10 @@
         (Tile(x=442, y=744, z=11), WEBMERC_TMS),
         id="curvilinear_central_us_z11(11/442/744)",
     ),
+    # TODO: uncomment
+    # pytest.param(
+    #     (Tile(x=3, y=5, z=4), WEBMERC_TMS), id="curvilinear_hycom_east_z4(4/3/5)"
+    # ),
 ]
 
 # South America benchmark tiles (for Sentinel dataset)

tests/conftest.py

@@ -10,6 +10,7 @@
     CURVILINEAR,
     EU3035,
     EU3035_HIRES,
+    GLOBAL_HYCOM,
     HRRR,
     REDGAUSS_N320,
     UTM33S,
@@ -119,6 +120,8 @@ def repo(pytestconfig):
 @pytest.fixture(
     params=tuple(map(",".join, product(["-90->90", "90->-90"], ["-180->180", "0->360"])))
     + ("reduced_gaussian_n320",)
+    # TODO: uncomment later
+    # + ("global_hycom",)
 )
 def global_datasets(request):
     param = request.param
@@ -129,6 +132,8 @@ def global_datasets(request):
 
     if param == "reduced_gaussian_n320":
         ds = REDGAUSS_N320.create()
+    elif param in {"global_hycom"}:
+        ds = GLOBAL_HYCOM.create()
     else:
         ds = create_global_dataset(lat_ascending=lat_ascending, lon_0_360=lon_0_360)
     ds.attrs["name"] = param

tests/test_grids.py

@@ -57,6 +57,7 @@
     PARA_HIRES,
     POPDS,
     REDGAUSS_N320,
+    REGIONAL_HYCOM,
     UTM33S_HIRES,
     UTM50S_HIRES,
     Dataset,
@@ -176,6 +177,32 @@
             ),
             id="roms",
         ),
+        pytest.param(
+            REGIONAL_HYCOM.create(),
+            "foo",
+            Curvilinear(
+                crs=CRS.from_user_input(4326),
+                bbox=BBox(
+                    south=0,
+                    north=80,
+                    east=-120,
+                    west=-180,
+                ),
+                X="longitude",
+                Y="latitude",
+                Xdim="X",
+                Ydim="Y",
+                indexes=(
+                    CurvilinearCellIndex(
+                        X=REGIONAL_HYCOM.create().longitude,
+                        Y=REGIONAL_HYCOM.create().latitude,
+                        Xdim="X",
+                        Ydim="Y",
+                    ),
+                ),
+            ),
+            id="hycom",
+        ),
         pytest.param(
             POPDS,
             "UVEL",
@@ -374,9 +401,11 @@ async def test_subset(global_datasets, tile, tms):
         slicer = next(iter(slicers["point"]))
         assert isinstance(slicer, UgridIndexer)
     else:
-        assert isinstance(slicers["latitude"], list)
-        assert isinstance(slicers["longitude"], list)
-        assert len(slicers["latitude"]) == 1  # Y dimension should always have one slice
+        assert isinstance(slicers.get("latitude", slicers.get("Y")), list)
+        assert isinstance(slicers.get("longitude", slicers.get("X")), list)
+        y_slicers = slicers.get("latitude", slicers.get("Y"))
+        assert y_slicers is not None
+        assert len(y_slicers) == 1  # Y dimension should always have one slice
 
     # Check that coordinates are within expected bounds (exact matching with controlled grid)
     actual = await apply_slicers(
@@ -760,6 +789,17 @@ def _create_test_dataset(
             coords={"x": np.arange(array_size), "y": np.arange(array_size)},
         )
         grid = Curvilinear.from_dataset(ds, CRS.from_epsg(4326), "lon", "lat")
+    elif grid_type == "curvilinear_hycom":
+        lon, lat = np.meshgrid(lon_coords, lat_coords)
+        ds = xr.Dataset(
+            {
+                "temp": (["y", "x"], data),
+                "lon": (["y", "x"], lon),
+                "lat": (["y", "x"], lat),
+            },
+            coords={"x": np.arange(array_size), "y": np.arange(array_size)},
+        )
+        grid = Curvilinear.from_dataset(ds, CRS.from_epsg(4326), "lon", "lat")
     elif grid_type == "raster_affine":
         pixel_size_x = (
             (lon_coords[-1] - lon_coords[0]) / (array_size - 1) if array_size > 1 else 1.0
@@ -831,7 +871,9 @@ def test_min_max_zoom_relationship(self, tms_id):
     @pytest.mark.parametrize(
         "tms_id", ["WebMercatorQuad", "WGS1984Quad", "WorldCRS84Quad"]
     )
-    @pytest.mark.parametrize("grid_type", ["rectilinear", "curvilinear", "raster_affine"])
+    @pytest.mark.parametrize(
+        "grid_type", ["rectilinear", "curvilinear", "curvilinear_hycom", "raster_affine"]
+    )
     @given(data=st.data())
     @settings(deadline=None)
     def test_max_zoom_matches_dataset_resolution(self, tms_id, grid_type, data):
@@ -853,7 +895,9 @@ def test_max_zoom_matches_dataset_resolution(self, tms_id, grid_type, data):
     @pytest.mark.parametrize(
         "tms_id", ["WebMercatorQuad", "WGS1984Quad", "WorldCRS84Quad"]
     )
-    @pytest.mark.parametrize("grid_type", ["rectilinear", "curvilinear", "raster_affine"])
+    @pytest.mark.parametrize(
+        "grid_type", ["rectilinear", "curvilinear", "curvilinear_hycom", "raster_affine"]
+    )
     @given(data=st.data())
     @settings(deadline=None)
     def test_min_zoom_matches_renderable_size_limit(self, tms_id, grid_type, data):