Changeset 2541 in lmdz_wrf for trunk/tools


Ignore:
Timestamp:
May 17, 2019, 5:21:17 AM (7 years ago)
Author:
lfita
Message:

Fixing `spacemean'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/nc_var_tools.py

    r2537 r2541  
    31423142
    31433143def spacemean(ncfile, varn):
    3144   """ Function to retrieve a space mean series from a multidimensional variable of a file
    3145   ncfile = netCDF file name
    3146   varn = variable name
    3147   """
    3148   import datetime as dt
    3149   import calendar as cal
    3150 
    3151   ofile = 'spacemean_' + varn + '.nc'
    3152   varfil=1.e20
    3153   statsn = ['min', 'max', 'mean', 'mean2', 'stdv', 'meanwgt', 'mean2wgt', 'stdvwgt', 'quant']
    3154   statslongn = ['minimum', 'maximum', 'mean', 'quadratic mean', 'standard deviation', 'weigthed mean',  'weigthed quadratic mean', 'weigthed standard deviation', 'quantiles']
    3155 
    3156   ncf = NetCDFFile(ncfile,'r')
    3157 
    3158   if not ncf.variables.has_key(varn):
    3159     print errormsg
    3160     print '    spacemean: File "' + ncfile + '" does not have variable "' + varn + '" !!!!'
    3161     print errormsg
    3162     ncf.close()
    3163     quit(-1)
    3164 
    3165   times = ncf.variables['time']
    3166   timevals = times[:]
    3167 
    3168   attvar = times.ncattrs()
    3169   if not gen.searchInlist(attvar, 'units'):
    3170     print errormsg
    3171     print '    spacemean: "time" does not have attribute: "units"'
    3172     ncf.close()
    3173     quit(-1)
    3174   else:
    3175     units = times.getncattr('units')
     3144    """ Function to retrieve a space mean series from a multidimensional variable of a file
     3145    ncfile = netCDF file name
     3146    varn = variable name
     3147    """
     3148    fname = 'spacemean'
     3149    import datetime as dt
     3150    import calendar as cal
     3151    tunitsavail = ['weeks', 'day', 'hours', 'minutes', 'seconds', 'miliseconds']
     3152
     3153    if varn == '-h':
     3154        print fname + '_____________________________________________________________'
     3155        print spacemean.__doc__
     3156        quit()
     3157
     3158    ofile = 'spacemean_' + varn + '.nc'
     3159    varfil = gen.fillValueF
     3160    statsn = ['min', 'max', 'mean', 'mean2', 'stdv', 'meanwgt', 'mean2wgt',          \
     3161      'stdvwgt', 'quant']
     3162    statslongn = ['minimum', 'maximum', 'mean', 'quadratic mean',                    \
     3163      'standard deviation', 'weigthed mean',  'weigthed quadratic mean',             \
     3164      'weigthed standard deviation', 'quantiles']
     3165
     3166    ncf = NetCDFFile(ncfile,'r')
     3167
     3168    if not ncf.variables.has_key(varn):
     3169        print errormsg
     3170        print '    ' + fname + ": File '" + ncfile + "' does not have variable '" +  \
     3171          varn + "' !!!!"
     3172        Varn = list(ncf.variables.keys())
     3173        Varn.sort()
     3174        print '    available ones:', Varn
     3175        ncf.close()
     3176        quit(-1)
     3177
     3178    times = ncf.variables['time']
     3179    timevals = times[:]
     3180
     3181    attvar = times.ncattrs()
     3182    if not gen.searchInlist(attvar, 'units'):
     3183        print errormsg
     3184        print '    spacemean: "time" does not have attribute: "units"'
     3185        ncf.close()
     3186        quit(-1)
     3187    else:
     3188        units = times.getncattr('units')
    31763189 
    3177   txtunits = units.split(' ')
    3178   tunits = txtunits[0]
    3179   Srefdate = txtunits[len(txtunits) - 1]
    3180 # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
    3181 ##
    3182   timeval = Srefdate.find(':')
    3183 
    3184   if not timeval == -1:
    3185     print '      spacemean: refdate with time!'
    3186     refdate = datetimeStr_datetime(txtunits[len(txtunits) - 2] + '_' + Srefdate)
    3187   else:
    3188     refdate = dateStr_date(Srefdate)
    3189 
    3190   dimt = len(timevals)
    3191   realdates = np.zeros((dimt, 6), dtype=int)
    3192   print realdates.shape
     3190    timeuinf = gen.CFtimeU_inf(units)
     3191
     3192    dimt = len(timevals)
     3193    realdates = np.zeros((dimt, 6), dtype=int)
    31933194
    31943195## Not in timedelta
    3195 #  if tunits == 'years':
    3196 #    for it in range(dimt):
    3197 #      realdate = refdate + dt.timedelta(years=float(times[it]))
    3198 #      realdates[it] = int(realdate.year)
    3199 #  elif tunits == 'months':
    3200 #    for it in range(dimt):
    3201 #      realdate = refdate + dt.timedelta(months=float(times[it]))
    3202 #      realdates[it] = int(realdate.year)
    3203   if tunits == 'weeks':
    3204     for it in range(dimt):
    3205       realdate = refdate + dt.timedelta(weeks=float(times[it]))
    3206       realdates[it,0] = int(realdate.year)
    3207       realdates[it,1] = int(realdate.month)
    3208       realdates[it,2] = int(realdate.day)
    3209       realdates[it,3] = int(realdate.hour)
    3210       realdates[it,4] = int(realdate.second)
    3211       realdates[it,5] = int(realdate.minute)
    3212   elif tunits == 'days':
    3213     for it in range(dimt):
    3214       realdate = refdate + dt.timedelta(days=float(times[it]))
    3215       realdates[it,0] = int(realdate.year)
    3216       realdates[it,1] = int(realdate.month)
    3217       realdates[it,2] = int(realdate.day)
    3218       realdates[it,3] = int(realdate.hour)
    3219       realdates[it,4] = int(realdate.second)
    3220       realdates[it,5] = int(realdate.minute)
    3221   elif tunits == 'hours':
    3222     for it in range(dimt):
    3223       realdate = refdate + dt.timedelta(hours=float(times[it]))
    3224       realdates[it,0] = int(realdate.year)
    3225       realdates[it,1] = int(realdate.month)
    3226       realdates[it,2] = int(realdate.day)
    3227       realdates[it,3] = int(realdate.hour)
    3228       realdates[it,4] = int(realdate.second)
    3229       realdates[it,5] = int(realdate.minute)
    3230   elif tunits == 'minutes':
    3231     for it in range(dimt):
    3232       realdate = refdate + dt.timedelta(minutes=float(times[it]))
    3233       realdates[it,0] = int(realdate.year)
    3234       realdates[it,1] = int(realdate.month)
    3235       realdates[it,2] = int(realdate.day)
    3236       realdates[it,3] = int(realdate.hour)
    3237       realdates[it,4] = int(realdate.second)
    3238       realdates[it,5] = int(realdate.minute)
    3239   elif tunits == 'seconds':
    3240     for it in range(dimt):
    3241       realdate = refdate + dt.timedelta(seconds=float(times[it]))
    3242       realdates[it,0] = int(realdate.year)
    3243       realdates[it,1] = int(realdate.month)
    3244       realdates[it,2] = int(realdate.day)
    3245       realdates[it,3] = int(realdate.hour)
    3246       realdates[it,4] = int(realdate.second)
    3247       realdates[it,5] = int(realdate.minute)
    3248   elif tunits == 'milliseconds':
    3249     for it in range(dimt):
    3250       realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
    3251       realdates[it,0] = int(realdate.year)
    3252       realdates[it,1] = int(realdate.month)
    3253       realdates[it,2] = int(realdate.day)
    3254       realdates[it,3] = int(realdate.hour)
    3255       realdates[it,4] = int(realdate.second)
    3256       realdates[it,5] = int(realdate.minute)
    3257   else:
    3258     print errormsg
    3259     print '    spacemean: time units "' + tunits + '" not ready!!!!'
    3260     ncf.close()
    3261     quit(-1)
     3196#    if tunits == 'years':
     3197#        for it in range(dimt):
     3198#            realdate = refdate + dt.timedelta(years=float(times[it]))
     3199#            realdates[it] = int(realdate.year)
     3200#    elif tunits == 'months':
     3201#        for it in range(dimt):
     3202#            realdate = refdate + dt.timedelta(months=float(times[it]))
     3203#            realdates[it] = int(realdate.year)
     3204    if timeuinf.Tunits == 'weeks':
     3205        for it in range(dimt):
     3206            realdate = timeuinf.refdateDT + dt.timedelta(weeks=float(times[it]))
     3207            realdates[it,0] = int(realdate.year)
     3208            realdates[it,1] = int(realdate.month)
     3209            realdates[it,2] = int(realdate.day)
     3210            realdates[it,3] = int(realdate.hour)
     3211            realdates[it,4] = int(realdate.second)
     3212            realdates[it,5] = int(realdate.minute)
     3213    elif timeuinf.Tunits == 'days':
     3214        for it in range(dimt):
     3215            realdate = timeuinf.refdateDT + dt.timedelta(days=float(times[it]))
     3216            realdates[it,0] = int(realdate.year)
     3217            realdates[it,1] = int(realdate.month)
     3218            realdates[it,2] = int(realdate.day)
     3219            realdates[it,3] = int(realdate.hour)
     3220            realdates[it,4] = int(realdate.second)
     3221            realdates[it,5] = int(realdate.minute)
     3222    elif timeuinf.Tunits == 'hours':
     3223        for it in range(dimt):
     3224            realdate = timeuinf.refdateDT + dt.timedelta(hours=float(times[it]))
     3225            realdates[it,0] = int(realdate.year)
     3226            realdates[it,1] = int(realdate.month)
     3227            realdates[it,2] = int(realdate.day)
     3228            realdates[it,3] = int(realdate.hour)
     3229            realdates[it,4] = int(realdate.second)
     3230            realdates[it,5] = int(realdate.minute)
     3231    elif timeuinf.Tunits == 'minutes':
     3232        for it in range(dimt):
     3233            realdate = timeuinf.refdateDT + dt.timedelta(minutes=float(times[it]))
     3234            realdates[it,0] = int(realdate.year)
     3235            realdates[it,1] = int(realdate.month)
     3236            realdates[it,2] = int(realdate.day)
     3237            realdates[it,3] = int(realdate.hour)
     3238            realdates[it,4] = int(realdate.second)
     3239            realdates[it,5] = int(realdate.minute)
     3240    elif timeuinf.Tunits == 'seconds':
     3241        for it in range(dimt):
     3242            realdate = timeuinf.refdateDT + dt.timedelta(seconds=float(times[it]))
     3243            realdates[it,0] = int(realdate.year)
     3244            realdates[it,1] = int(realdate.month)
     3245            realdates[it,2] = int(realdate.day)
     3246            realdates[it,3] = int(realdate.hour)
     3247            realdates[it,4] = int(realdate.second)
     3248            realdates[it,5] = int(realdate.minute)
     3249    elif timeuinf.Tunits == 'milliseconds':
     3250        for it in range(dimt):
     3251            realdate = timeuinf.refdateDT + dt.timedelta(milliseconds=float(times[it]))
     3252            realdates[it,0] = int(realdate.year)
     3253            realdates[it,1] = int(realdate.month)
     3254            realdates[it,2] = int(realdate.day)
     3255            realdates[it,3] = int(realdate.hour)
     3256            realdates[it,4] = int(realdate.second)
     3257            realdates[it,5] = int(realdate.minute)
     3258    else:
     3259        print errormsg
     3260        print '    ' + ": time units '" + tunits + "' not ready!!!!"
     3261        print '    available ones', tunitsavail
     3262        ncf.close()
     3263        quit(-1)
    32623264
    32633265# Variable values (assuming time as first dimension)
    32643266##
    3265   var = ncf.variables[varn]
    3266   varinf = variable_inf(var)
    3267   if gen.searchInlist(varinf.attributes, 'scale_factor'):
    3268     scalefact = var.getncattr('scale_factor')
    3269     print '      spacemean: data has a scale factor of ', scalefact
    3270     varinf.dtype=type(np.float(1.))
    3271   else:
    3272     scalefact = 1.
    3273 
    3274   if gen.searchInlist(varinf.attributes, 'add_offset'):
    3275     offset = var.getncattr('add_offset')
    3276     print '      spacemean: data has an offset of ', offset
    3277   else:
    3278     offset = 0.
     3267    var = ncf.variables[varn]
     3268    varinf = variable_inf(var)
     3269    if gen.searchInlist(varinf.attributes, 'scale_factor'):
     3270        scalefact = var.getncattr('scale_factor')
     3271        print '      ' + fname + ": data has a scale factor of ", scalefact
     3272        varinf.dtype=type(np.float(1.))
     3273    else:
     3274        scalefact = 1.
     3275
     3276    if gen.searchInlist(varinf.attributes, 'add_offset'):
     3277        offset = var.getncattr('add_offset')
     3278        print '      ' + fname + ": data has an offset of ", offset
     3279    else:
     3280        offset = 0.
    32793281
    32803282#  varvals = var[:]
    3281   vardims = var.shape
    3282   varshape = len(vardims)
    3283   vardimns = var.dimensions
     3283    vardims = var.shape
     3284    varshape = len(vardims)
     3285    vardimns = var.dimensions
    32843286##  print '      spacemean: Shape of data: ',varshape, ':' , vardims
    3285   dimt=vardims[0]
    3286   vardnames = []
     3287    dimt=vardims[0]
     3288    vardnames = []
    32873289
    32883290# Spatial average spatial weigthed
    32893291##
    3290   lonvar = ncf.variables['lon']
    3291   latvar = ncf.variables['lat']
    3292 
    3293   if not len(lonvar.shape) == 2:
    3294     lonv = lonvar[:]
    3295     latv = latvar[:]
    3296     dx=lonvar.shape[0]
    3297     dy=latvar.shape[0]
    3298     lonval = np.zeros((dy, dx), dtype=np.float64)   
    3299     latval = np.zeros((dy, dx), dtype=np.float64)   
    3300     for iy in range(dy):
    3301       lonval[iy,:] = lonv
    3302     for ix in range(dx):
    3303       latval[:,ix] = latv
    3304   else:
    3305     lonval = lonvar[:]
    3306     latval = latvar[:]
    3307 
    3308   weightsv = abs(np.cos(latval*np.pi/180.))
    3309 
    3310   if varinf.Ndims == 1:
    3311     print errormsg
    3312     print '    spacemean: You can not compute a space mean for a ', varinf.Ndims, 'D var!!!'
    3313     ncf.close()
    3314     quit(-1)
    3315   elif varinf.Ndims== 2:
    3316     print errormsg
    3317     print '    spacemean: You can not compute a space mean for a ', varinf.Ndims, 'D var!!!'
    3318     ncf.close()
    3319     quit(-1)
    3320   elif varinf.Ndims == 3:
    3321     varstats = np.ones((dimt, 8), dtype=np.float64)
    3322     varstats = varstats*varfil
    3323     varquant = np.ones((dimt, 21), dtype=np.float64)
    3324     varquant = varquant*varfil
    3325     vardnames.append(vardimns[0])
    3326     dy=vardims[1]
    3327     dx=vardims[2]
    3328     for it in range(dimt):
    3329       if not varinf.FillValue is None:
    3330         varvl = var[it,:,:]
    3331         varval = np.where(varvl == varinf.FillValue, None, varvl)
    3332       else:
    3333         varval = var[it,:,:]
    3334       percendone(it,dimt,5,'computed space statistics')
    3335       varst = gen.statsVal(varval)
    3336       varstwgt  = statsValWeigthed(varval, weightsv)
    3337       varstats[it,0] = varst.minv
    3338       varstats[it,1] = varst.maxv
    3339       varstats[it,2] = varst.meanv
    3340       varstats[it,3] = varst.mean2v
    3341       varstats[it,4] = varst.stdv
    3342       varquant[it,:] = varst.quantilesv
    3343       varstats[it,5] = varstwgt.meanv
    3344       varstats[it,6] = varstwgt.mean2v
    3345       varstats[it,7] = varstwgt.stdv
    3346   elif varshape == 4:
    3347     varstats = np.ones((dimt, vardims[1], 8), dtype=np.float64)
    3348     varstats = varstats*varfil
    3349     varquant = np.ones((dimt, vardims[1], 21), dtype=np.float64)
    3350     varquant = varquant*varfil
    3351     vardimnames = (str(vardimns[0]), str(vardimns[1]))
    3352     vardnames.append(vardimns[0])
    3353     vardnames.append(vardimns[1])
    3354     dy=vardims[2]
    3355     dx=vardims[3]
    3356     for it in range(dimt):
    3357       for ik in range(vardims[1]):
    3358         if not varinf.FillValue is None:
    3359           varvl = var[it,ik,:,:]
    3360           varval = np.where(varvl == varinf.FillValue, None, varvl)
    3361         else:
    3362           varval = var[it,ik,:,:]
    3363         percendone(it*vardims[1]+ik,dimt*vardims[1],5,'computed space statistics')
    3364         varst = gen.statsVal(varval)
    3365         varstwgt  = statsValWeigthed(varval, weightsv)
    3366         varstats[it,ik,0] = varst.minv
    3367         varstats[it,ik,1] = varst.maxv
    3368         varstats[it,ik,2] = varst.meanv
    3369         varstats[it,ik,3] = varst.mean2v
    3370         varstats[it,ik,4] = varst.stdv
    3371         varquant[it,ik,:] = varst.quantilesv
    3372         varstats[it,ik,5] = varstwgt.meanv
    3373         varstats[it,ik,6] = varstwgt.mean2v
    3374         varstats[it,ik,7] = varstwgt.stdv
    3375   elif varshape == 5:
    3376     varstats = np.ones((dimt, vardims[1], vardims[2], 8), dtype=np.float64)
    3377     varstats = varstats*varfil
    3378     varquant = np.ones((dimt,  vardims[1], vardims[2], 21), dtype=np.float64)
    3379     varquant = varquant*varfil
    3380     vardnames.append(vardimns[0])
    3381     vardnames.append(vardimns[1])
    3382     vardnames.append(vardimns[2])
    3383     dy=vardims[3]
    3384     dx=vardims[4]
    3385     for it in range(dimt):
    3386       for ik in range(vardims[1]):
    3387         for il in range(vardims[2]):
    3388           if not varinf.FillValue is None:
    3389             varvl = var[it,ik,il,:,:]
    3390             varval = np.where(varvl == varinf.FillValue, None, varvl)
    3391           else:
    3392             varval = var[it,ik,il,:,:]
    3393 
    3394           percendone(it*vardims[1]*vardims[2]+ik*vardims[1]+il,dimt*vardims[1]*vardims[2],5,'computed space statistics')
    3395           varst = gen.statsVal(varval)
    3396           varstwgt  = statsValWeigthed(varval, weightsv)
    3397           varstats[it,ik,il,0] = varst.minv
    3398           varstats[it,ik,il,1] = varst.maxv
    3399           varstats[it,ik,il,2] = varst.meanv
    3400           varstats[it,ik,il,3] = varst.mean2v
    3401           varstats[it,ik,il,4] = varst.stdv
    3402           varquant[it,ik,il,:] = varst.quantilesv
    3403           varstats[it,ik,il,5] = varstwgt.meanv
    3404           varstats[it,ik,il,6] = varstwgt.mean2v
    3405           varstats[it,ik,il,7] = varstwgt.stdv
    3406   else:
    3407      print errormsg
    3408      print '    spacemean: ', varshape, ' shape of matrix not prepared !'
    3409      ncf.close()
    3410      quit(-1)
    3411 
    3412   vardimnames = tuple(vardnames)
     3292    olonvar = ncf.variables['lon']
     3293    olatvar = ncf.variables['lat']
     3294
     3295    lonval, latval = gen.lonlat2D(olonvar[:],olatvar[:])
     3296
     3297    weightsv = abs(np.cos(latval*np.pi/180.))
     3298
     3299    if varinf.Ndims == 1:
     3300        print errormsg
     3301        print '    ' + fname + ": You can not compute a space mean for a ",          \
     3302          varinf.Ndims, 'D var!!!'
     3303        ncf.close()
     3304        quit(-1)
     3305    elif varinf.Ndims== 2:
     3306        print errormsg
     3307        print '    ' + fname + ': You can not compute a space mean for a ',          \
     3308          varinf.Ndims, 'D var!!!'
     3309        ncf.close()
     3310        quit(-1)
     3311    elif varinf.Ndims == 3:
     3312        varstats = np.ones((dimt, 8), dtype=np.float64)
     3313        varstats = varstats*varfil
     3314        varquant = np.ones((dimt, 21), dtype=np.float64)
     3315        varquant = varquant*varfil
     3316        vardnames.append(vardimns[0])
     3317        dy=vardims[1]
     3318        dx=vardims[2]
     3319        for it in range(dimt):
     3320            if not varinf.FillValue is None:
     3321                varvl = var[it,:,:]
     3322                varval0 = np.where(varvl == varinf.FillValue, None, varvl)
     3323            else:
     3324                varval0 = var[it,:,:]
     3325            varval = varval0.flatten()
     3326            gen.percendone(it,dimt,5,'computed space statistics')
     3327            varst = gen.statsVal(varval)
     3328            varstwgt  = gen.statsValWeigthed(varval, weightsv.flatten())
     3329            varstats[it,0] = varst.minv
     3330            varstats[it,1] = varst.maxv
     3331            varstats[it,2] = varst.meanv
     3332            varstats[it,3] = varst.mean2v
     3333            varstats[it,4] = varst.stdv
     3334            varquant[it,:] = varst.quantilesv
     3335            varstats[it,5] = varstwgt.meanv
     3336            varstats[it,6] = varstwgt.mean2v
     3337            varstats[it,7] = varstwgt.stdv
     3338    elif varshape == 4:
     3339        varstats = np.ones((dimt, vardims[1], 8), dtype=np.float64)
     3340        varstats = varstats*varfil
     3341        varquant = np.ones((dimt, vardims[1], 21), dtype=np.float64)
     3342        varquant = varquant*varfil
     3343        vardimnames = (str(vardimns[0]), str(vardimns[1]))
     3344        vardnames.append(vardimns[0])
     3345        vardnames.append(vardimns[1])
     3346        dy=vardims[2]
     3347        dx=vardims[3]
     3348        for it in range(dimt):
     3349            for ik in range(vardims[1]):
     3350                if not varinf.FillValue is None:
     3351                    varvl = var[it,ik,:,:]
     3352                    varval0 = np.where(varvl == varinf.FillValue, None, varvl)
     3353                else:
     3354                    varval0 = var[it,ik,:,:]
     3355                varval = varval0.flatten()
     3356                gen.percendone(it*vardims[1]+ik,dimt*vardims[1],5,'computed space statistics')
     3357                varst = gen.statsVal(varval)
     3358                varstwgt  = gen.statsValWeigthed(varval, weightsv)
     3359                varstats[it,ik,0] = varst.minv
     3360                varstats[it,ik,1] = varst.maxv
     3361                varstats[it,ik,2] = varst.meanv
     3362                varstats[it,ik,3] = varst.mean2v
     3363                varstats[it,ik,4] = varst.stdv
     3364                varquant[it,ik,:] = varst.quantilesv
     3365                varstats[it,ik,5] = varstwgt.meanv
     3366                varstats[it,ik,6] = varstwgt.mean2v
     3367                varstats[it,ik,7] = varstwgt.stdv
     3368    elif varshape == 5:
     3369        varstats = np.ones((dimt, vardims[1], vardims[2], 8), dtype=np.float64)
     3370        varstats = varstats*varfil
     3371        varquant = np.ones((dimt,  vardims[1], vardims[2], 21), dtype=np.float64)
     3372        varquant = varquant*varfil
     3373        vardnames.append(vardimns[0])
     3374        vardnames.append(vardimns[1])
     3375        vardnames.append(vardimns[2])
     3376        dy=vardims[3]
     3377        dx=vardims[4]
     3378        for it in range(dimt):
     3379            for ik in range(vardims[1]):
     3380                for il in range(vardims[2]):
     3381                    if not varinf.FillValue is None:
     3382                        varvl = var[it,ik,il,:,:]
     3383                        varval = np.where(varvl == varinf.FillValue, None, varvl)
     3384                    else:
     3385                        varval = var[it,ik,il,:,:]
     3386
     3387                    gen.percendone(it*vardims[1]*vardims[2]+ik*vardims[1]+il,dimt*\
     3388                      vardims[1]*vardims[2],5,'computed space statistics')
     3389                    varst = gen.statsVal(varval)
     3390                    varstwgt = gen.statsValWeigthed(varval, weightsv)
     3391                    varstats[it,ik,il,0] = varst.minv
     3392                    varstats[it,ik,il,1] = varst.maxv
     3393                    varstats[it,ik,il,2] = varst.meanv
     3394                    varstats[it,ik,il,3] = varst.mean2v
     3395                    varstats[it,ik,il,4] = varst.stdv
     3396                    varquant[it,ik,il,:] = varst.quantilesv
     3397                    varstats[it,ik,il,5] = varstwgt.meanv
     3398                    varstats[it,ik,il,6] = varstwgt.mean2v
     3399                    varstats[it,ik,il,7] = varstwgt.stdv
     3400    else:
     3401        print errormsg
     3402        print '    ' + fname + ": shape ", varshape, " of matrix not prepared !!"
     3403        ncf.close()
     3404        quit(-1)
     3405
     3406    vardimnames = tuple(vardnames)
    34133407#  print '  shape of desired values: ', vardesiredvalues.shape
    34143408# Creation of file
    34153409##
    34163410
    3417   ncfo = NetCDFFile( ofile, 'w')
    3418 
    3419   vartype = var.dtype
    3420   varattr = var.ncattrs()
     3411    ncfo = NetCDFFile( ofile, 'w')
     3412
     3413    vartype = var.dtype
     3414    varattr = var.ncattrs()
    34213415
    34223416# Checking dimensions
    34233417##
    3424   newdims = ncfo.dimensions
    3425   Nvardnames = len(vardimnames)
    3426   for idim in range(Nvardnames):
    3427       rdim = vardimnames[idim]
    3428       if not gen.searchInlist(newdims, rdim):
    3429           if not rdim == 'time':
    3430               print '      spacemean: Adding dimension ' + rdim
    3431               ncfo.sync()
    3432               ncf.close()
    3433               ncfo.close()
    3434               fdimadd(ncfile + ',' + rdim, ofile)
    3435               ncf = NetCDFFile(ncfile,'r')
    3436               ncfo = NetCDFFile(ofile,'a')
    3437           else:
    3438               ncfo.createDimension('time', None)
     3418    newdims = ncfo.dimensions
     3419    Nvardnames = len(vardimnames)
     3420    for idim in range(Nvardnames):
     3421        rdim = vardimnames[idim]
     3422        if not gen.searchInlist(newdims, rdim):
     3423            if not rdim == 'time':
     3424                print '      spacemean: Adding dimension ' + rdim
     3425                ncfo.sync()
     3426                ncf.close()
     3427                ncfo.close()
     3428                fdimadd(ncfile + ',' + rdim, ofile)
     3429                ncf = NetCDFFile(ncfile,'r')
     3430                ncfo = NetCDFFile(ofile,'a')
     3431            else:
     3432                ncfo.createDimension('time', None)
    34393433
    34403434# Checking fill value
    34413435##
    3442   if gen.searchInlist(varattr, '_FillValue'):
    3443       varfil = var._FillValue
    3444   else:
    3445       varfil = False
    3446 
    3447   Nstats = len(statsn)
    3448   for ist in range(Nstats):
    3449     if statsn[ist] == 'quant':
    3450       print ist, statsn[ist]##, ': ', varquant[int(dimt/2),10]
    3451       newdim = ncfo.createDimension('quant', 21)
    3452       newvardnames = list(vardnames)
    3453       newvardnames.append('quant')
    3454       newvar = ncfo.createVariable(varn + statsn[ist], varinf.dtype, tuple(newvardnames), fill_value=varfil)
    3455       newvar[:] = varquant
     3436    if gen.searchInlist(varattr, '_FillValue'):
     3437        varfil = var._FillValue
    34563438    else:
    3457       print ist, statsn[ist]##, ': ', varstats[int(dimt/2),ist]
    3458       newvar = ncfo.createVariable(varn + statsn[ist], varinf.dtype, vardimnames, fill_value=varfil)
    3459       if varshape == 3:
    3460         newvar[:] = varstats[:,ist]*1.
    3461       elif varshape == 4:
    3462         newvar[:] = varstats[:,:,ist]*1.
    3463       elif varshape == 5:
    3464         newvar[:] = varstats[:,:,:,ist]*1.
    3465 
    3466     newvar = ncfo.variables[varn + statsn[ist]]
    3467     for attr in varattr:
    3468          newvarattrs = newvar.ncattrs()
    3469          attrv = var.getncattr(attr)
    3470          if not gen.searchInlist(newvarattrs, attr):
    3471              if attr != 'scale_factor' and attr != 'add_offset' and attr != 'valid_range' \
    3472                 and attr != 'unpacked_valid_range' and attr != 'actual_range' :
     3439        varfil = False
     3440
     3441    Nstats = len(statsn)
     3442    for ist in range(Nstats):
     3443        if statsn[ist] == 'quant':
     3444            print ist, statsn[ist]##, ': ', varquant[int(dimt/2),10]
     3445            newdim = ncfo.createDimension('quant', 21)
     3446            newvardnames = list(vardnames)
     3447            newvardnames.append('quant')
     3448            newvar = ncfo.createVariable(varn + statsn[ist], varinf.dtype, tuple(newvardnames), fill_value=varfil)
     3449            newvar[:] = varquant
     3450        else:
     3451            print ist, statsn[ist]##, ': ', varstats[int(dimt/2),ist]
     3452            newvar = ncfo.createVariable(varn + statsn[ist], varinf.dtype, vardimnames, fill_value=varfil)
     3453            if varshape == 3:
     3454                newvar[:] = varstats[:,ist]*1.
     3455            elif varshape == 4:
     3456                newvar[:] = varstats[:,:,ist]*1.
     3457            elif varshape == 5:
     3458                newvar[:] = varstats[:,:,:,ist]*1.
     3459
     3460        newvar = ncfo.variables[varn + statsn[ist]]
     3461        for attr in varattr:
     3462             newvarattrs = newvar.ncattrs()
     3463             attrv = var.getncattr(attr)
     3464             if not gen.searchInlist(newvarattrs, attr):
     3465                 if attr != 'scale_factor' and attr != 'add_offset' and attr != 'valid_range' \
     3466                    and attr != 'unpacked_valid_range' and attr != 'actual_range' :
    34733467#                 if not statsn[ist] == 'stdv' and not statsn[ist] == 'stdvwgt':
    3474                      print 'attr:', attr
    3475                      newvar.setncattr(attr, attrv)
     3468                         print 'attr:', attr
     3469                         newvar.setncattr(attr, attrv)
    34763470#             else:
    34773471#                 newvar.setncattr(attr, attrv)
    34783472   
    3479     newvar.setncattr('cell_methods', 'space ' + statslongn[ist] + ' all domain ' + str(dy) + 'x' + str(dx))
     3473        newvar.setncattr('cell_methods', 'space ' + statslongn[ist] + ' all domain '+\
     3474          str(dy) + 'x' + str(dx))
    34803475
    34813476##  print '  Adding time variable'
    3482   vartdims = times.dimensions
    3483   vartype = times.dtype
    3484   varattr = times.ncattrs()
    3485 
    3486   newvar = ncfo.createVariable('time', vartype, vartdims, fill_value=varfil)
    3487   newvar = ncfo.variables['time']
    3488   newvar[:] = times
    3489 
    3490   ncf.close()
    3491   ncfo.sync()
    3492   ncfo.close()
    3493   fattradd('time', ncfile + ':time', ofile)
    3494   fvaradd(ncfile + ',lon', ofile)
    3495   fvaradd(ncfile + ',lat', ofile)
    3496 
    3497   ncfo = NetCDFFile(ofile,'a')
    3498   newvar = ncfo.variables['time']
    3499   newvarattrs = newvar.ncattrs()
    3500   if gen.searchInlist(newvarattrs, 'bounds'):
    3501       if newvar.getncattr('bounds') == 'time_bnds':
    3502           ncf = NetCDFFile(ncfile,'r')
    3503           tbnds = ncf.variables['time_bnds']
    3504           vardims = tbnds.dimensions
    3505           vartype = tbnds.dtype
    3506           varattr = tbnds.ncattrs()
    3507           ncfo.createDimension('bnds', 2)
    3508           newvar = ncfo.createVariable('time_bnds', vartype, vardims, fill_value=varfil)
    3509           newvar[:] = tbnds[:]
    3510 
    3511           ncf.close()
    3512           ncfo.sync()
    3513           ncfo.close()
    3514           fattradd('time_bnds', ncfile + ':time_bnds', ofile)
    3515       else:
    3516           ncfo.close()
    3517   else:
    3518       ncfo.close()
    3519 
    3520   fgaddattr(ncfile, ofile)
    3521 
    3522   print '    spacemean: File "' + ofile + '" as space mean of "' + varn + '" has been created'
    3523 
    3524   return
     3477    vartdims = times.dimensions
     3478    vartype = times.dtype
     3479    varattr = times.ncattrs()
     3480
     3481    newvar = ncfo.createVariable('time', vartype, vartdims, fill_value=varfil)
     3482    #newvar = ncfo.variables['time']
     3483    newvar[:] = times[:]
     3484
     3485    ncf.close()
     3486    ncfo.sync()
     3487    ncfo.close()
     3488    fattradd('time', ncfile + ':time', ofile)
     3489    fvaradd(ncfile + ',lon', ofile)
     3490    fvaradd(ncfile + ',lat', ofile)
     3491
     3492    ncfo = NetCDFFile(ofile,'a')
     3493    newvar = ncfo.variables['time']
     3494    newvarattrs = newvar.ncattrs()
     3495    if gen.searchInlist(newvarattrs, 'bounds'):
     3496        if newvar.getncattr('bounds') == 'time_bnds':
     3497            ncf = NetCDFFile(ncfile,'r')
     3498            tbnds = ncf.variables['time_bnds']
     3499            vardims = tbnds.dimensions
     3500            vartype = tbnds.dtype
     3501            varattr = tbnds.ncattrs()
     3502            ncfo.createDimension('bnds', 2)
     3503            newvar = ncfo.createVariable('time_bnds', vartype, vardims, fill_value=varfil)
     3504            newvar[:] = tbnds[:]
     3505
     3506            ncf.close()
     3507            ncfo.sync()
     3508            ncfo.close()
     3509            fattradd('time_bnds', ncfile + ':time_bnds', ofile)
     3510        else:
     3511            ncfo.close()
     3512    else:
     3513        ncfo.close()
     3514
     3515    fgaddattr(ncfile, ofile)
     3516    ncfo = NetCDFFile(ofile, 'a')
     3517    add_global_PyNCplot(ncfo, 'nc_var_tools', fname, '1.1')
     3518
     3519    print '    ' + fname + ": File '" + ofile + "' as space mean of '" + varn +      \
     3520      "' has been created"
     3521
     3522    return
    35253523
    35263524def timemean(values, ncfile, varn):
    … …  
    37843782                        varval = var[:,iz,iy,ix]
    37853783     
    3786                     gen.percendone(iz*dy*dx+iy*dx+ix,dz*dx*dy,5,'computed time statistics')
     3784                    gen.percendone(iz*dy*dx+iy*dx+ix,dz*dx*dy, 5,                    \
     3785                      'computed time statistics')
    37873786                    varst = gen.statsVal(varval)
    37883787                    var2st = gen.stats2Val(timesv, varval, powerv)
    … …  
    38293828                            varval = var[:,iN,iz,iy,ix]
    38303829     
    3831                         gen.percendone(iN*dy*dx*dz+iz*dy*dx+iy*dx+ix,dn*dz*dx*dy,5,'computed time statistics')
     3830                        gen.percendone(iN*dy*dx*dz+iz*dy*dx+iy*dx+ix,dn*dz*dx*dy,5,  \
     3831                          'computed time statistics')
    38323832                        varst = gen.statsVal(varval)
    38333833                        var2st = gen.stats2Val(timesv, varval, powerv)
    … …  
    1266512665                for i in range(Ngrid2,dh1-Ngrid2+1):
    1266612666                    for j in range(Ngrid2,dh2-Ngrid2+1):
    12667                         percendone(i*(dh1-Ngrid)*1.+j,(dh1-Ngrid)*(dh2-Ngrid)*1., 1.,  \
    12668                           vn + ' filtered: ' + str(i) + ', ' + str(j))
     12667                        gen.percendone(i*(dh1-Ngrid)*1.+j,(dh1-Ngrid)*(dh2-Ngrid)*1., \
     12668                          1., vn + ' filtered: ' + str(i) + ', ' + str(j))
    1266912669                        varslice = []
    1267012670                        varslicesum = []
Note: See TracChangeset for help on using the changeset viewer.