- 积分
- 1626
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-1-5
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
奇异值分解的物理意义不大理解,只是按照魏凤英老师‘“现代气候统计诊断与预测技术”一书的实例进行了计算。用1951-1990年中国夏季降水和夏季北美陆地气温进行Svd分析。由于没有英国哈德来中心的气候资料,所以以美国再分析资料来代替。导致分析结果有所不同,不知是数据问题,还是编程计算问题。也请各位帮忙指点。脚本如下:
from mipylib.numeric import stats
#左变量场X
f = addfile(r'F:\DataSet\air.mon.mean.nc')
fn = 'F:/DataSet/r160%d.txt'
stfn = 'F:/DataSet/st_160.csv'
Lfields = f['air'][0,'2.5:82.5','227.5:337.5']
tempJJA = [] #1951-1990年北美夏季温度
for i in range(1951,1991):
tidx = (i-1948)*12 #每年一月
t = f.gettime(tidx)
data = f['air'][tidx+5,'2.5:82.5','227.5:337.5']+\
f['air'][tidx+6,'2.5:82.5','227.5:337.5']+\
f['air'][tidx+7,'2.5:82.5','227.5:337.5']
tempJJA.append(data[::2,::2]) #间隔5度
nx = Lfields[::2,::2] .dimlen(1)
ny = Lfields[::2,::2] .dimlen(0)
var = Lfields[::2,::2]
n = len(tempJJA)
p = nx*ny
#print n,p
xmartix = zeros([p,n])
k = 0
for i in arange(0,nx):
for j in arange(0,ny):
x = []
k = j+i*ny
for t in arange(1951,1991):
tidx = t -1951
x.append(tempJJA[tidx][j,i])
xm = sum(x)/n
s2 = stats.covariance(x,x,bias = True)
#print xm,s2
#资料标准化
for t in range(0,n):
x[t] = (x[t] - xm)/sqrt(s2)
xmartix[k,:] = x[:]
#print xmartix[0,:]
#右变量场Y
datafn = 'F:/DataSet/r1606.txt'
stn = 160
#Read pre data of 160 stations from data file
pre = asciiread(datafn, shape=(n,stn))
mons = []
for i in range(6, 9):
mons.append(i)
monn = len(mons)
stn = 160
#Read preciptation data of 160 stations from data file
prep = zeros((3,n,stn))
prepJJA = zeros((stn,n))
for i in range(6, 9):
datafn = fn%i
#print datafn
data = asciiread(datafn, shape=(n,stn))
#print data
prep[i-6,:,:] =data
for i in range(0, n):
for j in range(0,stn):
prepJJA[j,i] = sum(prep[:,i,j])
ymartix = zeros([stn,n])
for i in arange(0,stn):
y = []
for j in arange(0,n):
y.append( prepJJA[i,j])
ym = sum(y)/n
s2 = stats.covariance(y,y,bias = True)
#资料标准化
for t in range(0,n):
y[t] = (y[t] - ym)/sqrt(s2)
ymartix[i,:] = y[:]
#交叉协方差矩阵
S = np.dot(xmartix,ymartix.T)/n #p*stn
#奇异值分解
L,lamda,R = np.linalg.svd(S)
#print lamda[0]*100/sum(lamda)
gridL = zeros([ny,nx])
for i in arange(0,nx):
for j in arange(0,ny):
gridL[j,i] = L[j+i*ny,0]
gridL = dim_array(gridL,var.dims)
#Plot
subplot(2,1,1)#中国夏季降水
#Read station name and lon/lat
table = readtable(stfn, delimiter=',', format='%i%2s%2f')
stnames = table['Name']
lat = table['LAT']
lon = table['LON']
#To grid data
x = arange(70, 140, 0.5)
y = arange(15, 58, 0.5)
gridR,gx,gy = griddata((lon, lat), R.T[:,0], xi=(x, y), method='idw', pointnum=4)
#Plot
axesm(newaxes=False)
mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
geoshow(mlayer,edgecolor=(115,115,115))
layer = contourfm(x, y, gridR,20)
colorbar(layer,orientation='vertical')
title(u'中国夏季降水(上)与北美夏季陆地温度(下)第一空间分布型',fontname=u'宋体', fontsize=16)
subplot(2,1,2)#北美夏季陆地温度
axesm(newaxes=False)
#axesm(tickfontsize=20)
#mlayer = shaperead('F:\MeteoInfo\MeteoInfo\Map\country1.shp')
geoshow(mlayer,edgecolor=(115,115,115))
layer = contourfm(gridL,20)
colorbar(layer,orientation='vertical')
lamda[lamda<0] = 0
print u'第一空间分布型解释方差:(%)',lamda[0]*100/sum(lamda)
U = np.dot(L[:,0],xmartix)
V = np.dot(R.T[:,0],ymartix)
slope, intercept, r, p, std,number = linregress(U, V)
print u'第一空间分布型相关系数:',r
|
-
|