脚本示例：线性相关分析Nino3.4海温指数与同期降水

qxsw2016 · 发表于 2018-4-27 20:18:07

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

您需要登录才可以下载或查看，没有帐号？立即注册

x

线性相关分析中用同期Nino3.4海温指数与同期降水/气温场做相关，计算相关系数前去掉了它们的线性趋势：
table = readtable('f:/nino34.csv', delimiter=',', format='%s%f')
yy = table['yyJJA']
nino34 = table['nino34']
#print nino34
f = addfile('f:/precip.mon.mean.2.5x2.5.nc')
#prep = f['precip'][0::12,'0:60','60:150']
var  = f['precip'][0,'0:60','60:150']
prepJJA = []
for i in range(1981,2016):
tidx = (i-1948)*12
t = f.gettime(tidx)
data = f['precip'][tidx+5,'0:60','60:150']+\
         f['precip'][tidx+6,'0:60','60:150']+\
         f['precip'][tidx+7,'0:60','60:150']
prepJJA.append(data)
nt = len(prepJJA)
#trend
for i in arange(0,nx):
for j in arange(0,ny):
      xx = []
      yy = []
      for t in syears:
         tidx = t -1981
         if(prepJJA[tidx][j,i]>-9.96921E36):
            xx.append(t)
            yy.append(prepJJA[tidx][j,i])
      y = zeros([nt])
      if(len(xx)>1):
         slope, intercept, r, p, std,number = linregress(xx, yy)
         #print slope, intercept, r
         y = polyval([slope, intercept], syears)
      else:
         y = nan
      #print i,j,y
      #1981-2015
      for t in syears:
         tidx = t -1981
         if(prepJJA[tidx][j,i]>-9.96921E36):
            prepJJA[tidx][j,i] = prepJJA[tidx][j,i]-y[tidx]
         else:
            prepJJA[tidx][j,i] = nan

syears = arange(1981,2016)
nx = var.dimlen(1)
ny = var.dimlen(0)
cor = zeros([ny,nx])
cor = dim_array(cor,var.dims)
for i in arange(0,nx):
for j in arange(0,ny):
      x = []
      y = []
      for t in syears:
         tidx = t -1981
         if(prepJJA[tidx][j,i]>-9.96921E36):
            x.append(nino34[t-1981])
            y.append(prepJJA[tidx][j,i])
      #scatter(xx, yy, s=2, color='k', label='PM')
      if(len(x)>1):
         slope, intercept, r, p, std,number = linregress(x, y)
         cor[j,i] = r
      else:
         cor[j,i] = nan
#print cor
axesm()
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
clevs =[-0.53,-0.43,-0.33,-0.28,0,0.28,0.33,0.43,0.53]
cols=[(86,49,5),(138,80,0),(186,145,42),(216,179,107),(245,236,212),\
(217,238,235),(117,197,186),(46,145,137),(0,99,91),(0,59,47)]
#layer = contourfm(cor,20)
layer = contourfm(cor,clevs,colors=cols)
colorbar(layer,orientation='horizontal')
xlim(60,150)
ylim(0,60)
title('Correlation:nino34(JJA)-prep(JJA) 1981-2015')

qxsw2016 · 发表于 2018-4-29 14:40:24

增加了简单的注释，请批评指正，参考脚本；table = readtable('f:/nino34.csv', delimiter=',', format='%s%f')
yy = table['yyJJA']
nino34 = table['nino34']
#print nino34
f = addfile('f:/precip.mon.mean.2.5x2.5.nc')
#prep = f['precip'][0::12,'0:60','60:150']
var  = f['precip'][0,'0:60','60:150']
syears = arange(1981,2016)
nx = var.dimlen(1)
ny = var.dimlen(0)
prepJJA = []  #1981-2015年夏季降水序列
syears = arange(1981,2016)
for i in range(1981,2016):
tidx = (i-1948)*12 #每年1月
t = f.gettime(tidx)
data = f['precip'][tidx+5,'0:60','60:150']+\
         f['precip'][tidx+6,'0:60','60:150']+\
         f['precip'][tidx+7,'0:60','60:150']
prepJJA.append(data)
nt = len(prepJJA)
#trend固定某一格点，得到年份和降水两个时间序列，计算两者线性回归拟合值
for i in arange(0,nx):
for j in arange(0,ny):
      xx = []
      yy = []
      for t in syears:
         tidx = t -1981
         if(prepJJA[tidx][j,i]>-9.96921E36):
            xx.append(t)
            yy.append(prepJJA[tidx][j,i])
      y = zeros([nt])
      if(len(xx)>1):
         slope, intercept, r, p, std,number = linregress(xx, yy)
         #print slope, intercept, r
         y = polyval([slope, intercept], syears)
      else:
         y = nan
      #print i,j,y
      #1981-2015去掉了降水的线性拟合
      for t in syears:
         tidx = t -1981
         if(prepJJA[tidx][j,i]>-9.96921E36):
            prepJJA[tidx][j,i] = prepJJA[tidx][j,i]-y[tidx]
         else:
            prepJJA[tidx][j,i] = nan
