Example 3 - stripy interpolation#
SRFPACK is a Fortran 77 software package that constructs a smooth interpolatory or approximating surface to data values associated with arbitrarily distributed points. It employs automatically selected tension factors to preserve shape properties of the data and avoid overshoot and undershoot associated with steep gradients.
The next three examples demonstrate the interface to SRFPACK provided through stripy
Notebook contents#
The next example is Ex4-Gradients
Define two different meshes#
Create a fine and a coarse mesh without common points
import stripy as stripy
xmin = 0.0
xmax = 10.0
ymin = 0.0
ymax = 10.0
extent = [xmin, xmax, ymin, ymax]
spacingX = 1.0
spacingY = 1.0
cmesh = stripy.cartesian_meshes.elliptical_mesh(extent, spacingX, spacingY, refinement_levels=1)
fmesh = stripy.cartesian_meshes.elliptical_mesh(extent, spacingX, spacingY, refinement_levels=3)
print("coarse mesh points = {}".format(cmesh.npoints))
print("fine mesh points = {}".format(fmesh.npoints))
coarse mesh points = 53
fine mesh points = 437
help(cmesh.interpolate)
Help on method interpolate in module stripy.cartesian:
interpolate(xi, yi, zdata, order=1, grad=None, sigma=None, threads=1) method of stripy.cartesian_meshes.elliptical_mesh instance
Base class to handle nearest neighbour, linear, and cubic interpolation.
Given a triangulation of a set of nodes and values at the nodes,
this method interpolates the value at the given xi,yi coordinates.
Args:
xi : float / array of floats, shape (l,)
x Cartesian coordinate(s)
yi : float / array of floats, shape (l,)
y Cartesian coordinate(s)
zdata : array of floats, shape (n,)
value at each point in the triangulation
must be the same size of the mesh
order : int (default=1)
order of the interpolatory function used:
- `order=0` = nearest-neighbour
- `order=1` = linear
- `order=3` = cubic
sigma : array of floats, shape (6n-12)
precomputed array of spline tension factors from
`get_spline_tension_factors(zdata, tol=1e-3, grad=None)`
(only used in cubic interpolation)
threads : int or 'auto', optional; default : 1
Number of threads to use for interpolation.
By default, only a single thread will be used. Use
`threads='auto'` to attempt to automatically determine how
many threads to use based on the size of the input and output
data.
Negative values count backwards, such that -1 is equivalent to
`multiprocessing.cpu_count()`, -2 to `cpu_count() - 1`, etc.
Returns:
zi : float / array of floats, shape (l,)
interpolates value(s) at (xi, yi)
err : int / array of ints, shape (l,)
whether interpolation (0), extrapolation (1) or error (other)
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
def mesh_fig(mesh, meshR, name):
fig = plt.figure(figsize=(10, 10), facecolor="none")
ax = plt.subplot(111)
ax.axis('off')
generator = mesh
refined = meshR
x0 = generator.x
y0 = generator.y
xR = refined.x
yR = refined.y
ax.scatter(x0, y0, color="Red", marker="o", s=150.0)
ax.scatter(xR, yR, color="DarkBlue", marker="o", s=50.0)
ax.triplot(xR, yR, refined.simplices, color="black", linewidth=0.5)
fig.savefig(name, dpi=250, transparent=True)
return
mesh_fig(cmesh, fmesh, "Two grids" )
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
3 import matplotlib.pyplot as plt
4 import numpy as np
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'
Analytic function#
Define a relatively smooth function that we can interpolate from the coarse mesh to the fine mesh and analyse
def analytic(xs, ys, k1, k2):
return np.cos(k1*xs) * np.sin(k2*ys)
coarse_afn = analytic(cmesh.x, cmesh.y, 0.1, 1.0)
fine_afn = analytic(fmesh.x, fmesh.y, 0.1, 1.0)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 4
1 def analytic(xs, ys, k1, k2):
2 return np.cos(k1*xs) * np.sin(k2*ys)
----> 4 coarse_afn = analytic(cmesh.x, cmesh.y, 0.1, 1.0)
5 fine_afn = analytic(fmesh.x, fmesh.y, 0.1, 1.0)
Cell In[4], line 2, in analytic(xs, ys, k1, k2)
1 def analytic(xs, ys, k1, k2):
----> 2 return np.cos(k1*xs) * np.sin(k2*ys)
NameError: name 'np' is not defined
The analytic function on the different samplings#
It is helpful to be able to view a mesh to verify that it is an appropriate choice. Here, for example, we visualise the analytic function on the elliptical mesh.
def mesh_field_fig(mesh, field, name):
fig = plt.figure(figsize=(10, 10), facecolor="none")
ax = plt.subplot(111)
ax.axis('off')
ax.tripcolor(mesh.x, mesh.y, mesh.simplices, field)
fig.savefig(name, dpi=250, transparent=True)
return
mesh_field_fig(cmesh, coarse_afn, "coarse analytic")
mesh_field_fig(fmesh, fine_afn, "fine analytic")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 13
9 fig.savefig(name, dpi=250, transparent=True)
11 return
---> 13 mesh_field_fig(cmesh, coarse_afn, "coarse analytic")
14 mesh_field_fig(fmesh, fine_afn, "fine analytic")
NameError: name 'coarse_afn' is not defined
Interpolation from coarse to fine#
The interpolate method of the Triangulation takes arrays of x, y points and an array of data on the mesh vertices.
It returns an array of interpolated values and a status array that states whether each value represents an interpolation, extrapolation or neither (an error condition).
The interpolation can be:
nearest-neighbour (order=0)
linear (order=1)
cubic spline (order=3)
interp_c2f1, err = cmesh.interpolate(fmesh.x, fmesh.y, order=1, zdata=coarse_afn)
interp_c2f3, err = cmesh.interpolate(fmesh.x, fmesh.y, order=3, zdata=coarse_afn)
err_c2f1 = interp_c2f1-fine_afn
err_c2f3 = interp_c2f3-fine_afn
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 1
----> 1 interp_c2f1, err = cmesh.interpolate(fmesh.x, fmesh.y, order=1, zdata=coarse_afn)
2 interp_c2f3, err = cmesh.interpolate(fmesh.x, fmesh.y, order=3, zdata=coarse_afn)
4 err_c2f1 = interp_c2f1-fine_afn
NameError: name 'coarse_afn' is not defined
def axis_mesh_field(ax, mesh, field, label):
ax.axis('off')
x0 = mesh.x
y0 = mesh.y
im = ax.tripcolor(x0, y0, mesh.simplices, field)
ax.set_title(str(label))
fig.colorbar(im, ax=ax)
return
fig, ax = plt.subplots(2,2, figsize=(10,8))
axis_mesh_field(ax[0,0], fmesh, interp_c2f1, "interp1")
axis_mesh_field(ax[0,1], fmesh, interp_c2f3, "interp3")
axis_mesh_field(ax[1,0], fmesh, err_c2f1, "interp_err1")
axis_mesh_field(ax[1,1], fmesh, err_c2f3, "interp_err3")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 14
10 fig.colorbar(im, ax=ax)
11 return
---> 14 fig, ax = plt.subplots(2,2, figsize=(10,8))
16 axis_mesh_field(ax[0,0], fmesh, interp_c2f1, "interp1")
17 axis_mesh_field(ax[0,1], fmesh, interp_c2f3, "interp3")
NameError: name 'plt' is not defined
Interpolate to grid#
Interpolating to a grid is useful for exporting maps of a region. The interpolate_to_grid method interpolates mesh data to a regular grid defined by the user. Values outside the convex hull are extrapolated.
interpolate_to_gridis a convenience function that yields identical results to interpolating over a meshed grid using theinterpolatemethod.
resX = resY = 100
grid_x = np.linspace(xmin, xmax, resX)
grid_y = np.linspace(ymin, ymax, resY)
grid_z1 = fmesh.interpolate_to_grid(grid_x, grid_y, interp_c2f3)
# compare with `interpolate` method
grid_xcoords, grid_ycoords = np.meshgrid(grid_x, grid_y)
grid_z2, ierr = fmesh.interpolate(grid_xcoords.ravel(), grid_ycoords.ravel(), interp_c2f3, order=3)
grid_z2 = grid_z2.reshape(resY,resX)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 3
1 resX = resY = 100
----> 3 grid_x = np.linspace(xmin, xmax, resX)
4 grid_y = np.linspace(ymin, ymax, resY)
6 grid_z1 = fmesh.interpolate_to_grid(grid_x, grid_y, interp_c2f3)
NameError: name 'np' is not defined
fig, (ax1,ax2) = plt.subplots(1,2, figsize=(10,4))
im1 = ax1.imshow(grid_z1)
im2 = ax2.imshow(grid_z2)
ax1.set_title("grid_z1")
ax2.set_title("grid_z2")
fig.colorbar(im1, ax=ax1)
fig.colorbar(im2, ax=ax2)
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 fig, (ax1,ax2) = plt.subplots(1,2, figsize=(10,4))
3 im1 = ax1.imshow(grid_z1)
4 im2 = ax2.imshow(grid_z2)
NameError: name 'plt' is not defined
The next example is Ex4-Gradients