python跟grads画出的图不一致

linzy

我在算Elnino海温指数（热带太平洋的区域平均）一个时间序列与这些年份对应的全球200hPa位势高度场的一元回归的时候发现，两者画出的图完全不一致，我找了两天了，完全没发现问题出在哪。grads和python用到的公式都是课本上的，我觉得应该没道理会画出两张不一样的图呀。拜托各位前辈帮我看一看。
python：
from scipy import stats
from xgrads import *
import numpy as np  # 调用numpy
import xarray as xr
from math import *
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter


def createmap():
    ###############################################生成地图##########################################################
    box = [-180, 180, -90, 90]  # 经度维度
    scale = '110m'  # 地图分辨率
    xstep = 20  # 下面标注经纬度的步长
    ystep = 10
    proj = ccrs.PlateCarree(central_longitude=180)  # 确定地图投影
    fig = plt.figure(figsize=(9, 6))  # dpi=150)###生成底图
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  # 确定子图，与grads的类似
    ##海岸线
    ax.coastlines(scale)
    # 标注坐标轴
    ax.set_xticks(np.arange(box[0], box[1] + xstep, xstep), crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(box[2], box[3] + ystep, ystep), crs=ccrs.PlateCarree())
    # 经纬度格式，把0经度设置不加E和W
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ############################################################################################################
    return ax, fig


############NCEP_Z200_30y_Wt转换成nc文件############
hgt = open_CtlDataset('D:\\grads\\TongJi\\NCEP_Z200_30y_Wt.ctl')
hgt.attrs['pdef'] = 'None'
hgt.to_netcdf('D:\\grads\\TongJi\\NCEP_Z200_30y_Wt.nc')
hgt = xr.open_dataset('D:\\grads\\TongJi\\NCEP_Z200_30y_Wt.nc')
############NCEP_TPSST_30y_Wt转换成nc文件###########
air = open_CtlDataset('D:\\grads\\TongJi\\NCEP_TPSST_30y_Wt.ctl')
air.attrs['pdef'] = 'None'
air.to_netcdf('D:\\grads\\TongJi\\NCEP_TPSST_30y_Wt.nc')
air = xr.open_dataset('D:\\grads\\TongJi\\NCEP_TPSST_30y_Wt.nc')
###########NCEP_IOSST_30y_Wt转换成nc文件###########
airindia = open_CtlDataset('D:\\grads\\TongJi\\NCEP_IOSST_30y_Wt.ctl')
airindia.attrs['pdef'] = 'None'
airindia.to_netcdf('D:\\grads\\TongJi\\NCEP_IOSST_30y_Wt.nc')
airindia = xr.open_dataset('D:\\grads\\TongJi\\NCEP_IOSST_30y_Wt.nc')

#####热带太平洋海温
lon = air['lon'][:]
lat = air['lat'][:]
st1 = air['st'][4, :, :]
st2 = air['st'][8, :, :]
st3 = air['st'][13, :, :]
st4 = air['st'][19, :, :]
st5 = air['st'][24, :, :]

for i in range(len(lon)):
    for j in range(len(lat)):
        if float(st1[j, i]) > 9999:
            st1[j, i] = np.NAN
        if float(st2[j, i]) > 9999:
            st2[j, i] = np.NAN
        if float(st3[j, i]) > 9999:
            st3[j, i] = np.NAN
        if float(st4[j, i]) > 9999:
            st4[j, i] = np.NAN
        if float(st5[j, i]) > 9999:
            st5[j, i] = np.NAN

avest1 = np.nanmean(st1)
avest2 = np.nanmean(st2)
avest3 = np.nanmean(st3)
avest4 = np.nanmean(st4)
avest5 = np.nanmean(st5)
###印度太平洋海温###
lon2 = airindia['lon'][:]
lat2 = airindia['lat'][:]
st_india1 = airindia['st'][4, :, :]
st_india2 = airindia['st'][8, :, :]
st_india3 = airindia['st'][13, :, :]
st_india4 = airindia['st'][19, :, :]
st_india5 = airindia['st'][24, :, :]

for i in range(len(lon2)):
    for j in range(len(lat2)):
        if float(st_india1[j, i]) > 9999:
            st_india1[j, i] = np.NAN
        if float(st_india2[j, i]) > 9999:
            st_india2[j, i] = np.NAN
        if float(st_india3[j, i]) > 9999:
            st_india3[j, i] = np.NAN
        if float(st_india4[j, i]) > 9999:
            st_india4[j, i] = np.NAN
        if float(st_india5[j, i]) > 9999:
            st_india5[j, i] = np.NAN
avest_india1 = np.nanmean(st_india1)
avest_india2 = np.nanmean(st_india2)
avest_india3 = np.nanmean(st_india3)
avest_india4 = np.nanmean(st_india4)
avest_india5 = np.nanmean(st_india5)
###200hPa位势高度场
time = hgt['time'][:]
lon1 = hgt['lon'][:]
lat1 = hgt['lat'][:]
hgt1 = hgt['hgt'][4, :, :]
hgt2 = hgt['hgt'][8, :, :]
hgt3 = hgt['hgt'][13, :, :]
hgt4 = hgt['hgt'][19, :, :]
hgt5 = hgt['hgt'][24, :, :]
avehgt = (hgt1 + hgt2 + hgt3 + hgt4 + hgt5) / 5