#构建相关系数二维场（年份序列和去掉了它们的线性趋势的降水序列）
cor = zeros([ny,nx])
cor = dim_array(cor,var.dims)
for i in arange(0,nx):
for j in arange(0,ny):
      x = []
      y = []
      for t in syears:
         tidx = t -1981
         if(prepJJA[tidx][j,i]>-9.96921E36):
            x.append(nino34[t-1981])
            y.append(prepJJA[tidx][j,i])
      #scatter(xx, yy, s=2, color='k', label='PM')
      if(len(x)>1):
         slope, intercept, r, p, std,number = linregress(x, y)
         cor[j,i] = r
      else:
         cor[j,i] = nan
#print cor
axesm()
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
geoshow(mlayer,edgecolor=(115,115,115))
clevs =[-0.53,-0.43,-0.33,-0.28,0,0.28,0.33,0.43,0.53]
cols=[(86,49,5),(138,80,0),(186,145,42),(216,179,107),(245,236,212),\
(217,238,235),(117,197,186),(46,145,137),(0,99,91),(0,59,47)]
#layer = contourfm(cor,20)
layer = contourfm(cor,clevs,colors=cols)
colorbar(layer,orientation='horizontal')
xlim(60,150)
ylim(0,60)
title('Correlation:nino34(JJA)-prep(JJA) 1981-2015')

MeteoInfo · 发表于 2018-4-29 11:10:53

很好的帖子！再增加些说明就更好了，还有图里怎么没有画出地图？

cgc_ustc · 发表于 2018-4-29 14:50:52

有道理

p_yf100 · 发表于 2018-5-3 17:45:10

学习学习一下，很好，课题研究正好需要借鉴，谢谢楼主分享。

wkszlhd · 发表于 2018-5-4 10:22:49

学习了谢谢

MeteoInfo · 发表于 2018-5-4 11:13:40

qxsw2016 发表于 2018-4-29 14:40
增加了简单的注释，请批评指正，参考脚本；table = readtable('f:/nino34.csv', delimiter=',', format='%s ...

已放入脚本汇总贴中。

qxsw2016 · 发表于 2018-10-3 09:03:53

本帖最后由 qxsw2016 于 2018-10-3 09:07 编辑

关于contour函数参数使用问题，这个函数定义如下：
def contour(*args, **kwargs):
"""
Plot contours.

:param x: (*array_like*) Optional. X coordinate array.
:param y: (*array_like*) Optional. Y coordinate array.
:param z: (*array_like*) 2-D z value array.
:param levs: (*array_like*) Optional. A list of floating point numbers indicating the level curves
      to draw, in increasing order.
:param cmap: (*string*) Color map string.
:param colors: (*list*) If None (default), the colormap specified by cmap will be used. If a
      string, like 鈥榬鈥?or 鈥榬ed鈥? all levels will be plotted in this color. If a tuple of matplotlib
      color args (string, float, rgb, etc), different levels will be plotted in different colors in
      the order specified.
:param smooth: (*boolean*) Smooth countour lines or not.

:returns: (*VectoryLayer*) Contour VectoryLayer created from array data.

我写的脚本，等值线设置不起作用，还用就是好像是连个图叠加是怎么回事？如果格点数据有nan值怎么判断。脚本如下：
f = addfile('F:/DataSet/ncep/hgt.mon.mean.nc')
var = f['hgt'][0,'500','0:60','100:180']
print var
axesm()
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
geoshow(mlayer)
clevs =[5600,5640,5680,5720,5760,5800,5840,5880]
cols=[(182,106,40),(205,133,63),(225,165,100),(245,205,132),\
(255,245,186),(205,255,205),(153,240,178),(83,189,159)]
layer = contour(var,clevs, colors=cols)
xlim(100,180)
ylim(0,60)
clabel(layer)

MeteoInfo · 发表于 2018-10-3 10:21:24

qxsw2016 发表于 2018-10-3 09:03
关于contour函数参数使用问题，这个函数定义如下：
def contour(*args, **kwargs):
"""

在地图上绘制等值线要用 contourm 函数，地图上的绘图函数通常会有 m 作为结尾以和一般坐标系区分。

qxsw2016 · 发表于 2018-10-3 17:27:10

多谢老师，假期愉快！还有一个问题，有没有别的解决方法。就是达到某个阈值填色，其他不填色，效果如图。我写了个脚本：f = addfile('F:/DataSet/ncep/hgt.mon.mean.nc')
var = f['hgt'][0,'500',:,:]
axesm()
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
geoshow(mlayer)
clevs =[5440,5480,5520,5560,5600,5640,5680,5720,5760,5800,5840,5880]
#cols=[(182,106,40),(205,133,63),(225,165,100),(245,205,132),\
# (255,245,186),(205,255,205),(153,240,178),(83,189,159)]
cols=[(0,0,0),(0,0,0),(0,0,0),(0,0,0),(0,0,0),(0,0,0),\
(0,0,0),(0,0,0),(0,0,0),(0,0,0),(0,0,0),(0,0,0)]
layer = contourm(var,clevs, colors=cols)
clabel(layer)
cols2=[(0,0,255),(255,255,255),(255,0,0)]
clevs2 =[5440,5840]
layer2 = contourfm(var,clevs2, colors=cols2)
#xlim(100,180)
#ylim(0,50)
colorbar(layer2)

		自动登录	找回密码
密码			立即注册