source: trunk/UTIL/PYTHON/mcs.py @ 464

Last change on this file since 464 was 464, checked in by acolaitis, 13 years ago

Changed name of the gcm utilities wrapper to gcm_transformations.py. Added hrecast to it, though there are no options yet in pp.py. hrecast will be used in mcs.py and in a future tes.py. It can be implemented in pp.py easily. Added --latreverse option to mcs.py to output files with a south_north latitude axis, as in the data.

  • Property svn:executable set to *
File size: 12.8 KB
Line 
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)
14if __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")
Note: See TracBrowser for help on using the repository browser.