python-为磁场强度数据返回对称轮廓图的底图

为没有发布图片表示歉意，这篇报道纯粹是为了这个问题。我的目标是得到一个场强图，看起来像其中一个：其他图片_skf/41586 _2018 _468 _Fig1 _HTML。巴布亚新几内亚

我得到的情节在南半球似乎是正确的，但随后在北半球得到了反映。不确定这是我把等式编错了还是用了Basemap错了，但用下面的代码我得到了下面的图。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pyshtools.legendre.legendre_functions import legendre_lm
#fixing system specific problem
import os
os.environ["PROJ_LIB"] = "C:\\Utilities\\Python\\Anaconda\\Library\\share"; #fixr
#end of fix
from mpl_toolkits.basemap import Basemap
data = pd.read_csv('Juno data.csv')
#create two data sets for g and h coefficients
gdata = data[(data['g h'] == 'g')]
hdata = data[(data['g h'] == 'h')]
nmax = data['n'].max()
mmax =data['m'].max()
r = 7.1492E6
Rs = 7E6
#latitude
th = np.linspace(-np.pi/2,np.pi/2,360)
#longitude
ph = np.linspace(-np.pi,np.pi,360)
lon,lat = np.meshgrid(ph*180/np.pi,th*180/np.pi)
#calculating potential
prompt = input("Enter the maximum harmonic degree/n: ")
deg = int(prompt)
V = 0
for n in range(1,deg+1):
    for m in range(0,n):
        #creating g and h gauss coefficients as variables
        g = gdata.loc[(gdata.n==n)&(gdata.m==m),'Value'].values[0]
        try: 
            hdata.loc[(hdata.n==n)&(hdata.m==m),'Value'].values[0]
        except:
            h = 0
        else:
            h = hdata.loc[(hdata.n==n)&(hdata.m==m),'Value'].values[0]
        V = V + ((Rs/r)**(n+1))*legendre_lm(n,m,np.cos(np.deg2rad(lat)),['schmidt'])*(g*np.cos(m*np.deg2rad(lon))+h*np.sin(m*np.deg2rad(lon)))
B = ((nmax+1)/r)*V*Rs
fig = plt.figure(figsize=(12,8))
m = Basemap(projection = 'moll',lon_0 = 0)
m.drawmeridians(np.arange(-180,180,60),labels=[1,0,0,1])
m.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])

x,y = m(lon,lat)
m.contourf(x,y,B)
m.colorbar(location='right')

#data south of the equator seems to be reflected in the equator???

我试图绘制的方程是纬度和经度的函数