Example 5 - stripy smoothing operations#
SSRFPACK is a Fortran 77 software package that constructs a smooth interpolatory or approximating surface to data values associated with arbitrarily distributed points on the surface of a sphere. It employs automatically selected tension factors to preserve shape properties of the data and avoid overshoot and undershoot associated with steep gradients.
Here we demonstrate how to access SSRFPACK smoothing through the stripy interface
Notebook contents#
The next example is Ex6-Scattered-Data
Define a computational mesh#
Use the (usual) icosahedron with face points included.
import stripy as stripy
mesh = stripy.spherical_meshes.icosahedral_mesh(refinement_levels=4, include_face_points=True)
print(mesh.npoints)
7682
Analytic function with noise and short wavelengths#
Define a relatively smooth function that we can interpolate from the coarse mesh to the fine mesh and analyse
import numpy as np
def analytic(lons, lats, k1, k2):
return np.cos(k1*lons) * np.sin(k2*lats)
def analytic_noisy(lons, lats, k1, k2, noise, short):
return np.cos(k1*lons) * np.sin(k2*lats) + short * (np.cos(k1*5.0*lons) * np.sin(k2*5.0*lats)) + noise * np.random.random(lons.shape)
# def analytic_ddlon(lons, lats, k1, k2):
# return -k1 * np.sin(k1*lons) * np.sin(k2*lats) / np.cos(lats)
# def analytic_ddlat(lons, lats, k1, k2):
# return k2 * np.cos(k1*lons) * np.cos(k2*lats)
analytic_sol = analytic(mesh.lons, mesh.lats, 5.0, 2.0)
analytic_sol_n = analytic_noisy(mesh.lons, mesh.lats, 5.0, 2.0, 0.1, 0.0)
%matplotlib inline
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10, 10), facecolor="none")
ax = plt.subplot(111, projection=ccrs.Orthographic(central_longitude=0.0, central_latitude=0.0, globe=None))
ax.coastlines(color="#999999", linewidth=2.0)
ax.set_global()
lons0 = np.degrees(mesh.lons)
lats0 = np.degrees(mesh.lats)
ax.scatter(lons0, lats0,
marker="o", s=10.0, transform=ccrs.PlateCarree(), c=analytic_sol_n-analytic_sol, cmap=plt.cm.RdBu)
pass
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
3 import cartopy
4 import cartopy.crs as ccrs
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'
Smoothing operations#
The sTriangulation.smoothing method directly wraps the SSRFPack smoother that smooths a surface f described
by values on the mesh vertices to find a new surface f’ (also described on the mesh vertices) by choosing nodal function values and gradients to minimize the linearized curvature of F subject to a bound on the deviation from the data values.
help(mesh.smoothing)
smoothing(self, f, w, sm, smtol, gstol)
method of stripy.spherical_meshes.icosahedral_mesh instance
Smooths a surface f by choosing nodal function values and gradients to
minimize the linearized curvature of F subject to a bound on the
deviation from the data values. This is more appropriate than interpolation
when significant errors are present in the data.
Parameters
----------
f : array of floats, shape (n,)
field to apply smoothing on
w : array of floats, shape (n,)
weights associated with data value in f
w[i] = 1/sigma_f^2 is a good rule of thumb.
sm : float
positive parameter specifying an upper bound on Q2(f).
generally n-sqrt(2n) <= sm <= n+sqrt(2n)
smtol : float
specifies relative error in satisfying the constraint
sm(1-smtol) <= Q2 <= sm(1+smtol) between 0 and 1.
gstol : float
tolerance for convergence.
gstol = 0.05*mean(sigma_f)^2 is a good rule of thumb.
Returns
-------
f_smooth : array of floats, shape (n,)
smoothed version of f
(dfdx, dfdy, dfdz) : tuple of floats, tuple of 3 shape (n,) arrays
first derivatives of f_smooth in the x, y, and z directions
err : error indicator
stripy_smoothed, dds, err = mesh.smoothing(analytic_sol_n, np.ones_like(analytic_sol_n), 10.0, 0.1, 0.01)
stripy_smoothed2, dds, err = mesh.smoothing(analytic_sol_n, np.ones_like(analytic_sol_n), 1.0, 0.1, 0.01)
stripy_smoothed3, dds, err = mesh.smoothing(analytic_sol_n, np.ones_like(analytic_sol_n), 50.0, 0.1, 0.01)
delta_n = analytic_sol_n - stripy_smoothed
delta_ns = analytic_sol - stripy_smoothed
delta_n2 = analytic_sol_n - stripy_smoothed2
delta_ns2 = analytic_sol - stripy_smoothed2
delta_n3 = analytic_sol_n - stripy_smoothed3
delta_ns3 = analytic_sol - stripy_smoothed3
Results of smoothing with different value of sm#
import k3d
plot = k3d.plot(camera_auto_fit=False, grid_visible=False,
menu_visibility=True, axes_helper=False )
indices = mesh.simplices.astype(np.uint32)
points = np.column_stack(mesh.points.T).astype(np.float32)
mesh_viewer = k3d.mesh(points, indices, wireframe=False, attribute=stripy_smoothed,
color_map=k3d.colormaps.basic_color_maps.CoolWarm,
name="original",
flat_shading=False, opacity=1.0 )
plot += mesh_viewer
plot += k3d.points(points, point_size=0.01,color=0x779977)
plot.display()
## ## ##
from ipywidgets import interact, interactive
import ipywidgets as widgets
choices = { "unsmoothed": analytic_sol_n,
"smooth1": stripy_smoothed,
"smooth2": stripy_smoothed2,
"smooth3": stripy_smoothed3,
"Delta smooth1": analytic_sol_n - stripy_smoothed,
"Delta smooth2": analytic_sol_n - stripy_smoothed2,
"Delta smooth3": analytic_sol_n - stripy_smoothed3, }
@interact(choice=choices.keys())
def chooser(choice):
mesh_viewer.attribute = choices[choice].astype(np.float32)
range = np.sqrt((choices[choice]**2).mean()) * 0.5
mesh_viewer.color_range = [-range, range]
return
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[5], line 1
----> 1 import k3d
3 plot = k3d.plot(camera_auto_fit=False, grid_visible=False,
4 menu_visibility=True, axes_helper=False )
6 indices = mesh.simplices.astype(np.uint32)
ModuleNotFoundError: No module named 'k3d'
The next example is Ex6-Scattered-Data