| 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 |
|---|
| 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 | |
|---|
| 33 | ############################# |
|---|
| 34 | ### Get options and variables |
|---|
| 35 | (opt,args) = parser.parse_args() |
|---|
| 36 | |
|---|
| 37 | ############################# |
|---|
| 38 | ### Load and check data |
|---|
| 39 | |
|---|
| 40 | # Files |
|---|
| 41 | |
|---|
| 42 | nc=Dataset(opt.file) |
|---|
| 43 | ncmcs=Dataset(opt.mcsfile) |
|---|
| 44 | |
|---|
| 45 | # Dimensions |
|---|
| 46 | |
|---|
| 47 | lon=nc.variables["longitude"][:] |
|---|
| 48 | lat=nc.variables["latitude"][:] |
|---|
| 49 | alt=nc.variables["altitude"][:] |
|---|
| 50 | time=nc.variables["Time"][:] # in fraction of sols |
|---|
| 51 | controle=nc.variables["controle"][:] |
|---|
| 52 | day_ini=controle[3]%669 |
|---|
| 53 | time[:]=time[:]+day_ini |
|---|
| 54 | nx=len(lon) |
|---|
| 55 | ny=len(lat) |
|---|
| 56 | nz=len(alt) |
|---|
| 57 | nt=len(time) |
|---|
| 58 | lstime=sol2ls(time) |
|---|
| 59 | |
|---|
| 60 | # MCS |
|---|
| 61 | |
|---|
| 62 | dtimemintmp=ncmcs.variables["dtimemin"][:,:,:] |
|---|
| 63 | dtimemaxtmp=ncmcs.variables["dtimemax"][:,:,:] |
|---|
| 64 | ntimemintmp=ncmcs.variables["ntimemin"][:,:,:] |
|---|
| 65 | ntimemaxtmp=ncmcs.variables["ntimemax"][:,:,:] |
|---|
| 66 | lonmcs=ncmcs.variables["longitude"][:] |
|---|
| 67 | latmcs=ncmcs.variables["latitude"][:] |
|---|
| 68 | timemcs=ncmcs.variables["time"][:]%360 # IN LS |
|---|
| 69 | |
|---|
| 70 | dtimemin=np.ma.masked_where(dtimemintmp < 0.,dtimemintmp) |
|---|
| 71 | dtimemin.set_fill_value([np.NaN]) |
|---|
| 72 | dtimemax=np.ma.masked_where(dtimemaxtmp < 0.,dtimemaxtmp) |
|---|
| 73 | dtimemax.set_fill_value([np.NaN]) |
|---|
| 74 | ntimemin=np.ma.masked_where(ntimemintmp < 0.,ntimemintmp) |
|---|
| 75 | ntimemin.set_fill_value([np.NaN]) |
|---|
| 76 | ntimemax=np.ma.masked_where(ntimemaxtmp < 0.,ntimemaxtmp) |
|---|
| 77 | ntimemax.set_fill_value([np.NaN]) |
|---|
| 78 | |
|---|
| 79 | # Variables to treat |
|---|
| 80 | |
|---|
| 81 | varz=[] |
|---|
| 82 | varznames=separatenames(opt.var[0]) |
|---|
| 83 | n=0 |
|---|
| 84 | for zn in varznames: |
|---|
| 85 | varz.append(getfield(nc,zn)) |
|---|
| 86 | print "Found: "+zn+" with dimensions: " |
|---|
| 87 | print np.array(varz[n]).shape |
|---|
| 88 | n=n+1 |
|---|
| 89 | |
|---|
| 90 | nzvar=len(varz) |
|---|
| 91 | dimensions={} |
|---|
| 92 | vv=0 |
|---|
| 93 | for var in varz: |
|---|
| 94 | a=len(np.array(var).shape) |
|---|
| 95 | if a == 3: dimensions[vv]=a |
|---|
| 96 | elif a == 4: dimensions[vv]=a |
|---|
| 97 | else: |
|---|
| 98 | print "Warning, only 3d and 4d variables are supported for time-recasting" |
|---|
| 99 | exit() |
|---|
| 100 | vv=vv+1 |
|---|
| 101 | |
|---|
| 102 | # Variables to save but not treated (only along z, or phisinit-like files) |
|---|
| 103 | |
|---|
| 104 | aps=nc.variables["aps"][:] |
|---|
| 105 | bps=nc.variables["bps"][:] |
|---|
| 106 | fullnames=["aps","bps"] |
|---|
| 107 | for name in varznames: |
|---|
| 108 | fullnames.append("d"+name) |
|---|
| 109 | fullnames.append("n"+name) |
|---|
| 110 | print "Will output: " |
|---|
| 111 | print fullnames |
|---|
| 112 | ############################# |
|---|
| 113 | ### Building |
|---|
| 114 | ############################# |
|---|
| 115 | |
|---|
| 116 | ### We loop over chunks of gcm data corresponding to MCS time dimension |
|---|
| 117 | ### Bin sizes for mcs data is 5 degrees ls centered on value |
|---|
| 118 | varday=np.zeros([len(timemcs),nz,ny,nx]) |
|---|
| 119 | varnight=np.zeros([len(timemcs),nz,ny,nx]) |
|---|
| 120 | vardayout=np.zeros([nzvar,len(timemcs),nz,ny,nx]) |
|---|
| 121 | varnightout=np.zeros([nzvar,len(timemcs),nz,ny,nx]) |
|---|
| 122 | vardayout=np.ma.masked_invalid(vardayout) |
|---|
| 123 | varnightout=np.ma.masked_invalid(varnightout) |
|---|
| 124 | i=0 |
|---|
| 125 | for ls in timemcs: |
|---|
| 126 | lsstart=ls-2.5 |
|---|
| 127 | lsstop=ls+2.5 |
|---|
| 128 | istart=find_nearest(lstime,lsstart,strict=True) |
|---|
| 129 | istop=find_nearest(lstime,lsstop,strict=True) |
|---|
| 130 | varchk=0 |
|---|
| 131 | if ((istart is np.NaN) or (istop is np.NaN)): |
|---|
| 132 | vardayout[:,i,:,:,:]=np.NaN |
|---|
| 133 | varnightout[:,i,:,:,:]=np.NaN |
|---|
| 134 | print "Time interval skipped. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' |
|---|
| 135 | i=i+1 |
|---|
| 136 | continue |
|---|
| 137 | print "--->> Processing Data. Ls MCS: (",lsstart,';',lsstop,')',"// Ls diagfi: (",lstime.min(),';',lstime.max(),')' |
|---|
| 138 | # 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. |
|---|
| 139 | # 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 |
|---|
| 140 | print " .initialisation." |
|---|
| 141 | |
|---|
| 142 | |
|---|
| 143 | varchk=np.zeros([nzvar,istop-istart+1,nz,ny,nx]) |
|---|
| 144 | vv=0 |
|---|
| 145 | for variable in varz: |
|---|
| 146 | if dimensions[vv] is 3: |
|---|
| 147 | varchk[vv,:,0,:,:]=variable[istart:istop+1,:,:] |
|---|
| 148 | else: |
|---|
| 149 | varchk[vv,:,:,:,:]=variable[istart:istop+1,:,:,:] |
|---|
| 150 | vv=vv+1 |
|---|
| 151 | varchktime=time[istart:istop+1] |
|---|
| 152 | ndays=np.floor(varchktime[len(varchktime)-1])-np.floor(varchktime[0]) |
|---|
| 153 | dtmichk=dtimemin[i,:,:] |
|---|
| 154 | dtmachk=dtimemax[i,:,:] |
|---|
| 155 | ntmichk=ntimemin[i,:,:] |
|---|
| 156 | ntmachk=ntimemax[i,:,:] |
|---|
| 157 | dtmichk.set_fill_value(np.NaN) |
|---|
| 158 | dtmachk.set_fill_value(np.NaN) |
|---|
| 159 | ntmichk.set_fill_value(np.NaN) |
|---|
| 160 | ntmachk.set_fill_value(np.NaN) |
|---|
| 161 | dtmichk=dtmichk.filled() |
|---|
| 162 | dtmachk=dtmachk.filled() |
|---|
| 163 | ntmichk=ntmichk.filled() |
|---|
| 164 | ntmachk=ntmachk.filled() |
|---|
| 165 | |
|---|
| 166 | ### We iterate for each day in the chunk, on each grid point we find |
|---|
| 167 | ### the closest corresponding MCS grid point and the index of the |
|---|
| 168 | ### time in the chunk closest to the time in the closest MCS grid point. |
|---|
| 169 | ### (yea it's complicated) |
|---|
| 170 | |
|---|
| 171 | vartmpnight=np.zeros([nzvar,ndays,nz,ny,nx]) |
|---|
| 172 | vartmpday=np.zeros([nzvar,ndays,nz,ny,nx]) |
|---|
| 173 | nd=0 |
|---|
| 174 | print " .time indices MCS grid -> diagfi grid." |
|---|
| 175 | while nd < ndays: |
|---|
| 176 | |
|---|
| 177 | daystart=find_nearest(varchktime-varchktime[0],nd) |
|---|
| 178 | daystop=find_nearest(varchktime-varchktime[0],nd+1) |
|---|
| 179 | # varchktime_lon=np.zeros([daystop-daystart+1,len(lon)]) |
|---|
| 180 | ix=0 |
|---|
| 181 | for x in lon: |
|---|
| 182 | |
|---|
| 183 | varchktime_lon = 24.*(varchktime[daystart:daystop+1]-varchktime[daystart]) + x/15. |
|---|
| 184 | |
|---|
| 185 | iy=0 |
|---|
| 186 | for y in lat: |
|---|
| 187 | niy=find_nearest(latmcs,y) |
|---|
| 188 | nix=find_nearest(lonmcs,x) |
|---|
| 189 | nitdtmichk=find_nearest(varchktime_lon,dtmichk[niy,nix]) |
|---|
| 190 | nitdtmachk=find_nearest(varchktime_lon,dtmachk[niy,nix]) |
|---|
| 191 | nitntmichk=find_nearest(varchktime_lon,ntmichk[niy,nix]) |
|---|
| 192 | nitntmachk=find_nearest(varchktime_lon,ntmachk[niy,nix]) |
|---|
| 193 | for vv in np.arange(nzvar): |
|---|
| 194 | if ((nitdtmichk is np.NaN) or (nitdtmachk is np.NaN)): |
|---|
| 195 | vartmpday[vv,nd,:,iy,ix]=np.NaN |
|---|
| 196 | elif nitdtmichk > nitdtmachk: |
|---|
| 197 | vartmpday[vv,nd,:,iy,ix]=(np.mean(varchk[vv,daystart+nitdtmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[vv,daystart:daystart+nitdtmachk+1,:,iy,ix],axis=0))/2. |
|---|
| 198 | else: |
|---|
| 199 | vartmpday[vv,nd,:,iy,ix]=np.mean(varchk[vv,daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0) |
|---|
| 200 | if ((nitntmichk is np.NaN) or (nitntmachk is np.NaN)): |
|---|
| 201 | vartmpnight[vv,nd,:,iy,ix]=np.NaN |
|---|
| 202 | elif nitntmichk > nitntmachk: |
|---|
| 203 | vartmpnight[vv,nd,:,iy,ix]=(np.mean(varchk[vv,daystart+nitntmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[vv,daystart:daystart+nitntmachk+1,:,iy,ix],axis=0))/2. |
|---|
| 204 | else: |
|---|
| 205 | vartmpnight[vv,nd,:,iy,ix]=np.mean(varchk[vv,daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0) |
|---|
| 206 | iy=iy+1 |
|---|
| 207 | ix=ix+1 |
|---|
| 208 | nd=nd+1 |
|---|
| 209 | |
|---|
| 210 | print " .creating bins." |
|---|
| 211 | |
|---|
| 212 | vartmpdaymasked=np.ma.masked_invalid(vartmpday) |
|---|
| 213 | vartmpnightmasked=np.ma.masked_invalid(vartmpnight) |
|---|
| 214 | for vv in np.arange(nzvar): |
|---|
| 215 | vardayout[vv,i,:,:,:]=np.ma.mean(vartmpdaymasked[vv,:,:,:,:],axis=0) |
|---|
| 216 | varnightout[vv,i,:,:,:]=np.ma.mean(vartmpnightmasked[vv,:,:,:,:],axis=0) |
|---|
| 217 | print " ."+varznames[vv]+".done" |
|---|
| 218 | i=i+1 |
|---|
| 219 | |
|---|
| 220 | print "--->> Preparing Data for ncdf. Missing value is NaN." |
|---|
| 221 | |
|---|
| 222 | vardayout=np.ma.masked_invalid(vardayout) |
|---|
| 223 | varnightout=np.ma.masked_invalid(varnightout) |
|---|
| 224 | vardayout.set_fill_value([np.NaN]) |
|---|
| 225 | varnightout.set_fill_value([np.NaN]) |
|---|
| 226 | |
|---|
| 227 | all=[aps,bps] |
|---|
| 228 | for vv in np.arange(nzvar): |
|---|
| 229 | if dimensions[vv] == 3: |
|---|
| 230 | all.append(vardayout[vv,:,0,:,:].filled()) |
|---|
| 231 | all.append(varnightout[vv,:,0,:,:].filled()) |
|---|
| 232 | elif dimensions[vv] == 4: |
|---|
| 233 | all.append(vardayout[vv,:,:,:,:].filled()) |
|---|
| 234 | all.append(varnightout[vv,:,:,:,:].filled()) |
|---|
| 235 | |
|---|
| 236 | make_gcm_netcdf (zfilename="diagfi_MCS.nc", \ |
|---|
| 237 | zdescription="Temperatures from diagfi reworked to match MCS format", \ |
|---|
| 238 | zlon=lon, \ |
|---|
| 239 | zlat=lat, \ |
|---|
| 240 | zalt=alt, \ |
|---|
| 241 | ztime=timemcs, \ |
|---|
| 242 | zvariables=all, \ |
|---|
| 243 | znames=fullnames) |
|---|
