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

axs, gs = plot_basis_functions(Y_nm, sampling)