Y_nm = spharpy.SphericalHarmonics(
    n_max=n_max, coordinates=sampling,
    basis_type='real',
    normalization='N3D',
    condon_shortley='auto'
).basis

axs, gs = plot_basis_functions(Y_nm, sampling)