请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 164|回复: 3

[求助] python跟grads画出的图不一致

[复制链接]

新浪微博达人勋

发表于 2021-11-22 18:48:41 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
我在算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经度设置不加EW
    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])

###Z200TPSST的一元回归
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()


密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2021-11-22 18:50:46 | 显示全部楼层
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'
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-11-22 18:52:48 | 显示全部楼层
热带太平洋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
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-11-22 19:00:58 | 显示全部楼层
拜托了!!!

grads

grads

python

python

公式

公式
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表