Delaunay Triangulation on sphere (unstructured long lat grid)

Question

I work with an unstructured grid of weather data, and I'm trying to plot it. For that I use the Delaunay triangulation.

def triangulate(vertices, x="Longitude", y="Latitude"):
    """
    Generate a triangular mesh for the given x,y,z vertices, using Delaunay triangulation.
    For large n, typically results in about double the number of triangles as vertices.
    """
    triang = Delaunay(vertices[[x, y]].values)
    return pd.DataFrame(triang.simplices, columns=['v0', 'v1', 'v2'])

x = ((np.rad2deg(data.clon) - 180) % 360) - 180
y = np.rad2deg(data.clat)
    
value = data[attr][t][i]
pts = np.stack((x, y, value)).T
verts = pd.DataFrame(pts, columns=['Longitude', 'Latitude', 'value'])
triangulate(verts)

I use geoviews to plot the trimesh in a NorthPolarStereo projection. There I have the issue that the triangles don't overlap due to the periodic behavior of the sin/cos.

Is there a way to modify the Delaunay triangulation that the triangles overlap at the 0° / 360° area?

The plotting is done by the following:

trimesh = rasterize(gv.TriMesh((g.tris, hv.Points(g.verts, vdims='value')), crs=ccrs.PlateCarree()))
    
(trimesh).opts(
                cmap='Greys',
                fig_size=180,
                aspect=1,
                projection=ccrs.NorthPolarStereo())

Thanks for any advise. In case more code is needed let me know :)

Answer 1

Some libraries implement direct support for a triangulation over the surface of a sphere. But absent that, one work-around solution is to create a planar domain that extends larger than 180w and 180e, say from -360 and +360. You could include multiple vertices for your set of data points. A point located at longitude -90 would also be added as a vertex at longitude 270. You could then conduct whatever processing you need to do over the entire Delaunay triangulation and clip the results.

I will attach a picture from my employer's Dispatch Weather Client (DWC) product that shows this approach for surface temperature data from weather-reporting stations. This example extends the domain across longitudes. For your purposes, you could use the same idea extending data beyond the polar latitudes.

Updated Answer 21 Oct 2022

I am not familiar with the API you are using. So, reading your code, I am not completely sure what you're doing. But based on your picture, it looks like you might be constructing your Delaunay mesh from latitude/longitude coordinates and then applying the Polar Stereographic projection over the results. If that is a correct assessment, then I think you will get better results if you perform the map projection (from lat/lon to Cartesian x/y) first and then build your mesh. This approach will automatically get rid of the seam that you are seeing in your plot. And, better yet, it will convert the non-isotropic latitude/longitude coordinates to a more conformal Cartesian coordinate system. Results from the Delaunay tend to be more accurate when its input coordinates are uniform in their representation of distance. Here's a second picture. Because the lat/lon coordinates are projected to the Polar Stereographic, the second mesh I produced doesn't even need the redundant vertices.

Here's a bit of code for the Polar Stereographic, spherical case:

public void transform(
    boolean southPole, double centralMeridian,
    double latitude, double longitude,
    double[] output) {

double earthRadius = 6378137.0;  // WGS84 semi-major axis
double phi = Math.toRadians(latitude);
double delta = Math.toRadians(longitude - centralMeridian);
if (southPole) {
    delta = -delta;
    phi = - phi;
}
double p = 2 * tan(PI / 4 - phi / 2) * earthRadius;
output[0] = p * sin(delta);   // X coordinate
output[1] = -p * cos(delta);  // Y coordinate

}