-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathssh_mean_shf.py
111 lines (88 loc) · 3.2 KB
/
ssh_mean_shf.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
107
108
109
110
111
# coding: utf-8
# In[1]:
import os
import numpy as np
import pandas as pd
import xarray as xr
import glob
import matplotlib.pyplot as plt
#get_ipython().magic(u'matplotlib inline')
#from matplotlib import animation
#import cartopy.crs as ccrs
from scipy.io import netcdf
from scipy import signal
from netCDF4 import num2date
from netCDF4 import Dataset as NetCDFFile
import matplotlib.dates as mdates
#from mpl_toolkits.basemap import Basemap
import datetime
#from wavelets import WaveletAnalysis
# In[2]:
path='/ccc/scratch/cont003/gen0727/garciagi/'
#mfiles= 'GOM025-GSL301.050_y1993-2012.1d_gridT.nc','GOM025-GSL301.001_y1993-2012.1d_gridT.nc'
#filenames=sorted(glob.glob('/home/users/garcia8ix/workd/ENTROPY/5-days-avg/*.nc'))
#print(len(filenames))
#meanfile='/home/users/garcia8ix/workd/ENTROPY/stats/Smean.nc'
# In[3]:
mfiles=sorted(glob.glob('/ccc/scratch/cont003/gen0727/garciagi/FILTER/TEMP/GOM025.GSL301_m*.nc'))
trfiles=sorted(glob.glob('/ccc/scratch/cont003/gen0727/garciagi/outtrends/test/splinesH*'))
#mfiles=mfiles[0:25]
#dsh = xr.open_mfdataset(mfiles,concat_dim='ensemble')
# In[4]:
dsl=xr.open_dataset('/ccc/scratch/cont003/gen0727/garciagi/1d/NATL025.GSL301_m024.1d_gridT.nc')
#mds=xr.open_dataset(meanfile)
lats=dsl.nav_lat
lons=dsl.nav_lon
time2=dsl.time_centered
time3 = pd.date_range('1993-01-01','2012-12-26' , freq='5D')
ifens=0
ntime=len(time2[:])
nrow=lats.shape[0]
ncols=lats.shape[1]
print(len(mfiles))
meanS= NetCDFFile('/ccc/scratch/cont003/gen0727/garciagi/SSHmeanHF.nc', 'w', format='NETCDF3_64BIT')
meanS.createDimension('time_counter', ntime)
meanS.createDimension('y', nrow)
meanS.createDimension('x', ncols)
meanS.createDimension('member',len(mfiles))
time = meanS.createVariable('time', 'f4', ('time_counter',))
lats = meanS.createVariable('lons','f4',('y','x'))
lons = meanS.createVariable('lats','f4',('y','x'))
meanssh = meanS.createVariable('meanssh', 'f4', ('time_counter', 'y','x'))
varssh = meanS.createVariable('varssh', 'f4', ('time_counter', 'y','x'))
varssh2 = meanS.createVariable('varssh2', 'f4', ('time_counter', 'y','x'))
varsshi = meanS.createVariable('varsshi', 'f4', ('member', 'y','x'))
filread = meanS.createVariable('filread', 'f4', ('member'))
mlo=[1]
stm=np.zeros((len(time2[:]),lats.shape[0],lats.shape[1]))
stdm=np.zeros((len(time2[:]),lats.shape[0],lats.shape[1]))
sti=np.zeros((len(mfiles),lats.shape[0],lats.shape[1]))
fils=np.zeros(len(mfiles))
for item in mfiles:
mean_file=mfiles[ifens]
print(mean_file)
#lines=np.load(trfiles[ifens])
dsem = xr.open_dataset(mean_file)
for item in mlo:
lo=50#mlo[idt]
la=90#mla[idt]
xh=dsem.ssht[:,:,:]
print(xh.shape)
stm=stm+xh
stdm=stdm+(xh*xh)
sti[ifens,:,:]=np.var(xh,axis=0)
print('ensemble %s' %ifens)
#dm["meanE{0}".format(inm)]=stm
#ds["stdE{0}".format(inm)]=stdm
fils[ifens]=ifens
ifens=ifens+1
dsem=None
time[:]=dsl.time_centered[:].data
lons[:,:]=dsl.nav_lon.data
lats[:,:]=dsl.nav_lat.data
meanssh[:,:,:]=stm.data/50.
varssh[:,:,:] =stdm.data/50.
varssh2[:,:,:] =(stdm.data/50.)-((stm.data/50.)*(stm.data/50.))
varsshi[:,:,:]=sti
filread[:]=fils
meanS.close()