From 26856bfc139f6bf3e891591f420c77f5faac61fc Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 27 Apr 2026 10:57:31 -0400 Subject: [PATCH 1/2] [#2521] Correct PIC check for curvilinear grids with mesh.type="spherical" When the mesh type is "spherical" for a curvilinear structured grid, we now use a PIC check that maps the curvilinear cell to x,y,z space on the unit sphere. Within the cell, we construct an orthonormal basis for the cell using vectors that point between edge centers in the xsi and eta directions. Particle positions are converted to x,y,z space and projected onto the cell plane. This approach is identical to what is done for unstructure "spherical" grids; it is robust across the polar caps and across the antimeridians. This change removes the need for consistency checking in the _xinterpolator.py for curvilinear grids for points in the vicinity of the antimeridian --- src/parcels/_core/index_search.py | 165 +++++++++++++++---- src/parcels/interpolators/_xinterpolators.py | 8 - 2 files changed, 130 insertions(+), 43 deletions(-) diff --git a/src/parcels/_core/index_search.py b/src/parcels/_core/index_search.py index 5e2c223d0..5e84eee2e 100644 --- a/src/parcels/_core/index_search.py +++ b/src/parcels/_core/index_search.py @@ -92,56 +92,151 @@ def _search_time_index(field: Field, time: np.ndarray): def curvilinear_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray, xi: np.ndarray): - xsi = eta = -1.0 * np.ones(len(x), dtype=float) - invA = np.array( - [ - [1, 0, 0, 0], - [-1, 1, 0, 0], - [-1, 0, 0, 1], - [1, -1, 1, -1], - ] - ) + """Locate points within 2D curvilinear grid cells. - px = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]]) + For spherical grids, the bilinear inverse runs in a tangent plane on a sphere + through each cell's centroid — robust across the antimeridian and the + polar cap. For flat grids it runs directly in (lon, lat). In both + cases the returned (xsi, eta) are bilinear weights in the space the inverse + was solved in; downstream interpolators must agree on that same space. + """ + clon = np.asarray( + [grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]], + dtype=float, + ) + clat = np.asarray( + [grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]], + dtype=float, + ) if grid._mesh == "spherical": - # Map grid and particle longitudes to [-180,180) - px = ((px + 180.0) % 360.0) - 180.0 - x = ((x + 180.0) % 360.0) - 180.0 + xsi, eta = _bilinear_inverse_tangent_plane(clon, clat, x, y) + is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0) + else: + xsi, eta = _bilinear_inverse_latlon(clon, clat, x, y, apply_antimeridian_shift=False) + is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0) - # Create a mask for antimeridian cells - lon_span = px.max(axis=0) - px.min(axis=0) - antimeridian_cell = lon_span > 180.0 + return is_in_cell, np.column_stack((xsi, eta)) - if np.any(antimeridian_cell): - # For any antimeridian cell ... - # If particle longitude is closer to 180.0, then negative cell longitudes need to be bumped by +360 - mask = (px < 0.0) & antimeridian_cell[np.newaxis, :] & (x[np.newaxis, :] > 0.0) - px[mask] += 360.0 - # If particle longitude is closer to -180.0, then positive cell longitudes need to be bumped by -360 - mask = (px > 0.0) & antimeridian_cell[np.newaxis, :] & (x[np.newaxis, :] < 0.0) - px[mask] -= 360.0 - py = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]]) +_invA = np.array( + [ + [1, 0, 0, 0], + [-1, 1, 0, 0], + [-1, 0, 0, 1], + [1, -1, 1, -1], + ] +) + - a, b = np.dot(invA, px), np.dot(invA, py) +def _bilinear_inverse(px, py, xq, yq): + """Solve for (xsi, eta) s.t. bilinear blend of (px, py) corners equals (xq, yq).""" + eta_init = -1.0 * np.ones(len(xq), dtype=float) + a, b = np.dot(_invA, px), np.dot(_invA, py) aa = a[3] * b[2] - a[2] * b[3] - bb = a[3] * b[0] - a[0] * b[3] + a[1] * b[2] - a[2] * b[1] + x * b[3] - y * a[3] - cc = a[1] * b[0] - a[0] * b[1] + x * b[1] - y * a[1] + bb = a[3] * b[0] - a[0] * b[3] + a[1] * b[2] - a[2] * b[1] + xq * b[3] - yq * a[3] + cc = a[1] * b[0] - a[0] * b[1] + xq * b[1] - yq * a[1] det2 = bb * bb - 4 * aa * cc with np.errstate(divide="ignore", invalid="ignore"): - det = np.where(det2 > 0, np.sqrt(det2), eta) - eta = np.where(abs(aa) < 1e-12, -cc / bb, np.where(det2 > 0, (-bb + det) / (2 * aa), eta)) - + det = np.where(det2 > 0, np.sqrt(det2), eta_init) + eta = np.where(abs(aa) < 1e-12, -cc / bb, np.where(det2 > 0, (-bb + det) / (2 * aa), eta_init)) xsi = np.where( abs(a[1] + a[3] * eta) < 1e-12, - ((y - py[0]) / (py[1] - py[0]) + (y - py[3]) / (py[2] - py[3])) * 0.5, - (x - a[0] - a[2] * eta) / (a[1] + a[3] * eta), + ((yq - py[0]) / (py[1] - py[0]) + (yq - py[3]) / (py[2] - py[3])) * 0.5, + (xq - a[0] - a[2] * eta) / (a[1] + a[3] * eta), ) - is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0) + return xsi, eta - return is_in_cell, np.column_stack((xsi, eta)) + +def _bilinear_inverse_latlon(clon, clat, x, y, apply_antimeridian_shift=True): + """Bilinear inverse in raw lat/lon, with antimeridian corner shifting.""" + px = clon.copy() + xq = np.asarray(x, dtype=float).copy() + if apply_antimeridian_shift: + px = ((px + 180.0) % 360.0) - 180.0 + xq = ((xq + 180.0) % 360.0) - 180.0 + lon_span = px.max(axis=0) - px.min(axis=0) + antimeridian_cell = lon_span > 180.0 + if np.any(antimeridian_cell): + mask = (px < 0.0) & antimeridian_cell[np.newaxis, :] & (xq[np.newaxis, :] > 0.0) + px[mask] += 360.0 + mask = (px > 0.0) & antimeridian_cell[np.newaxis, :] & (xq[np.newaxis, :] < 0.0) + px[mask] -= 360.0 + py = clat + yq = np.asarray(y, dtype=float) + return _bilinear_inverse(px, py, xq, yq) + + +def _bilinear_inverse_tangent_plane(clon, clat, x, y): + """Bilinear inverse after tangent-plane projection through the cell centroid. + + Avoids lat/lon singularities at poles and on the antimeridian. + """ + px, py, xq, yq = _spherical_project_cell_and_query(clon, clat, x, y) + return _bilinear_inverse(px, py, xq, yq) + + +def _spherical_project_cell_and_query(clon, clat, x, y): + """Project 4 cell corners and the query onto a plane defined by the cell. + + The plane's basis (e_u, e_v) is built from the cell's edge-midpoint + difference vectors, then orthonormalized using Gramm-Schmidt + + Inputs + ------ + clon, clat : (4, N) arrays of corner lon/lat in degrees + x, y : (N,) query lon/lat in degrees + + Returns + ------- + px, py : (4, N) projected corner coordinates (u, v) + xq, yq : (N,) projected query coordinates in the same (u, v) + + Notes + ----- + Orthogonal projection is 2-to-1 on the full sphere: a query and its + antipode map to the same (u, v). In practice this is harmless because the + downstream bilinear gate ``(xsi, eta) in [0, 1]`` only admits queries whose + projected magnitude is ≲ the cell's angular size. + """ + cX, cY, cZ = _latlon_rad_to_xyz(np.deg2rad(clat), np.deg2rad(clon)) # (4, N) each + qX, qY, qZ = _latlon_rad_to_xyz(np.deg2rad(np.asarray(y, dtype=float)), np.deg2rad(np.asarray(x, dtype=float))) + + # Cell-intrinsic tangent basis from edge-midpoint difference vectors. + # Corner order is CCW: c0=(yi,xi), c1=(yi,xi+1), c2=(yi+1,xi+1), c3=(yi+1,xi). + # u-axis (along xsi): right-edge midpoint - left-edge midpoint = (c1+c2) - (c0+c3) + # v-axis (along eta): top-edge midpoint - bottom-edge midpoint = (c2+c3) - (c0+c1) + ux = (cX[1] + cX[2]) - (cX[0] + cX[3]) + uy = (cY[1] + cY[2]) - (cY[0] + cY[3]) + uz = (cZ[1] + cZ[2]) - (cZ[0] + cZ[3]) + u_norm = np.sqrt(ux * ux + uy * uy + uz * uz) + u_norm = np.where(u_norm == 0.0, 1.0, u_norm) + e_ux, e_uy, e_uz = ux / u_norm, uy / u_norm, uz / u_norm + + vx = (cX[2] + cX[3]) - (cX[0] + cX[1]) + vy = (cY[2] + cY[3]) - (cY[0] + cY[1]) + vz = (cZ[2] + cZ[3]) - (cZ[0] + cZ[1]) + # Gram-Schmidt: make e_v perpendicular to e_u so the bilinear inverse + # in (u, v) is well-conditioned even on skewed cells. + v_dot_eu = vx * e_ux + vy * e_uy + vz * e_uz + vx = vx - v_dot_eu * e_ux + vy = vy - v_dot_eu * e_uy + vz = vz - v_dot_eu * e_uz + v_norm = np.sqrt(vx * vx + vy * vy + vz * vz) + v_norm = np.where(v_norm == 0.0, 1.0, v_norm) + e_vx, e_vy, e_vz = vx / v_norm, vy / v_norm, vz / v_norm + + # Orthogonal projection onto span(e_u, e_v): two dot products per point. + # Same pattern as uxgrid_point_in_cell (drop the component normal to the plane). + def _project(vx, vy, vz): + u = vx * e_ux + vy * e_uy + vz * e_uz + v = vx * e_vx + vy * e_vy + vz * e_vz + return u, v + + px_u, px_v = _project(cX, cY, cZ) # (4, N) each + xq_u, xq_v = _project(qX, qY, qZ) # (N,) each + return px_u, px_v, xq_u, xq_v def _search_indices_curvilinear_2d( diff --git a/src/parcels/interpolators/_xinterpolators.py b/src/parcels/interpolators/_xinterpolators.py index 7c0eebdcb..725422d02 100644 --- a/src/parcels/interpolators/_xinterpolators.py +++ b/src/parcels/interpolators/_xinterpolators.py @@ -177,7 +177,6 @@ def CGrid_Velocity( yi, eta = grid_positions["Y"]["index"], grid_positions["Y"]["bcoord"] zi, zeta = grid_positions["Z"]["index"], grid_positions["Z"]["bcoord"] ti, tau = grid_positions["T"]["index"], grid_positions["T"]["bcoord"] - lon = particle_positions["lon"] U = vectorfield.U.data V = vectorfield.V.data @@ -305,13 +304,6 @@ def _compute_corner_data(data, selection_dict) -> np.ndarray: u /= conversion v /= conversion - # check whether the grid conversion has been applied correctly - xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3] - dlon = xx - lon - if grid._mesh == "spherical": - dlon = ((dlon + 180.0) % 360.0) - 180.0 - u = np.where(np.abs(dlon / lon) > 1e-4, np.nan, u) - if vectorfield.W: W = vectorfield.W.data From f59efdbceaae4e2c37a4849cc8db8b190b8c7a38 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 27 Apr 2026 11:22:41 -0400 Subject: [PATCH 2/2] [#2521] Verify (xsi, eta) reconstruct query position via 3D corner blend After the curvilinear PIC moved to a tangent-plane bilinear inverse, the old "lat/lon-bilinear blend of corner lons reconstructs query lon" check no longer holds. Replace it with a blend of corner 3D unit-sphere vectors, renormalized onto the sphere, and compare the recovered lat/lon to the query. Tolerance 1e-5 gives ~30x headroom over the worst observed error (~3e-7) on the 1/4 degree NEMO test grid. --- tests/test_index_search.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/tests/test_index_search.py b/tests/test_index_search.py index bb7ec3f3b..b85e6e3b6 100644 --- a/tests/test_index_search.py +++ b/tests/test_index_search.py @@ -4,7 +4,7 @@ import parcels.tutorial from parcels import Field, XGrid -from parcels._core.index_search import _search_indices_curvilinear_2d +from parcels._core.index_search import _latlon_rad_to_xyz, _search_indices_curvilinear_2d from parcels._datasets.structured.generic import datasets from parcels.interpolators import XLinear @@ -73,6 +73,25 @@ def test_indexing_nemo_curvilinear(): for i in range(lons.shape[0]): np.testing.assert_allclose(px[1, i], px[:, i], atol=5) - # Reconstruct lons values from cornerpoints - xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3] - np.testing.assert_allclose(xx, lons, atol=1e-6) + # Each query should have been located inside some cell: xsi, eta in [0, 1] + assert np.all((xsi >= 0) & (xsi <= 1)), f"xsi out of [0,1]: {xsi}" + assert np.all((eta >= 0) & (eta <= 1)), f"eta out of [0,1]: {eta}" + + # Reconstruct query lat/lon by bilinear-blending the 4 corner 3D unit-sphere + # vectors with (xsi, eta) and renormalizing onto the unit sphere. Tolerance + # reflects the residual spherical curvature for a 1/4° NEMO cell. + clat = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]]) + clon = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]]) + cX, cY, cZ = _latlon_rad_to_xyz(np.deg2rad(clat), np.deg2rad(clon)) + w = np.array([(1 - xsi) * (1 - eta), xsi * (1 - eta), xsi * eta, (1 - xsi) * eta]) + x = np.sum(w * cX, axis=0) + y = np.sum(w * cY, axis=0) + z = np.sum(w * cZ, axis=0) + n = np.sqrt(x * x + y * y + z * z) + x, y, z = x / n, y / n, z / n + lat_recon = np.rad2deg(np.arcsin(z)) + lon_recon = np.rad2deg(np.arctan2(y, x)) + lons_wrapped = ((lons + 180) % 360) - 180 + dlon = ((lon_recon - lons_wrapped + 180) % 360) - 180 + np.testing.assert_allclose(lat_recon, lats, atol=1e-5) + np.testing.assert_allclose(dlon, 0.0, atol=1e-5)