-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathssh_demap.py
executable file
·106 lines (98 loc) · 3.36 KB
/
ssh_demap.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
import os
import numpy as np
import pandas as pd
import xarray as xr
import glob
import sys
import statsmodels.api as sm
lowess = sm.nonparametric.lowess
import scipy.interpolate as scpi
import time
from netCDF4 import num2date
from netCDF4 import Dataset as NetCDFFile
pathi='/ccc/scratch/cont003/gen0727/garciagi/'
patho='/ccc/scratch/cont003/gen0727/garciagi/outtrends/'
mfiles=sorted(glob.glob(pathi+'1d/*.nc'))
trhfiles=sorted(glob.glob(pathi+'outtrends/SPLINES/splinesH*'))
trtfiles=sorted(glob.glob(pathi+'outtrends/SPLINES/splinesT*'))
print(len(mfiles),len(trtfiles),len(trhfiles))
mem=sys.argv[1]
membs=[int(mem)]
for im in range(len(membs)):
ifens=membs[im]
idmem=ifens+1
#print(ifens)
mean_file=mfiles[ifens]
print(mean_file)
dsem = xr.open_dataset(mean_file)
linest=np.load(trtfiles[ifens])
linesh=np.load(trhfiles[ifens])
print(trtfiles[ifens],trhfiles[ifens])
time2=dsem.time_centered
timet = pd.date_range('1993-01-01','2012-12-26' , freq='D')
#print(time2.shape,timet.shape,timet2.shape)
nt,ni,nj=dsem.ssh.shape
####Declare the NETCDF file
meanM= NetCDFFile(patho+'outdetrend/GOM025.GSL301_m0%s_1d_dtrend.nc' %(idmem), 'w', format='NETCDF3_64BIT')
meanM.createDimension('time_counter', nt)
meanM.createDimension('y', ni)
meanM.createDimension('x', nj)
#time_center = meanM.createVariable('time_center', 'f4', ('time_counter',))
lons = meanM.createVariable('lons','f4',('y','x'))
lats = meanM.createVariable('lats','f4',('y','x'))
sshd = meanM.createVariable('sshd', 'f4', ('time_counter', 'y','x'))
sstd = meanM.createVariable('sstd', 'f4', ('time_counter', 'y','x'))
timet2 = pd.date_range('1993-01-01','2012-12-26' , freq='D')
nt,ni,nj=dsem.ssh.shape
print(nt,ni,nj)
nt2=nt-(292*5+1)
sshtr=np.ones([nt,ni,nj])
sshtr[:]=np.NaN
ssttr=np.ones([nt,ni,nj])
ssttr[:]=np.NaN
#tremh=np.ones([nt,ni,nj])
#start1=time.clock()
for ii in range(ni):
print "Start ii:", ii
for jj in range(nj):
#print "jj:", jj
xh=dsem.ssh[:,ii,jj]
xt=dsem.sst[:,ii,jj]
#treh=linesh[ii,jj]
#tret=linest[ii,jj]
#xh=dsem.ssh[292*5+1:,ii,jj]
#xt=dsem.sst[:,ii,jj]
start2=time.clock()
if np.isnan(xh[100]):
tremh=np.ones(len(timet2))
tremh[:]=np.NaN
tremt=np.ones(len(timet2))
tremt[:]=np.NaN
pass
else:
#start2=time.clock()
#shm=xh.resample('M',dim='time_counter',how='mean')
#print(time.clock()-start2
#print(xt.shape,xh.shape,timet2.shape)
treh=linesh[ii,jj]
#print(ii,jj,'treh ',time.clock()-start2)
tret=linest[ii,jj]
#print('tret',time.clock()-start2)
tremh=treh(timet2.to_julian_date())
#print('trem',time.clock()-start2)
tremt=tret(timet2.to_julian_date())
#print('tremt',time.clock()-start2)
#start2=time.clock()
sshtr[:,ii,jj]=xh-tremh
#print(ii,jj,'sshtr',time.clock()-start2)
ssttr[:,ii,jj]=xt-tremt
#print('sstr',time.clock()-start2)
#print(time.clock()-start1)
#print(ii)
lons[:,:]=dsem.nav_lon.data
lats[:,:]=dsem.nav_lat.data
#time_center[:]=time2
sshd[:,:,:]=sshtr
sstd[:,:,:]=ssttr
#print(ifens)
meanM.close()