Example 1 - Triangulation of arbitrary points#
stripy provides a python interfact to TRIPACK and SRFPACK (Renka 1997c,d) as a triangulation class that would typically be used as follows:
import stripy
triangulation = stripy.Triangulation(x=vertices_x, y=vertices_y)
areas = triangulation.areas()
The methods of the Triangulation class include interpolation, smoothing and gradients (from SRFPACK), triangle areas, point location by simplex and nearest vertex, refinement operations by edge or centroid, and neighbourhood search / distance computations through a k-d tree algorithm suited to points on the surface of a unit sphere. stripy also includes template triangulated meshes with refinement operations.
In this notebook we introduce the Triangulation class itself.
Notebook contents#
References#
Renka, R. J. (1997), Algorithm 772: STRIPACK: Delaunay triangulation and Voronoi diagram on the surface of a sphere, ACM Transactions on Mathematical Software (TOMS).
Renka, R. J. (1997), Algorithm 773: SSRFPACK: interpolation of scattered data on the surface of a sphere with a surface under tension, ACM Transactions on Mathematical Software (TOMS), 23(3), 435–442, doi:10.1145/275323.275330.
Renka, R. J. (1996), Algorithm 751; TRIPACK: a constrained two-dimensional Delaunay triangulation package, ACM Transactions on Mathematical Software, 22(1), 1–8, doi:10.1145/225545.225546.
Renka, R. J. (1996), Algorithm 752; SRFPACK: software for scattered data fitting with a constrained surface under tension, ACM Transactions on Mathematical Software, 22(1), 9–17, doi:10.1145/225545.225547.
The next example is Ex2-SphericalGrids
Triangulate the vertices of a square mesh#
import stripy as stripy
import numpy as np
# Vertices of a square mesh
vertices_xy = np.array([[0.0, 0.0],
[1.0, 0.0],
[0.0, 1.0],
[1.0, 1.0]])
vertices_x = vertices_xy.T[0]
vertices_y = vertices_xy.T[1]
triangulation = stripy.Triangulation(x=vertices_x, y=vertices_y, permute=False)
print(triangulation.areas())
print(triangulation.npoints)
[0.5 0.5]
4
This creates a triangulation object constructed using the wrapped fortran code of Renka (1997). The triangulation object has a number of useful methods and attached data which can be listed with
help(triangulation)
refined_triangulation = stripy.Triangulation(x=vertices_x, y=vertices_y, refinement_levels=4, permute=False)
print(refined_triangulation.npoints)
289
Making a plot of the triangulation#
We can make a plot of the two grids and the most straightforward way to display the information is through a standard map projection of the sphere to the plane.
(Here we superimpose the points on a global map of coastlines using the cartopy map library and use the Mollweide projection.
Other projections to try include Robinson, Orthographic, PlateCarree)
%matplotlib inline
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20, 10), facecolor="none", sharey=True)
## Plot the vertices and the edges for the original mesh
x = triangulation.x
y = triangulation.y
simplices = triangulation.simplices
ax1.triplot(x, y, simplices, linewidth=0.5, color='black')
ax1.scatter(x, y, color='Red', alpha=0.5, marker='o')
## Plot the vertices and the edges for the refined mesh
x = refined_triangulation.x
y = refined_triangulation.y
simplices = refined_triangulation.simplices
ax2.triplot(x, y, simplices, linewidth=0.5, color='black')
ax2.scatter(x, y, color='Red', alpha=0.5, marker='o')
plt.show()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[4], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
3 import matplotlib.pyplot as plt
6 fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20, 10), facecolor="none", sharey=True)
File ~/.local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:2482, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
2480 kwargs['local_ns'] = self.get_local_scope(stack_depth)
2481 with self.builtin_trap:
-> 2482 result = fn(*args, **kwargs)
2484 # The code below prevents the output from being displayed
2485 # when using magics with decorator @output_can_be_silenced
2486 # when the last Python token in the expression is a ';'.
2487 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File ~/.local/lib/python3.12/site-packages/IPython/core/magics/pylab.py:103, in PylabMagics.matplotlib(self, line)
98 print(
99 "Available matplotlib backends: %s"
100 % _list_matplotlib_backends_and_gui_loops()
101 )
102 else:
--> 103 gui, backend = self.shell.enable_matplotlib(args.gui)
104 self._show_matplotlib_backend(args.gui, backend)
File ~/.local/lib/python3.12/site-packages/IPython/core/interactiveshell.py:3667, in InteractiveShell.enable_matplotlib(self, gui)
3664 import matplotlib_inline.backend_inline
3666 from IPython.core import pylabtools as pt
-> 3667 gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
3669 if gui != None:
3670 # If we have our first gui selection, store it
3671 if self.pylab_gui_select is None:
File ~/.local/lib/python3.12/site-packages/IPython/core/pylabtools.py:338, in find_gui_and_backend(gui, gui_select)
321 def find_gui_and_backend(gui=None, gui_select=None):
322 """Given a gui string return the gui and mpl backend.
323
324 Parameters
(...)
335 'WXAgg','Qt4Agg','module://matplotlib_inline.backend_inline','agg').
336 """
--> 338 import matplotlib
340 if _matplotlib_manages_backends():
341 backend_registry = matplotlib.backends.registry.backend_registry
ModuleNotFoundError: No module named 'matplotlib'
Predefined meshes#
One common use of stripy is in meshing x,y coordinates and, to this end, we provide pre-defined meshes for square and elliptical triangulations. A random mesh is included as a counterpoint to the regular meshes. The square mesh defined above can be created directly using:
triangulation = stripy.cartesian_meshes.square_mesh(extent, spacingX, spacingY, refinement_levels=0)
refined_triangulation = stripy.cartesian_meshes.square_mesh(extent, spacingX, spacingY, refinement_levels=3)
This capability is shown in a companion notebook Ex2-CartesianGrids