avest = (avest1 + avest2 + avest3 + avest4 + avest5) / 5
avest_india = (avest_india1 + avest_india2 + avest_india3 + avest_india4 + avest_india5) / 5

stTP = np.array([[avest1], [avest2], [avest3], [avest4], [avest5]])
stindia = np.array([[avest_india1], [avest_india2], [avest_india3], [avest_india4], [avest_india5]])
hgtZ200 = np.array([hgt1, hgt2, hgt3, hgt4, hgt5])

###Z200和TPSST的一元回归
slope = np.zeros((len(lat1), len(lon1)))
intercept = np.zeros((len(lat1), len(lon1)))
r_value = np.zeros((len(lat1), len(lon1)))
p_value = np.zeros((len(lat1), len(lon1)))
std_err = np.zeros((len(lat1), len(lon1)))
b = np.zeros((len(lat1), len(lon1)))

for i in range(len(lon1)):
    for j in range(len(lat1)):
        b[j, i] = (avest1 * hgt1[j, i] + avest2 * hgt2[j, i] + avest3 * hgt3[j, i] + avest4 * hgt4[j, i] + avest5 *
                   hgt5[j, i] - 5 * avest * avehgt[j, i]) / (
                          pow(avest1, 2) + pow(avest2, 2) + pow(avest3, 2) + pow(avest4, 2) + pow(avest5, 2) - 5 * pow(
                      avest, 2))
ax, fig = createmap()
ax.contourf(lon1, lat1, b, transform=ccrs.PlateCarree())
plt.show()

linzy

grads gs脚本：

'reinit'
'set display color white'
'c'
'open D:\grads\TongJi\NCEP_TPSST_30y_Wt.ctl'
'open D:\grads\TongJi\NCEP_IOSST_30y_Wt.ctl'
'open D:\grads\TongJi\NCEP_Z200_30y_Wt.ctl'
'set lon 0 360'
'set lat -90 90'
'set t 5'
'define TP1=aave(st.1,lon=120,lon=300,lat=-20,lat=20)'
'define IO1=aave(st.2,lon=35,lon=125,lat=-20,lat=20)'
'define Z1=hgt.3'
'set t 9'
'define TP2=aave(st.1,lon=120,lon=300,lat=-20,lat=20)'
'define IO2=aave(st.2,lon=35,lon=125,lat=-20,lat=20)'
'define Z2=hgt.3'
'set t 14'
'define TP3=aave(st.1,lon=120,lon=300,lat=-20,lat=20)'
'define IO3=aave(st.2,lon=35,lon=125,lat=-20,lat=20)'
'define Z3=hgt.3'
'set t 20'
'define TP4=aave(st.1,lon=120,lon=300,lat=-20,lat=20)'
'define IO4=aave(st.2,lon=35,lon=125,lat=-20,lat=20)'
'define Z4=hgt.3'
'set t 25'
'define TP5=aave(st.1,lon=120,lon=300,lat=-20,lat=20)'
'define IO5=aave(st.2,lon=35,lon=125,lat=-20,lat=20)'
'define Z5=hgt.3'
'define aveTP=(TP1+TP2+TP3+TP4+TP5)/5'
'define aveIO=(IO1+IO2+IO3+IO4+IO5)/5'
'define avehgt=(Z1+Z2+Z3+Z4+Z5)/5'
***Z200andTPSST***
'define b=(TP1*Z1+TP2*Z2+TP3*Z3+TP4*Z4+TP5*Z5-5*aveTP*avehgt)/(pow(TP1,2)+pow(TP2,2)+pow(TP3,2)+pow(TP4,2)+pow(TP5,2)-5*pow(aveTP,2))'
'set gxout shaded'
'd b'
'cbar 1 0'
'draw title Z200andTPSST'
'gxprint D:\grads\TongJi\ex4\Z200andTPSST.png'

linzy

热带太平洋30年海温场描述文件
dset D:\grads\TongJi\NCEP_TPSST_30y_Wt.dat
undef  1.0E+33
title  NCEP Data
xdef  32  linear  121.88  5.650
ydef   7  levels     -18.0950  -12.3808   -6.6666   -0.9524    4.7618   10.4761   16.1902
zdef   1  levels   1000
tdef  30  linear  0z01Jan1979  1yr
vars   1
st   0   99  TPSST
endvars
全球200hPa位势高度描述文件：
dset D:\grads\TongJi\NCEP_Z200_30y_Wt.dat
undef  1.0E+33
title  NCEP data
xdef  48  levels   5          12.5            20          27.5            35          42.5            50          57.5            65          72.5            80          87.5            95         102.5           110         117.5           125         132.5           140         147.5           155         162.5           170         177.5           185         192.5           200         207.5           215         222.5           230         237.5           245         252.5           260         267.5           275         282.5           290         297.5           305         312.5           320         327.5           335         342.5           350        356.25
ydef  24  levels   -87.5          -80        -72.5          -65        -57.5          -50        -42.5          -35        -27.5          -20        -12.5           -5          2.5           10         17.5           25         32.5           40         47.5           55         62.5           70         77.5           85
zdef   1  levels   1000
tdef  30  linear  0z01Jan1978  1yr
vars   1
hgt   0   99  Z200 Winter
endvars

linzy

拜托了！！！