How to use the pyshtools.SHGrid function in pyshtools

for th in np.linspace(0, np.pi, N, endpoint=False):
    lats = []
    for ph in np.linspace(0, 2 * np.pi, N, endpoint=False):
        p = spher2cart(rmax, th, ph)
        intersections = shape.intersectWithLine([0, 0, 0], p)  #
        if len(intersections):
            value = mag(intersections[0])
            lats.append(value)
            pts.append(intersections[0])
        else:
            lats.append(rmax)
            pts.append(p)
    agrid.append(lats)
agrid = np.array(agrid)

grid = pyshtools.SHGrid.from_array(agrid)
clm = grid.expand()
grid_reco = clm.expand(lmax=lmax)  # cut "high frequency" components


#############################################################
# interpolate to a finer grid
agrid_reco = grid_reco.to_array()

# adding 1 column
dim = (lmax * 2 + 2) + 1
gridTmp = np.zeros((dim - 1, dim))
gridTmp[:, 0:dim - 1] = agrid_reco
col0 = gridTmp[:, 0]
gridTmp[:, -1] = col0
agrid_reco = gridTmp

def sphrot_shtools(f, x, lmax=None, latcols=True):
    """ Rotate function on sphere f by Euler angles x (Z-Y-Z?)  """
    if 'pyshtools' not in sys.modules:
        raise ImportError('pyshtools not available!')

    if latcols:
        f = f.T
    c = SHGrid.from_array(f).expand()
    c_r = c.rotate(*x, degrees=False)
    f_r = c_r.expand(lmax=lmax).to_array()
    if latcols:
        f_r = f_r.T
    
    return f_r

shape1 = Sphere(alpha=0.2)
shape2 = vp.load(datadir + "icosahedron.vtk").normalize().lineWidth(1)

agrid1, actorpts1 = makeGrid(shape1, N)

vp.show(shape1, actorpts1, at=0)

agrid2, actorpts2 = makeGrid(shape2, N)
vp.show(shape2, actorpts2, at=1)

vp.camera.Zoom(1.2)
vp.interactive = False

clm1 = pyshtools.SHGrid.from_array(agrid1).expand()
clm2 = pyshtools.SHGrid.from_array(agrid2).expand()
# clm1.plot_spectrum2d()  # plot the value of the sph harm. coefficients
# clm2.plot_spectrum2d()

for t in arange(0, 1, 0.005):
    act21 = Points(morph(clm2, clm1, t, lmax), c="r", r=4)
    act12 = Points(morph(clm1, clm2, t, lmax), c="g", r=4)

    vp.show(act21, at=2, resetcam=0)
    vp.show(act12, at=3)
    vp.camera.Azimuth(2)

vp.show(interactive=1)

vp = Plotter(shape=[2, 2], verbose=0, axes=3, interactive=0)

shape1 = Sphere(alpha=0.2)
shape2 = vp.load(datadir + "icosahedron.vtk").normalize().lineWidth(1)

agrid1, actorpts1 = makeGrid(shape1, N)

vp.show(shape1, actorpts1, at=0)

agrid2, actorpts2 = makeGrid(shape2, N)
vp.show(shape2, actorpts2, at=1)

vp.camera.Zoom(1.2)
vp.interactive = False

clm1 = pyshtools.SHGrid.from_array(agrid1).expand()
clm2 = pyshtools.SHGrid.from_array(agrid2).expand()
# clm1.plot_spectrum2d()  # plot the value of the sph harm. coefficients
# clm2.plot_spectrum2d()

for t in arange(0, 1, 0.005):
    act21 = Points(morph(clm2, clm1, t, lmax), c="r", r=4)
    act12 = Points(morph(clm1, clm2, t, lmax), c="g", r=4)

    vp.show(act21, at=2, resetcam=0)
    vp.show(act12, at=3)
    vp.camera.Azimuth(2)

vp.show(interactive=1)