code: def morph(clm1, clm2, t, lmax):
      
      clm = (1 - t) * clm1 + t * clm2
      grid_reco = clm.expand(lmax=lmax)  # cut "high frequency" components
      agrid_reco = grid_reco.to_array()
      pts = []
      for i, longs in enumerate(agrid_reco):
          ilat = grid_reco.lats()[i]
          for j, value in enumerate(longs):
              ilong = grid_reco.lons()[j]
              th = (90 - ilat) / 57.3
              ph = ilong / 57.3
              r = value + rbias
              p = np.array([sin(th) * cos(ph), sin(th) * sin(ph), cos(th)]) * r
              pts.append(p)
      return pts
docs: Interpolate linearly the two sets of sph harm. coeeficients.