Spherical triangulations, tessellation

[DanielVandH]

• Proper testing of the triangulations/tessellations
• Test peek so that NaturalNeighbours might just work directly (?)
• Test refinement
• Doc examples
• Figure out better plotting methods for the doc examples especially... can't depend on Makie here (because it would be a circular dependency) so will have to paste some code into the examples.
• News

This PR introduces spherical triangulations and tessellations. Ghost vertices are interpreted as the north pole. Pretty experimental for now.

[codecov]

Codecov Report

Attention: Patch coverage is 58.10277% with 212 lines in your changes missing coverage. Please review.

Project coverage is 93.07%. Comparing base (ec3059e) to head (0621696).

Files with missing lines	Patch %	Lines
src/spherical/plotting.jl	17.27%	91 Missing ⚠️
src/spherical/triangulation.jl	58.41%	42 Missing ⚠️
src/spherical/refine.jl	0.00%	37 Missing ⚠️
src/spherical/voronoi.jl	0.00%	25 Missing ⚠️
src/spherical/utils.jl	44.00%	14 Missing ⚠️
src/spherical/geometry.jl	96.61%	2 Missing ⚠️
src/spherical/spherical_triangles.jl	98.46%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #187      +/-   ##
==========================================
- Coverage   94.82%   93.07%   -1.76%     
==========================================
  Files         103      111       +8     
  Lines       10364    10805     +441     
==========================================
+ Hits         9828    10057     +229     
- Misses        536      748     +212     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[DanielVandH]

@asinghvi17 This is a working implementation of spherical delaunay + voronoi. It also supports mesh refinement (the points inserted are the centroids [the ones that split the spherical triangles into equal-area regions, not the median] instead of circumcircles, seems to be more visually appealing). I think with this implementation the peek keyword argument might just work, so NaturalNeighbours.jl might be able to just use all of this directly (though it might need some adjustments to use spherical areas/distances, which I do have implemented in this PR, NaturalNeighbours.jl might just need some overloads for ::SphericalTriangulation).

Screen.Recording.2024-09-22.132518.mp4

Maybe you would like to try this branch out at some point to see what you think? The code is in theDelaunayTriangulation.SphericalDelaunay module. There are some examples at the bottom of the spherical.jl folder in the tests. I have some Makie overloads in there too to support plotting. The mesh recipes render terribly though.

[DanielVandH]

Some thoughts (for this PR and also for what needs to be done in NaturalNeighbours)

• find_triangle should, by default, return the triangle in projected space. This means that the result will not return the actual spherical triangle containing the point for points outside of the convex hull. find_spherical_triangle should be the replacement in this case. So I probably need to revert parts of the last commits. The fix for RepresentativeCoordinates should stay though, storing those as NaNs is bad and leads to errors like what was reported on Slack.
• Add tests that check that find_triangle returns the wrong triangle only for points outside of the convex hull. Also add warnings to find_triangle that the result may be wrong for spherical triangulations and to use find_spherical_triangle instead.

For NaturalNeighbours:

• Check that areas are being computed for spherical triangles and not their projection (spherical_triangle_area), which triangle_area(tri, T) would use. I think right now it won't be, since the areas are probably being computed using triangle_area(p, q, r) which knows nothing about the spherical points depending on if the points were already projected.
• Check that distances are being computed using spherical_distance? Maybe need some methods taking the triangulation as the first argument to allow dispatching on SphericalTriangulations.
• Check that the Voronoi cells involving the ghost vertex know to use the north pole. See spherical_voronoi for how this is being done already.
• For data outside of the convex hull, NaturalNeighbours.jl will project onto the solid edge of the ghost triangle containing the point, and the interpolant just averages the data at the two vertices of that edge. For the spherical case, maybe I need to demand that the data be known at the north pole and, for SphericalTriangulations, ensure that the natural neighbour coordinates are being used and the interpolant isn't doing extrapolation.

[asinghvi17]

That sounds good. I got triangulation and Voronoi tessellation running for my test cases, but have not yet tried NaturalNeighbours. Will hopefully get to that today.

Out of curiosity, have you considered implementing the GeoInterface for Voronoi tessellations (maybe as a feature collection)? It was a bit hard for me to parse how to get a list of points per polygon out, and having them be GeoInterface polygons would allow things like Makie plotting, writing to shapefile, etc. to "just work".

[DanielVandH]

That sounds good. I got triangulation and Voronoi tessellation running for my test cases, but have not yet tried NaturalNeighbours. Will hopefully get to that today.

Great!

Out of curiosity, have you considered implementing the GeoInterface for Voronoi tessellations (maybe as a feature collection)? It was a bit hard for me to parse how to get a list of points per polygon out, and having them be GeoInterface polygons would allow things like Makie plotting, writing to shapefile, etc. to "just work".

I've never used any of the GeoInterface stuff, so I've not ever looked into it. Is it something that can exist as a package extension? I don't really want to enforce any sort of interface that requires a dependency by default.

Also regarding getting the coordinates: There are a few different ways. What did you end up using? Maybe I need to improve the documentation here. VoronoiTessellations see a lot less love from me than Delaunay triangulations unfortunately, since I rarely actually need to use voronoi.

[asinghvi17]

I've never used any of the GeoInterface stuff, so I've not ever looked into it. Is it something that can exist as a package extension? I don't really want to enforce any sort of interface that requires a dependency by default.

I think GeoInterface could probably be implemented in a package extension, but there is an argument to be made that you can use it as a fallback for methods like getxy. The main point of GeoInterface is that it's an interface that defines (a) how to access the contents of a geometry and (b) how to convert geometries between packages. Both of these could be useful here in different ways.

Also regarding getting the coordinates: There are a few different ways. What did you end up using?

vpolys = [
	[get_polygon_point(vorn, i) for i in DelTri.get_polygon(vorn, j)] 
	for j in Iterators.filter(!in(DelaunayTriangulation.get_unbounded_polygons(vorn)), DelTri.each_polygon_index(vorn))
]

(seems pretty involved but I couldn't see a better way to do it)

[DanielVandH]

I think GeoInterface could probably be implemented in a package extension, but there is an argument to be made that you can use it as a fallback for methods like getxy. The main point of GeoInterface is that it's an interface that defines (a) how to access the contents of a geometry and (b) how to convert geometries between packages. Both of these could be useful here in different ways.

I get that it can be useful for something like that, but another main concern for me is that, if I make it the default fallback, it's a bit of extra maintenance on my part since I don't really ever use any of that stuff myself, if that makes sense? I'd be happy to maintain it as an extra extension but maybe not as a fallback.

--

For the coordinates.. I guess that is probably not a bad way.

[get_polygon_point(vorn, i) for i in DelTri.get_polygon(vorn, j)]

is the same as get_polygon_coordinates(vorn, i) for i in each_polygon_index(vorn), except that it won't work for unbounded polygons (as you have seen). I usually need both the finite and infinite polygons when I've extracted the polygons so I've never implemented a better method for filtering out just the finite ones - is there a particular reason to want to discard all the infinite ones instead of just clipping them?

Maybe a method each_bounded_polygon(vorn) would be useful that basically just does that filter you're doing. Usually I'm working with the vertices rather than vectors of coordinates so this stuff is probably underdeveloped.

[DanielVandH]

I'm still yet to come back to this but just curious if you got around to trying it with NaturalNeighbours @asinghvi17? No problem if not but if you have feel free to let comment any other issues I need to fix.