1 | #!/usr/bin/env python |
---|
2 | |
---|
3 | ### A. Colaitis |
---|
4 | |
---|
5 | ## |
---|
6 | # This routine transforms a diagfi.nc file into a diagfi_MCS.nc file where |
---|
7 | # the fields are directly comparable to those contained in MCS data, i.e. |
---|
8 | # fields are re-binned at times over the ranges specified in the MCS file. |
---|
9 | ### |
---|
10 | |
---|
11 | ########################################################################################### |
---|
12 | ########################################################################################### |
---|
13 | ### What is below relate to running the file as a command line executable (very convenient) |
---|
14 | if __name__ == "__main__": |
---|
15 | import sys |
---|
16 | from optparse import OptionParser ### to be replaced by argparse |
---|
17 | from netCDF4 import Dataset |
---|
18 | from os import system,path |
---|
19 | from times import sol2ls |
---|
20 | import numpy as np |
---|
21 | from mymath import find_nearest |
---|
22 | from myplot import getfield,separatenames |
---|
23 | from make_netcdf import make_gcm_netcdf |
---|
24 | from gcm_transformations import call_zrecast,call_hrecast |
---|
25 | parser = OptionParser() |
---|
26 | |
---|
27 | ############################# |
---|
28 | ### Options |
---|
29 | parser.add_option('-f', '--file', action='store',dest='file', type="string", default=None, help='[NEEDED] filename.') |
---|
30 | parser.add_option('-m', '--mfile', action='store',dest='mcsfile', type="string", default=None, help='[NEEDED] filename for MCS comparison.') |
---|
31 | parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='[NEEDED] Variables to process. (coma-separated list. aps and bps are always included.)') |
---|
32 | parser.add_option('-x', action='store_true',dest='recast', default=False, help='Force aps and bps to be ommited in output file (usefull if your file is already recasted along z) [False]') |
---|
33 | parser.add_option('-i', '--zrecast', action='store_true', dest='zrecast', default=False, help='Cast zrecast.e on diagfi file with MCS pressure levels. Will pass this operation is recasted file is already present, unless --override is specified. [False]') |
---|
34 | parser.add_option('--override', action='store_true', dest='override', default=False, help='Force zrecast.e to act even if recasted file is already present(will erase previous recasted file) [False]') |
---|
35 | parser.add_option('--ditch', action='store_true', dest='ditch', default=False, help='Ditch recasted file when interpolation is performed. [False]') |
---|
36 | parser.add_option('--latreverse', action='store_true', dest='latreverse', default=False, help='Reverse the latitude axis in output diagfi. [False]') |
---|
37 | |
---|
38 | ############################# |
---|
39 | ### Get options and variables |
---|
40 | (opt,args) = parser.parse_args() |
---|
41 | |
---|
42 | ############################# |
---|
43 | ### Load and check data |
---|
44 | |
---|
45 | if opt.var is None: |
---|
46 | print "You must specify at least a field to process with -v." |
---|
47 | exit() |
---|
48 | |
---|
49 | # Zrecast |
---|
50 | |
---|
51 | varznames=separatenames(opt.var[0]) |
---|
52 | |
---|
53 | if opt.zrecast: |
---|
54 | if (path.exists(opt.file[0:len(opt.file)-3]+"_P.nc") and (not opt.override)): |
---|
55 | print "--> "+opt.file[0:len(opt.file)-3]+"_P.nc" |
---|
56 | print "Recasted file is already there, skipping interpolation. [use --override to force interpolation]" |
---|
57 | filename=opt.file[0:len(opt.file)-3]+"_P.nc" |
---|
58 | else: |
---|
59 | print "--> "+opt.file[0:len(opt.file)-3]+"_P.nc" |
---|
60 | filename=call_zrecast ( interp_mode = 2, \ |
---|
61 | input_name = [opt.file], \ |
---|
62 | fields = varznames, \ |
---|
63 | predifined = 'mcs')[0] |
---|
64 | else:filename=opt.file |
---|
65 | # Files |
---|
66 | |
---|
67 | print "--> Loading diagfi dataset." |
---|
68 | |
---|
69 | nc=Dataset(filename) |
---|
70 | ncmcs=Dataset(opt.mcsfile) |
---|
71 | |
---|
72 | # Dimensions |
---|
73 | |
---|
74 | lon=nc.variables["longitude"][:] |
---|
75 | lat=nc.variables["latitude"][:] |
---|
76 | alt=nc.variables["altitude"][:] |
---|
77 | time=nc.variables["Time"][:] # in fraction of sols |
---|
78 | if "controle" in nc.variables: |
---|
79 | controle=nc.variables["controle"][:] |
---|
80 | day_ini=controle[3]%669 |
---|
81 | else: |
---|
82 | if opt.zrecast: |
---|
83 | nccontrol=Dataset(opt.file) |
---|
84 | if "controle" in nccontrol.variables: |
---|
85 | controle=nccontrol.variables["controle"][:] |
---|
86 | day_ini=controle[3]%669 |
---|
87 | else: |
---|
88 | print "Error: could not find controle variable in diagfi." |
---|
89 | day_ini=input("Please type initial sol number:")%669 |
---|
90 | else: |
---|
91 | print "Error: could not find controle variable in diagfi." |
---|
92 | day_ini=input("Please type initial sol number:")%669 |
---|
93 | time[:]=time[:]+day_ini |
---|
94 | nx=len(lon) |
---|
95 | ny=len(lat) |
---|
96 | nz=len(alt) |
---|
97 | nt=len(time) |
---|
98 | lstime=sol2ls(time) |
---|
99 | |
---|
100 | # MCS |
---|
101 | |
---|
102 | print "--> Loading and preparing MCS dataset." |
---|
103 | |
---|
104 | dtimemintmp=ncmcs.variables["dtimemin"][:,:,:] |
---|
105 | dtimemaxtmp=ncmcs.variables["dtimemax"][:,:,:] |
---|
106 | ntimemintmp=ncmcs.variables["ntimemin"][:,:,:] |
---|
107 | ntimemaxtmp=ncmcs.variables["ntimemax"][:,:,:] |
---|
108 | lonmcs=ncmcs.variables["longitude"][:] |
---|
109 | latmcs=ncmcs.variables["latitude"][:] |
---|
110 | timemcs=ncmcs.variables["time"][:]%360 # IN LS |
---|
111 | |
---|
112 | dtimemin=np.ma.masked_where(dtimemintmp < 0.,dtimemintmp) |
---|
113 | dtimemin.set_fill_value([np.NaN]) |
---|
114 | dtimemax=np.ma.masked_where(dtimemaxtmp < 0.,dtimemaxtmp) |
---|
115 | dtimemax.set_fill_value([np.NaN]) |
---|
116 | ntimemin=np.ma.masked_where(ntimemintmp < 0.,ntimemintmp) |
---|
117 | ntimemin.set_fill_value([np.NaN]) |
---|
118 | ntimemax=np.ma.masked_where(ntimemaxtmp < 0.,ntimemaxtmp) |
---|
119 | ntimemax.set_fill_value([np.NaN]) |
---|
120 | |
---|
121 | # Variables to treat |
---|
122 | |
---|
123 | print "--> Preparing diagfi dataset." |
---|
124 | |
---|
125 | varz=[] |
---|
126 | n=0 |
---|
127 | for zn in varznames: |
---|
128 | load=getfield(nc,zn) |
---|
129 | load=np.ma.masked_where(load < -1.e-20,load) |
---|
130 | load.set_fill_value([np.NaN]) |
---|
131 | load=load.filled() |
---|
132 | load=np.ma.masked_invalid(load) |
---|
133 | load.set_fill_value([np.NaN]) |
---|
134 | load=load.filled() |
---|
135 | varz.append(load) |
---|
136 | load=0. |
---|
137 | print "Found: "+zn+" with dimensions: " |
---|
138 | print np.array(varz[n]).shape |
---|
139 | n=n+1 |
---|
140 | |
---|
141 | nzvar=len(varz) |
---|
142 | dimensions={} |
---|
143 | vv=0 |
---|
144 | for var in varz: |
---|
145 | a=len(np.array(var).shape) |
---|
146 | if a == 3: dimensions[vv]=a |
---|
147 | elif a == 4: dimensions[vv]=a |
---|
148 | else: |
---|
149 | print "Warning, only 3d and 4d variables are supported for time-recasting" |
---|
150 | exit() |
---|
151 | vv=vv+1 |
---|
152 | |
---|
153 | # Variables to save but not treated (only along z, or phisinit-like files) |
---|
154 | |
---|
155 | aps=nc.variables["aps"][:] |
---|
156 | bps=nc.variables["bps"][:] |
---|
157 | fullnames=["aps","bps"] |
---|
158 | for name in varznames: |
---|
159 | fullnames.append("d"+name) |
---|
160 | fullnames.append("n"+name) |
---|
161 | print "Will output: " |
---|
162 | if opt.recast: print fullnames[2:] |
---|
163 | else: print fullnames |
---|
164 | ############################# |
---|
165 | ### Building |
---|
166 | ############################# |
---|
167 | |
---|
168 | ### We loop over chunks of gcm data corresponding to MCS time dimension |
---|
169 | ### Bin sizes for mcs data is 5 degrees ls centered on value |
---|
170 | varday=np.zeros([len(timemcs),nz,ny,nx]) |
---|
171 | varnight=np.zeros([len(timemcs),nz,ny,nx]) |
---|
172 | vardayout=np.zeros([nzvar,len(timemcs),nz,ny,nx]) |
---|
173 | varnightout=np.zeros([nzvar,len(timemcs),nz,ny,nx]) |
---|
174 | vardayout=np.ma.masked_invalid(vardayout) |
---|
175 | varnightout=np.ma.masked_invalid(varnightout) |
---|
176 | i=0 |
---|
177 | for ls in timemcs: |
---|
178 | lsstart=ls-2.5 |
---|
179 | lsstop=ls+2.5 |
---|
180 | istart=find_nearest(lstime,lsstart,strict=True) |
---|
181 | istop=find_nearest(lstime,lsstop,strict=True) |
---|
182 | varchk=0 |
---|
183 | if ((istart is np.NaN) or (istop is np.NaN)): |
---|
184 | vardayout[:,i,:,:,:]=np.NaN |
---|
185 | varnightout[:,i,:,:,:]=np.NaN |
---|
186 | print "Time interval skipped. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' |
---|
187 | i=i+1 |
---|
188 | continue |
---|
189 | print "--->> Processing Data. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' |
---|
190 | # warning, python's convention is that the second index of array[a:b] is the array index of element after the one being picked last. |
---|
191 | # for that reason, array[0:0] is nan and array[0:1] is only one value. Hence, len(array[a:b+1]) is b-a+1 and not b-a+2 |
---|
192 | print " .initialisation." |
---|
193 | |
---|
194 | |
---|
195 | varchk=np.zeros([nzvar,istop-istart+1,nz,ny,nx]) |
---|
196 | vv=0 |
---|
197 | for variable in varz: |
---|
198 | if dimensions[vv] is 3: |
---|
199 | varchk[vv,:,0,:,:]=variable[istart:istop+1,:,:] |
---|
200 | else: |
---|
201 | varchk[vv,:,:,:,:]=variable[istart:istop+1,:,:,:] |
---|
202 | vv=vv+1 |
---|
203 | varchk=np.ma.masked_invalid(varchk) |
---|
204 | varchk.set_fill_value([np.NaN]) |
---|
205 | varchktime=time[istart:istop+1] |
---|
206 | ndays=np.floor(varchktime[len(varchktime)-1])-np.floor(varchktime[0]) |
---|
207 | dtmichk=dtimemin[i,:,:] |
---|
208 | dtmachk=dtimemax[i,:,:] |
---|
209 | ntmichk=ntimemin[i,:,:] |
---|
210 | ntmachk=ntimemax[i,:,:] |
---|
211 | dtmichk.set_fill_value([np.NaN]) |
---|
212 | dtmachk.set_fill_value([np.NaN]) |
---|
213 | ntmichk.set_fill_value([np.NaN]) |
---|
214 | ntmachk.set_fill_value([np.NaN]) |
---|
215 | dtmichk=dtmichk.filled() |
---|
216 | dtmachk=dtmachk.filled() |
---|
217 | ntmichk=ntmichk.filled() |
---|
218 | ntmachk=ntmachk.filled() |
---|
219 | |
---|
220 | ### We iterate for each day in the chunk, on each grid point we find |
---|
221 | ### the closest corresponding MCS grid point and the index of the |
---|
222 | ### time in the chunk closest to the time in the closest MCS grid point. |
---|
223 | ### (yea it's complicated) |
---|
224 | |
---|
225 | vartmpnight=np.zeros([nzvar,ndays,nz,ny,nx]) |
---|
226 | vartmpday=np.zeros([nzvar,ndays,nz,ny,nx]) |
---|
227 | vartmpnight=np.ma.masked_invalid(vartmpnight) |
---|
228 | vartmpday=np.ma.masked_invalid(vartmpday) |
---|
229 | vartmpnight.set_fill_value([np.NaN]) |
---|
230 | vartmpday.set_fill_value([np.NaN]) |
---|
231 | |
---|
232 | nd=0 |
---|
233 | print " .time indices MCS grid -> diagfi grid." |
---|
234 | while nd < ndays: |
---|
235 | |
---|
236 | daystart=find_nearest(varchktime-varchktime[0],nd) |
---|
237 | daystop=find_nearest(varchktime-varchktime[0],nd+1) |
---|
238 | # varchktime_lon=np.zeros([daystop-daystart+1,len(lon)]) |
---|
239 | ix=0 |
---|
240 | for x in lon: |
---|
241 | |
---|
242 | varchktime_lon = 24.*(varchktime[daystart:daystop+1]-varchktime[daystart]) + x/15. |
---|
243 | |
---|
244 | iy=0 |
---|
245 | for y in lat: |
---|
246 | niy=find_nearest(latmcs,y) |
---|
247 | nix=find_nearest(lonmcs,x) |
---|
248 | nitdtmichk=find_nearest(varchktime_lon,dtmichk[niy,nix]) |
---|
249 | nitdtmachk=find_nearest(varchktime_lon,dtmachk[niy,nix]) |
---|
250 | nitntmichk=find_nearest(varchktime_lon,ntmichk[niy,nix]) |
---|
251 | nitntmachk=find_nearest(varchktime_lon,ntmachk[niy,nix]) |
---|
252 | for vv in np.arange(nzvar): |
---|
253 | if ((nitdtmichk is np.NaN) or (nitdtmachk is np.NaN)): |
---|
254 | vartmpday[vv,nd,:,iy,ix]=np.NaN |
---|
255 | elif nitdtmichk > nitdtmachk: |
---|
256 | vartmpday[vv,nd,:,iy,ix]=(np.ma.mean(varchk[vv,daystart+nitdtmichk:daystop+1,:,iy,ix],axis=0)+np.ma.mean(varchk[vv,daystart:daystart+nitdtmachk+1,:,iy,ix],axis=0))/2. |
---|
257 | else: |
---|
258 | vartmpday[vv,nd,:,iy,ix]=np.ma.mean(varchk[vv,daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0) |
---|
259 | if ((nitntmichk is np.NaN) or (nitntmachk is np.NaN)): |
---|
260 | vartmpnight[vv,nd,:,iy,ix]=np.NaN |
---|
261 | elif nitntmichk > nitntmachk: |
---|
262 | vartmpnight[vv,nd,:,iy,ix]=(np.ma.mean(varchk[vv,daystart+nitntmichk:daystop+1,:,iy,ix],axis=0)+np.ma.mean(varchk[vv,daystart:daystart+nitntmachk+1,:,iy,ix],axis=0))/2. |
---|
263 | else: |
---|
264 | vartmpnight[vv,nd,:,iy,ix]=np.ma.mean(varchk[vv,daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0) |
---|
265 | iy=iy+1 |
---|
266 | ix=ix+1 |
---|
267 | nd=nd+1 |
---|
268 | |
---|
269 | print " .creating bins." |
---|
270 | |
---|
271 | vartmpdaymasked=np.ma.masked_invalid(vartmpday) |
---|
272 | vartmpnightmasked=np.ma.masked_invalid(vartmpnight) |
---|
273 | vartmpdaymasked.set_fill_value([np.NaN]) |
---|
274 | vartmpnightmasked.set_fill_value([np.NaN]) |
---|
275 | for vv in np.arange(nzvar): |
---|
276 | vardayout[vv,i,:,:,:]=np.ma.mean(vartmpdaymasked[vv,:,:,:,:],axis=0) |
---|
277 | varnightout[vv,i,:,:,:]=np.ma.mean(vartmpnightmasked[vv,:,:,:,:],axis=0) |
---|
278 | print " ."+varznames[vv]+".done" |
---|
279 | i=i+1 |
---|
280 | |
---|
281 | print "--->> Preparing Data for ncdf. Missing value is NaN." |
---|
282 | |
---|
283 | vardayout=np.ma.masked_invalid(vardayout) |
---|
284 | varnightout=np.ma.masked_invalid(varnightout) |
---|
285 | vardayout.set_fill_value([np.NaN]) |
---|
286 | varnightout.set_fill_value([np.NaN]) |
---|
287 | |
---|
288 | if opt.latreverse: |
---|
289 | vardayout[:,:,:,:,:]=vardayout[:,:,:,::-1,:] |
---|
290 | varnightout[:,:,:,:,:]=varnightout[:,:,:,::-1,:] |
---|
291 | |
---|
292 | all=[aps,bps] |
---|
293 | for vv in np.arange(nzvar): |
---|
294 | if dimensions[vv] == 3: |
---|
295 | all.append(vardayout[vv,:,0,:,:].filled()) |
---|
296 | all.append(varnightout[vv,:,0,:,:].filled()) |
---|
297 | elif dimensions[vv] == 4: |
---|
298 | all.append(vardayout[vv,:,:,:,:].filled()) |
---|
299 | all.append(varnightout[vv,:,:,:,:].filled()) |
---|
300 | |
---|
301 | if opt.recast: |
---|
302 | all=all[2:] |
---|
303 | fullnames=fullnames[2:] |
---|
304 | |
---|
305 | if opt.latreverse: |
---|
306 | lat=lat[::-1] |
---|
307 | |
---|
308 | make_gcm_netcdf (zfilename="diagfi_MCS.nc", \ |
---|
309 | zdescription="Temperatures from diagfi reworked to match MCS format", \ |
---|
310 | zlon=lon, \ |
---|
311 | zlat=lat, \ |
---|
312 | zalt=alt, \ |
---|
313 | ztime=timemcs, \ |
---|
314 | zvariables=all, \ |
---|
315 | znames=fullnames) |
---|
316 | if opt.zrecast and opt.ditch: |
---|
317 | print "removing interpolated file" |
---|
318 | system("rm -f "+opt.file[0:len(opt.file)-3]+"_P.nc") |
---|