Icosphere Transfer

Demonstrates data transfer between an Icosphere and a UniformMesh using nearest-neighbor, weighted, and RBF interpolation.

from sphedron import Icosphere, UniformMesh
from sphedron.transfer import MeshTransfer
from sphedron.extra import get_mesh_landmask
from sphedron.extra import plot_2d_mesh
import matplotlib.pyplot as plt
import numpy as np

factor = 32
icosphere_mesh = Icosphere(refine_factor=factor)
uniform_mesh = UniformMesh(resolution=1)

print(icosphere_mesh)
print(uniform_mesh)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[2], line 2
      1 factor = 32
----> 2 icosphere_mesh = Icosphere(refine_factor=factor)
      3 uniform_mesh = UniformMesh(resolution=1)
      5 print(icosphere_mesh)

TypeError: Mesh.__init__() got an unexpected keyword argument 'refine_factor'

plot_2d_mesh(
    icosphere_mesh,
    title=f"Full Icosphere, depth={icosphere_mesh.metadata['factor']}",
    scatter=True,
    s=0.5,
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 2
      1 plot_2d_mesh(
----> 2     icosphere_mesh,
      3     title=f"Full Icosphere, depth={icosphere_mesh.metadata['factor']}",
      4     scatter=True,
      5     s=0.5,
      6 )

NameError: name 'icosphere_mesh' is not defined

ico_land_mask = get_mesh_landmask(icosphere_mesh)
icosphere_mesh.mask_nodes(ico_land_mask)

plot_2d_mesh(
    icosphere_mesh,
    title=f"Land-masked Icosphere, depth={icosphere_mesh.metadata['factor']}",
    scatter=True,
    s=0.5,
)

plot_2d_mesh(
    uniform_mesh,
    title=f"Uniform mesh, resolution={uniform_mesh.resolution}",
    scatter=True,
    s=0.5,
)

The mesh nodes are in xyz, let’s create a test map where the values are |x|.

ico_values = np.abs(icosphere_mesh.nodes[:,0])

Now we want to transfer this map to the uniform mesh. We’ll use a nearest neighbor transfer, where the value at each uniform node is the average of its 3 nearest nodes from the icosphere. MeshTransfer will only compute the neighbors once.

transfer = MeshTransfer(icosphere_mesh, uniform_mesh, n_neighbors=5)

def weighted_func(w):
    min_w = np.min(w,axis=-1,keepdims=True) + 1e-16
    weights =  np.exp(-0.2*(w/min_w)**2)
    return weights / np.sum(weights, axis=-1, keepdims=True)

%%time
uniform_values = transfer.transfer(ico_values)

CPU times: user 65.2 ms, sys: 3.25 ms, total: 68.4 ms
Wall time: 10.7 ms

%%time
weighed_uniform_values = transfer.weighted_transfer(ico_values,weight_func=weighted_func)

CPU times: user 4.88 ms, sys: 2.09 ms, total: 6.97 ms
Wall time: 6.27 ms

The uniform values are flattened, we can reshape them as follows to get a 2D arrray:

uniform_values = uniform_mesh.reshape(uniform_values)
weighed_uniform_values = uniform_mesh.reshape(weighed_uniform_values)

plt.matshow(uniform_values[::-1])
plt.show()
plt.matshow(weighed_uniform_values[::-1])

<matplotlib.image.AxesImage at 0x7fdf66f0e8f0>

Let’s now mask land on the uniform mesh as well: first we compute the mask, then reshape it

uniform_mask = uniform_mesh.reshape(get_mesh_landmask(uniform_mesh))

uniform_values[uniform_mask] = -0.2
weighed_uniform_values[uniform_mask] = -0.2

plt.matshow(uniform_values[::-1])
plt.show()
plt.matshow(weighed_uniform_values[::-1])

<matplotlib.image.AxesImage at 0x7fdf552473d0>

true_values = uniform_mesh.reshape(np.abs(uniform_mesh.nodes[:,0]))

np.linalg.norm((weighed_uniform_values-true_values)[~uniform_mask])

np.float64(0.7653098799034589)

np.linalg.norm((uniform_values-true_values)[~uniform_mask])

np.float64(2.103699462953429)

from scipy.interpolate import RBFInterpolator

ico_coords = icosphere_mesh.nodes
ico_values = np.abs(icosphere_mesh.nodes[:, 0])
uni_coords = uniform_mesh.nodes


interpolated_values = RBFInterpolator(
    ico_coords,
    ico_values,
    neighbors=5,
    smoothing=0.1,
)(uni_coords)

interpolated_values = uniform_mesh.reshape(interpolated_values)

np.linalg.norm((interpolated_values-true_values)[~uniform_mask])

np.float64(0.32493005606554143)

interpolated_values[uniform_mask] = -0.2

plt.matshow(interpolated_values[::-1])

<matplotlib.image.AxesImage at 0x7fdf65fcd360>