| 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 zrecast_wrapper import call_zrecast |
|---|
| 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 | ############################# |
|---|
| 37 | ### Get options and variables |
|---|
| 38 | (opt,args) = parser.parse_args() |
|---|
| 39 | |
|---|
| 40 | ############################# |
|---|
| 41 | ### Load and check data |
|---|
| 42 | |
|---|
| 43 | if opt.var is None: |
|---|
| 44 | print "You must specify at least a field to process with -v." |
|---|
| 45 | exit() |
|---|
| 46 | |
|---|
| 47 | # Zrecast |
|---|
| 48 | |
|---|
| 49 | varznames=separatenames(opt.var[0]) |
|---|
| 50 | |
|---|
| 51 | if opt.zrecast: |
|---|
| 52 | if (path.exists(opt.file[0:len(opt.file)-3]+"_P.nc") and (not opt.override)): |
|---|
| 53 | print "--> "+opt.file[0:len(opt.file)-3]+"_P.nc" |
|---|
| 54 | print "Recasted file is already there, skipping interpolation. [use --override to force interpolation]" |
|---|
| 55 | filename=opt.file[0:len(opt.file)-3]+"_P.nc" |
|---|
| 56 | else: |
|---|
| 57 | print "--> "+opt.file[0:len(opt.file)-3]+"_P.nc" |
|---|
| 58 | filename=call_zrecast ( interp_mode = 2, \ |
|---|
| 59 | input_name = [opt.file], \ |
|---|
| 60 | fields = varznames, \ |
|---|
| 61 | predifined = 'mcs')[0] |
|---|
| 62 | else:filename=opt.file |
|---|
| 63 | # Files |
|---|
| 64 | |
|---|
| 65 | print "--> Loading diagfi dataset." |
|---|
| 66 | |
|---|
| 67 | nc=Dataset(filename) |
|---|
| 68 | ncmcs=Dataset(opt.mcsfile) |
|---|
| 69 | |
|---|
| 70 | # Dimensions |
|---|
| 71 | |
|---|
| 72 | lon=nc.variables["longitude"][:] |
|---|
| 73 | lat=nc.variables["latitude"][:] |
|---|
| 74 | alt=nc.variables["altitude"][:] |
|---|
| 75 | time=nc.variables["Time"][:] # in fraction of sols |
|---|
| 76 | if "controle" in nc.variables: |
|---|
| 77 | controle=nc.variables["controle"][:] |
|---|
| 78 | day_ini=controle[3]%669 |
|---|
| 79 | else: |
|---|
| 80 | if opt.zrecast: |
|---|
| 81 | nccontrol=Dataset(opt.file) |
|---|
| 82 | if "controle" in nccontrol.variables: |
|---|
| 83 | controle=nccontrol.variables["controle"][:] |
|---|
| 84 | day_ini=controle[3]%669 |
|---|
| 85 | else: |
|---|
| 86 | print "Error: could not find controle variable in diagfi." |
|---|
| 87 | day_ini=input("Please type initial sol number:")%669 |
|---|
| 88 | time[:]=time[:]+day_ini |
|---|
| 89 | nx=len(lon) |
|---|
| 90 | ny=len(lat) |
|---|
| 91 | nz=len(alt) |
|---|
| 92 | nt=len(time) |
|---|
| 93 | lstime=sol2ls(time) |
|---|
| 94 | |
|---|
| 95 | # MCS |
|---|
| 96 | |
|---|
| 97 | print "--> Loading and preparing MCS dataset." |
|---|
| 98 | |
|---|
| 99 | dtimemintmp=ncmcs.variables["dtimemin"][:,:,:] |
|---|
| 100 | dtimemaxtmp=ncmcs.variables["dtimemax"][:,:,:] |
|---|
| 101 | ntimemintmp=ncmcs.variables["ntimemin"][:,:,:] |
|---|
| 102 | ntimemaxtmp=ncmcs.variables["ntimemax"][:,:,:] |
|---|
| 103 | lonmcs=ncmcs.variables["longitude"][:] |
|---|
| 104 | latmcs=ncmcs.variables["latitude"][:] |
|---|
| 105 | timemcs=ncmcs.variables["time"][:]%360 # IN LS |
|---|
| 106 | |
|---|
| 107 | dtimemin=np.ma.masked_where(dtimemintmp < 0.,dtimemintmp) |
|---|
| 108 | dtimemin.set_fill_value([np.NaN]) |
|---|
| 109 | dtimemax=np.ma.masked_where(dtimemaxtmp < 0.,dtimemaxtmp) |
|---|
| 110 | dtimemax.set_fill_value([np.NaN]) |
|---|
| 111 | ntimemin=np.ma.masked_where(ntimemintmp < 0.,ntimemintmp) |
|---|
| 112 | ntimemin.set_fill_value([np.NaN]) |
|---|
| 113 | ntimemax=np.ma.masked_where(ntimemaxtmp < 0.,ntimemaxtmp) |
|---|
| 114 | ntimemax.set_fill_value([np.NaN]) |
|---|
| 115 | |
|---|
| 116 | # Variables to treat |
|---|
| 117 | |
|---|
| 118 | print "--> Preparing diagfi dataset." |
|---|
| 119 | |
|---|
| 120 | varz=[] |
|---|
| 121 | n=0 |
|---|
| 122 | for zn in varznames: |
|---|
| 123 | load=getfield(nc,zn) |
|---|
| 124 | load=np.ma.masked_where(load < -1.e-20,load) |
|---|
| 125 | load.set_fill_value([np.NaN]) |
|---|
| 126 | load=load.filled() |
|---|
| 127 | load=np.ma.masked_invalid(load) |
|---|
| 128 | load.set_fill_value([np.NaN]) |
|---|
| 129 | load=load.filled() |
|---|
| 130 | varz.append(load) |
|---|
| 131 | load=0. |
|---|
| 132 | print "Found: "+zn+" with dimensions: " |
|---|
| 133 | print np.array(varz[n]).shape |
|---|
| 134 | n=n+1 |
|---|
| 135 | |
|---|
| 136 | nzvar=len(varz) |
|---|
| 137 | dimensions={} |
|---|
| 138 | vv=0 |
|---|
| 139 | for var in varz: |
|---|
| 140 | a=len(np.array(var).shape) |
|---|
| 141 | if a == 3: dimensions[vv]=a |
|---|
| 142 | elif a == 4: dimensions[vv]=a |
|---|
| 143 | else: |
|---|
| 144 | print "Warning, only 3d and 4d variables are supported for time-recasting" |
|---|
| 145 | exit() |
|---|
| 146 | vv=vv+1 |
|---|
| 147 | |
|---|
| 148 | # Variables to save but not treated (only along z, or phisinit-like files) |
|---|
| 149 | |
|---|
| 150 | aps=nc.variables["aps"][:] |
|---|
| 151 | bps=nc.variables["bps"][:] |
|---|
| 152 | fullnames=["aps","bps"] |
|---|
| 153 | for name in varznames: |
|---|
| 154 | fullnames.append("d"+name) |
|---|
| 155 | fullnames.append("n"+name) |
|---|
| 156 | print "Will output: " |
|---|
| 157 | if opt.recast: print fullnames[2:] |
|---|
| 158 | else: print fullnames |
|---|
| 159 | ############################# |
|---|
| 160 | ### Building |
|---|
| 161 | ############################# |
|---|
| 162 | |
|---|
| 163 | ### We loop over chunks of gcm data corresponding to MCS time dimension |
|---|
| 164 | ### Bin sizes for mcs data is 5 degrees ls centered on value |
|---|
| 165 | varday=np.zeros([len(timemcs),nz,ny,nx]) |
|---|
| 166 | varnight=np.zeros([len(timemcs),nz,ny,nx]) |
|---|
| 167 | vardayout=np.zeros([nzvar,len(timemcs),nz,ny,nx]) |
|---|
| 168 | varnightout=np.zeros([nzvar,len(timemcs),nz,ny,nx]) |
|---|
| 169 | vardayout=np.ma.masked_invalid(vardayout) |
|---|
| 170 | varnightout=np.ma.masked_invalid(varnightout) |
|---|
| 171 | i=0 |
|---|
| 172 | for ls in timemcs: |
|---|
| 173 | lsstart=ls-2.5 |
|---|
| 174 | lsstop=ls+2.5 |
|---|
| 175 | istart=find_nearest(lstime,lsstart,strict=True) |
|---|
| 176 | istop=find_nearest(lstime,lsstop,strict=True) |
|---|
| 177 | varchk=0 |
|---|
| 178 | if ((istart is np.NaN) or (istop is np.NaN)): |
|---|
| 179 | vardayout[:,i,:,:,:]=np.NaN |
|---|
| 180 | varnightout[:,i,:,:,:]=np.NaN |
|---|
| 181 | print "Time interval skipped. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' |
|---|
| 182 | i=i+1 |
|---|
| 183 | continue |
|---|
| 184 | print "--->> Processing Data. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' |
|---|
| 185 | # 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. |
|---|
| 186 | # 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 |
|---|
| 187 | print " .initialisation." |
|---|
| 188 | |
|---|
| 189 | |
|---|
| 190 | varchk=np.zeros([nzvar,istop-istart+1,nz,ny,nx]) |
|---|
| 191 | vv=0 |
|---|
| 192 | for variable in varz: |
|---|
| 193 | if dimensions[vv] is 3: |
|---|
| 194 | varchk[vv,:,0,:,:]=variable[istart:istop+1,:,:] |
|---|
| 195 | else: |
|---|
| 196 | varchk[vv,:,:,:,:]=variable[istart:istop+1,:,:,:] |
|---|
| 197 | vv=vv+1 |
|---|
| 198 | varchk=np.ma.masked_invalid(varchk) |
|---|
| 199 | varchk.set_fill_value([np.NaN]) |
|---|
| 200 | varchktime=time[istart:istop+1] |
|---|
| 201 | ndays=np.floor(varchktime[len(varchktime)-1])-np.floor(varchktime[0]) |
|---|
| 202 | dtmichk=dtimemin[i,:,:] |
|---|
| 203 | dtmachk=dtimemax[i,:,:] |
|---|
| 204 | ntmichk=ntimemin[i,:,:] |
|---|
| 205 | ntmachk=ntimemax[i,:,:] |
|---|
| 206 | dtmichk.set_fill_value([np.NaN]) |
|---|
| 207 | dtmachk.set_fill_value([np.NaN]) |
|---|
| 208 | ntmichk.set_fill_value([np.NaN]) |
|---|
| 209 | ntmachk.set_fill_value([np.NaN]) |
|---|
| 210 | dtmichk=dtmichk.filled() |
|---|
| 211 | dtmachk=dtmachk.filled() |
|---|
| 212 | ntmichk=ntmichk.filled() |
|---|
| 213 | ntmachk=ntmachk.filled() |
|---|
| 214 | |
|---|
| 215 | ### We iterate for each day in the chunk, on each grid point we find |
|---|
| 216 | ### the closest corresponding MCS grid point and the index of the |
|---|
| 217 | ### time in the chunk closest to the time in the closest MCS grid point. |
|---|
| 218 | ### (yea it's complicated) |
|---|
| 219 | |
|---|
| 220 | vartmpnight=np.zeros([nzvar,ndays,nz,ny,nx]) |
|---|
| 221 | vartmpday=np.zeros([nzvar,ndays,nz,ny,nx]) |
|---|
| 222 | vartmpnight=np.ma.masked_invalid(vartmpnight) |
|---|
| 223 | vartmpday=np.ma.masked_invalid(vartmpday) |
|---|
| 224 | vartmpnight.set_fill_value([np.NaN]) |
|---|
| 225 | vartmpday.set_fill_value([np.NaN]) |
|---|
| 226 | |
|---|
| 227 | nd=0 |
|---|
| 228 | print " .time indices MCS grid -> diagfi grid." |
|---|
| 229 | while nd < ndays: |
|---|
| 230 | |
|---|
| 231 | daystart=find_nearest(varchktime-varchktime[0],nd) |
|---|
| 232 | daystop=find_nearest(varchktime-varchktime[0],nd+1) |
|---|
| 233 | # varchktime_lon=np.zeros([daystop-daystart+1,len(lon)]) |
|---|
| 234 | ix=0 |
|---|
| 235 | for x in lon: |
|---|
| 236 | |
|---|
| 237 | varchktime_lon = 24.*(varchktime[daystart:daystop+1]-varchktime[daystart]) + x/15. |
|---|
| 238 | |
|---|
| 239 | iy=0 |
|---|
| 240 | for y in lat: |
|---|
| 241 | niy=find_nearest(latmcs,y) |
|---|
| 242 | nix=find_nearest(lonmcs,x) |
|---|
| 243 | nitdtmichk=find_nearest(varchktime_lon,dtmichk[niy,nix]) |
|---|
| 244 | nitdtmachk=find_nearest(varchktime_lon,dtmachk[niy,nix]) |
|---|
| 245 | nitntmichk=find_nearest(varchktime_lon,ntmichk[niy,nix]) |
|---|
| 246 | nitntmachk=find_nearest(varchktime_lon,ntmachk[niy,nix]) |
|---|
| 247 | for vv in np.arange(nzvar): |
|---|
| 248 | if ((nitdtmichk is np.NaN) or (nitdtmachk is np.NaN)): |
|---|
| 249 | vartmpday[vv,nd,:,iy,ix]=np.NaN |
|---|
| 250 | elif nitdtmichk > nitdtmachk: |
|---|
| 251 | 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. |
|---|
| 252 | else: |
|---|
| 253 | vartmpday[vv,nd,:,iy,ix]=np.ma.mean(varchk[vv,daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0) |
|---|
| 254 | if ((nitntmichk is np.NaN) or (nitntmachk is np.NaN)): |
|---|
| 255 | vartmpnight[vv,nd,:,iy,ix]=np.NaN |
|---|
| 256 | elif nitntmichk > nitntmachk: |
|---|
| 257 | 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. |
|---|
| 258 | else: |
|---|
| 259 | vartmpnight[vv,nd,:,iy,ix]=np.ma.mean(varchk[vv,daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0) |
|---|
| 260 | iy=iy+1 |
|---|
| 261 | ix=ix+1 |
|---|
| 262 | nd=nd+1 |
|---|
| 263 | |
|---|
| 264 | print " .creating bins." |
|---|
| 265 | |
|---|
| 266 | vartmpdaymasked=np.ma.masked_invalid(vartmpday) |
|---|
| 267 | vartmpnightmasked=np.ma.masked_invalid(vartmpnight) |
|---|
| 268 | vartmpdaymasked.set_fill_value([np.NaN]) |
|---|
| 269 | vartmpnightmasked.set_fill_value([np.NaN]) |
|---|
| 270 | for vv in np.arange(nzvar): |
|---|
| 271 | vardayout[vv,i,:,:,:]=np.ma.mean(vartmpdaymasked[vv,:,:,:,:],axis=0) |
|---|
| 272 | varnightout[vv,i,:,:,:]=np.ma.mean(vartmpnightmasked[vv,:,:,:,:],axis=0) |
|---|
| 273 | print " ."+varznames[vv]+".done" |
|---|
| 274 | i=i+1 |
|---|
| 275 | |
|---|
| 276 | print "--->> Preparing Data for ncdf. Missing value is NaN." |
|---|
| 277 | |
|---|
| 278 | vardayout=np.ma.masked_invalid(vardayout) |
|---|
| 279 | varnightout=np.ma.masked_invalid(varnightout) |
|---|
| 280 | vardayout.set_fill_value([np.NaN]) |
|---|
| 281 | varnightout.set_fill_value([np.NaN]) |
|---|
| 282 | |
|---|
| 283 | all=[aps,bps] |
|---|
| 284 | for vv in np.arange(nzvar): |
|---|
| 285 | if dimensions[vv] == 3: |
|---|
| 286 | all.append(vardayout[vv,:,0,:,:].filled()) |
|---|
| 287 | all.append(varnightout[vv,:,0,:,:].filled()) |
|---|
| 288 | elif dimensions[vv] == 4: |
|---|
| 289 | all.append(vardayout[vv,:,:,:,:].filled()) |
|---|
| 290 | all.append(varnightout[vv,:,:,:,:].filled()) |
|---|
| 291 | |
|---|
| 292 | if opt.recast: |
|---|
| 293 | all=all[2:] |
|---|
| 294 | fullnames=fullnames[2:] |
|---|
| 295 | |
|---|
| 296 | make_gcm_netcdf (zfilename="diagfi_MCS.nc", \ |
|---|
| 297 | zdescription="Temperatures from diagfi reworked to match MCS format", \ |
|---|
| 298 | zlon=lon, \ |
|---|
| 299 | zlat=lat, \ |
|---|
| 300 | zalt=alt, \ |
|---|
| 301 | ztime=timemcs, \ |
|---|
| 302 | zvariables=all, \ |
|---|
| 303 | znames=fullnames) |
|---|
| 304 | if opt.zrecast and opt.ditch: |
|---|
| 305 | print "removing interpolated file" |
|---|
| 306 | system("rm -f "+opt.file[0:len(opt.file)-3]+"_P.nc") |
|---|
