奇异值分解，脚本参考：

qxsw2016

奇异值分解的物理意义不大理解，只是按照魏凤英老师‘“现代气候统计诊断与预测技术”一书的实例进行了计算。用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

MeteoInfo

很好的工作！这种问题别人指望不上，自己好好钻研一下。