Example 4 - stripy gradients on the sphere

Example 4 - `stripy` gradients on the sphere#

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.

The next example is Ex5-Smoothing

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)

Analytic function#

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_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_ddlon = analytic_ddlon(mesh.lons, mesh.lats, 5.0, 2.0)
analytic_sol_ddlat = analytic_ddlat(mesh.lons, mesh.lats, 5.0, 2.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="lightgrey")
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, 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'

Derivatives of solution compared to analytic values#

The gradient_lonlat method of the sTriangulation takes a data array reprenting values on the mesh vertices and returns the lon and lat derivatives. There is an equivalent gradient_xyz method which returns the raw derivatives in Cartesian form. Although this is generally less useful, if you are computing the slope (for example) that can be computed in either coordinate system and may cross the pole, consider using the Cartesian form.

stripy_ddlon, stripy_ddlat = mesh.gradient_lonlat(analytic_sol)

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=analytic_sol,
                   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 = { "analytic": analytic_sol,
             "stripy ddlon": stripy_ddlon, 
             "stripy ddlat": stripy_ddlat, 
             "error ddlon":  stripy_ddlon-analytic_sol_ddlon, 
             "error ddlat":  stripy_ddlat-analytic_sol_ddlat }

@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
      2 plot = k3d.plot(camera_auto_fit=False, grid_visible=False, 
      3                 menu_visibility=True, axes_helper=False )
      5 indices = mesh.simplices.astype(np.uint32)

ModuleNotFoundError: No module named 'k3d'