S
Sudheer Joseph
HI,
I have been trying to compute cross correlation between a time series at a location f(1) and the timeseries of spatial data f(XYT) and saving the resulting correlation coefficients and lags in a 3 dimensional array which is of fairly big size. Though the code I made for this purpose works upto few iterations then it hangs due to apparent memory crunch. Can anybodysuggest a better way to handle this situation so that the computation and data storing can be done with out hangups. Finally I intend to save the data as netcdf file which is not implemented as of now. Below is the piece of code I wrote for this purpose.
from mpl_toolkits.basemap import Basemap as bm, shiftgrid, cm
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from math import pow, sqrt
import sys
from scipy.stats import t
indep=120
nlags=365
ncin = Dataset('qu_ru.nc', 'r')
lons = ncin.variables['LON421_600'][:]
lats = ncin.variables['LAT81_220'][:]
dep = ncin.variables['DEPTH1_29'][:]
adep=(dep==indep).nonzero()
didx=int(adep[0])
qu = ncin.variables['qu'][:,:,:]
#qv = ncin.variables['QV'][0,:,:]
ru = ncin.variables['ru'][:,didx,0,0]
ncin.close()
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# use major and minor sphere radii from WGS84 ellipsoid.
m = bm(projection='cyl', llcrnrlon=30, llcrnrlat=-40,urcrnrlon=120, urcrnrlat=30)
# transform to nx x ny regularly spaced 5km native projection grid
nx = int((m.xmax-m.xmin))+1; ny = int((m.ymax-m.ymin)+1)
q=ru[1:2190]
qmean=np.mean(q)
qstd=np.std(q)
qnorm=(q-qmean)/qstd
lags3d=np.arange(731*140*180).reshape(731,140,180)
r3d=np.arange(731*140*180).reshape(731,140,180)
for i in np.arange(len(lons)):
for j in np.arange(len(lats)):
print i,j
p=qu[1:2190,j,i].squeeze()
p.shape
pmean=np.mean(p)
pstd=np.std(p)
pnorm=(p-pmean)/pstd
n=len(p)
# fg=plt.figure()
c=plt.xcorr(p,q,usevlines=True,maxlags=nlags,normed=True,lw=2)
acp=plt.acorr(p,usevlines=True,maxlags=nlags,normed=True,lw=2)
acq=plt.acorr(q,usevlines=True,maxlags=nlags,normed=True,lw=2)
acp[1][nlags]=0
acq[1][nlags]=0
lags=c[0]
r=c[1]
lags3d[:,j,i]=lags
r3d[:,j,i]=r
I have been trying to compute cross correlation between a time series at a location f(1) and the timeseries of spatial data f(XYT) and saving the resulting correlation coefficients and lags in a 3 dimensional array which is of fairly big size. Though the code I made for this purpose works upto few iterations then it hangs due to apparent memory crunch. Can anybodysuggest a better way to handle this situation so that the computation and data storing can be done with out hangups. Finally I intend to save the data as netcdf file which is not implemented as of now. Below is the piece of code I wrote for this purpose.
from mpl_toolkits.basemap import Basemap as bm, shiftgrid, cm
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from math import pow, sqrt
import sys
from scipy.stats import t
indep=120
nlags=365
ncin = Dataset('qu_ru.nc', 'r')
lons = ncin.variables['LON421_600'][:]
lats = ncin.variables['LAT81_220'][:]
dep = ncin.variables['DEPTH1_29'][:]
adep=(dep==indep).nonzero()
didx=int(adep[0])
qu = ncin.variables['qu'][:,:,:]
#qv = ncin.variables['QV'][0,:,:]
ru = ncin.variables['ru'][:,didx,0,0]
ncin.close()
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# use major and minor sphere radii from WGS84 ellipsoid.
m = bm(projection='cyl', llcrnrlon=30, llcrnrlat=-40,urcrnrlon=120, urcrnrlat=30)
# transform to nx x ny regularly spaced 5km native projection grid
nx = int((m.xmax-m.xmin))+1; ny = int((m.ymax-m.ymin)+1)
q=ru[1:2190]
qmean=np.mean(q)
qstd=np.std(q)
qnorm=(q-qmean)/qstd
lags3d=np.arange(731*140*180).reshape(731,140,180)
r3d=np.arange(731*140*180).reshape(731,140,180)
for i in np.arange(len(lons)):
for j in np.arange(len(lats)):
print i,j
p=qu[1:2190,j,i].squeeze()
p.shape
pmean=np.mean(p)
pstd=np.std(p)
pnorm=(p-pmean)/pstd
n=len(p)
# fg=plt.figure()
c=plt.xcorr(p,q,usevlines=True,maxlags=nlags,normed=True,lw=2)
acp=plt.acorr(p,usevlines=True,maxlags=nlags,normed=True,lw=2)
acq=plt.acorr(q,usevlines=True,maxlags=nlags,normed=True,lw=2)
acp[1][nlags]=0
acq[1][nlags]=0
lags=c[0]
r=c[1]
lags3d[:,j,i]=lags
r3d[:,j,i]=r