Changeset 2228 in lmdz_wrf


Ignore:
Timestamp:
Nov 14, 2018, 9:08:39 PM (6 years ago)
Author:
lfita
Message:

Fixing `cyclevar_within'

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/generic_tools.py

    r2227 r2228  
    1402314023    >>> cyclevar_within(np.arange(12)+1, 12, 3, 5)
    1402414024    False
     14025    >>> cyclevar_within(np.arange(12)+1, 12, 2, 2)
     14026    True
    1402514027    """
    1402614028    fname = 'cyclevar_within'
    … …  
    1403614038    if eind < bind:
    1403714039        eind = eind + dcyc
     14040        vind = vind + dcyc
    1403814041
    1403914042    if vind < bind:
    1404014043        vind = vind + dcyc
    1404114044
    14042     if vind >= bind and vind < eind: within = True
     14045    #print bind, '<=',  vind, '<', eind
     14046
     14047    if vind >= bind and vind <= eind: within = True
    1404314048    else: within = False
    1404414049
    1404514050    return within
    14046 
    14047 def time_slices(tv, tu, cal, per, amount, temp_desc):
    14048     """ Function to return temporal slices of a series of times for a given amount
    14049         of periods
    14050       tv: time vallues (CF format)
    14051       tu: CF time-units
    14052       cal: calendar
    14053       per: period
    14054         'year': full year [01/01 00:00:00 - 12/31 23:59:59]
    14055         'month': full months [01 00:00:00 - [31,30,28/29] 23:59:59]
    14056         'day': full days [00:00:00 - 23:59:59]
    14057         'hour': full hours [00:00 - 59:59]
    14058         'minute': full minutes [00 - 59]
    14059         'aggmonth': aggregation mulitple months:
    14060           yy1/[1,...,12]/01 00:00:00 - yym/[1,...,12]/01 00:00:00
    14061         'aggday': aggregation mulitple days:
    14062           yy1/mm1/[01,...,31,30,28/29] 00:00:00 - yym/mm1/[01,...,31,30,28/29] 00:00:00
    14063         'agghour': aggregation mulitple hours:
    14064           yy1/mm1/dd1 [00,...,23]:00:00 - yym/mmm/ddm [00,...,23]:59:59
    14065       amount: amount of periods
    14066       temp_desc: disctionary with the description of the time_steps (from
    14067         `temporal_desc' function))
    14068     """
    14069     import datetime as dt
    14070     fname = 'time_slices'
    14071 
    14072     availper = ['year', 'month', 'day', 'hour', 'minute', 'mmonth', 'mday', 'mhour']
    14073 
    14074     mattimes = CFtimesvar_datetime(tv, tu, cal)
    14075 
    14076     if type(tv) == type(np.arange(2)):
    14077         dimt = tv.shape[0]
    14078     elif type(tv) == type(range(2)):
    14079         dimt = len(tv)
    14080     else:
    14081         print errormsg
    14082         print '  ' + fname + ': time values type ', type(tv), 'not ready !!'
    14083         print '    available ones:', type(np.arange(1)), type(range(2))
    14084         quit(-1)
    14085     #print 'mat times _______'   
    14086     #for it in range(dimt):
    14087     #    print '  ',it, mattimes[it,:]
    14088 
    14089     # Beginning and Ending of time data
    14090     it = 0
    14091     idate = dt.datetime(mattimes[it,0], mattimes[it,1], mattimes[it,2],             \
    14092       mattimes[it,3], mattimes[it,4], mattimes[it,5])
    14093     it = dimt-1
    14094     edate = dt.datetime(mattimes[it,0], mattimes[it,1], mattimes[it,2],             \
    14095       mattimes[it,3], mattimes[it,4], mattimes[it,5])
    14096 
    14097     istats = mattimes[0,:]
    14098     estats = mattimes[dimt-1,:]
    14099     #print 'idate:', idate, 'edate:', edate
    14100 
    14101     # Getting the dates
    14102     if per == 'year':
    14103         # Arranging accordingly beginning/ending statistics periods
    14104         istats[1] = 1
    14105         istats[2] = 1
    14106         istats[3:5] = 0
    14107         estats[0] = estats[0] + 1
    14108         estats[1] = 1
    14109         estats[2] = 1
    14110         estats[3:5] = 0
    14111 
    14112         istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
    14113           istats[5])
    14114         estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
    14115           estats[5])
    14116 
    14117         iper = 0
    14118         # Getting statistics periods
    14119         stmatdates = {iper: [istats]}
    14120         # Single periods per stats
    14121         idate = istdate
    14122         imdate = istats.copy()
    14123         while idate < estdate:
    14124             iper = iper + 1
    14125             imdate = imdate.copy()
    14126             imdate[0] = imdate[0] + amount
    14127             idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
    14128               imdate[5])
    14129             if idate < estdate: stmatdates[iper] = [imdate]
    14130             else: break
    14131         stmatdates[iper] = [imdate]
    14132 
    14133     elif per == 'month':
    14134         # Arranging accordingly beginning/ending statistics periods
    14135         istats[2] = 1
    14136         istats[3:5] = 0
    14137         estats[1] = estats[1]+1
    14138         if estats[1] > 12:
    14139             estats[0] = estats[0] + 1
    14140             estats[1] = 1
    14141         estats[2] = 1
    14142         estats[3:5] = 0
    14143 
    14144         istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
    14145           istats[5])
    14146         estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
    14147           estats[5])
    14148 
    14149         iper = 0
    14150         # Getting statistics periods
    14151         stmatdates = {iper: [istats]}
    14152         idate = istdate
    14153         imdate = istats.copy()
    14154         while idate < estdate:
    14155             iper = iper + 1
    14156             imdate = imdate.copy()
    14157             imdate[1] = imdate[1] + amount
    14158             if imdate[1] > 12:
    14159                 imdate[0] = imdate[0] + 1
    14160                 imdate[1] = 1
    14161             idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
    14162               imdate[5])
    14163             if idate < estdate: stmatdates[iper] = [imdate]
    14164             else: break
    14165         stmatdates[iper] = [imdate]
    14166 
    14167     elif per == 'day':
    14168         # Arranging accordingly beginning/ending statistics periods
    14169         istats[3:5] = 0
    14170         iper = 0
    14171         estats[2] = estats[2]+1
    14172         if estats[2] > dmon:
    14173             estats[1] = estats[1] + 1
    14174             if estats[1] > 12:
    14175                 estats[0] = estats[0] + 1
    14176                 estats[1] = 1
    14177         estats[3:5] = 0
    14178 
    14179         istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
    14180           istats[5])
    14181         estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
    14182           estats[5])
    14183 
    14184         iper = 0
    14185         # Getting statistics periods
    14186         stmatdates = {iper: [istats]}
    14187         idate = istdate
    14188         imdate = istats.copy()
    14189         while idate < estdate:
    14190             iper = iper + 1
    14191             imdate = imdate.copy()
    14192             dmon = days_month(imdate[0], imdate[1])
    14193             imdate[2] = imdate[2] + amount
    14194             if imdate[2] > dmon:
    14195                 imdate[2] = imdate[2] - dmon
    14196                 imdate[1] = imdate[1] + 1
    14197                 if imdate[1] > 12:
    14198                     imdate[0] = imdate[0] + 1
    14199                     imdate[1] = 1
    14200             idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
    14201               imdate[5])
    14202             if idate < estdate: stmatdates[iper] = [imdate]
    14203             else: break
    14204         stmatdates[iper] = [imdate]
    14205 
    14206     elif per == 'aggmonth':
    14207         if temp_desc is None:
    14208             print errormsg
    14209             print '  ' + fname + ": For '" + per + "' temporal description is required"
    14210             quit(-1)
    14211 
    14212         amount = 1
    14213         # Arranging accordingly beginning/ending statistics periods
    14214         istats[2] = 1
    14215         istats[3:5] = 0
    14216         estats[1] = estats[1]+1
    14217         if estats[1] > 12:
    14218             estats[0] = estats[0] + 1
    14219             estats[1] = 1
    14220         estats[2] = 1
    14221         estats[3:5] = 0
    14222 
    14223         istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
    14224           istats[5])
    14225         estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
    14226           estats[5])
    14227 
    14228         iper = 0
    14229         # Getting statistics periods. We got 12 periods
    14230         stmatdates = {istats[1]: [istats]}
    14231         idate = istdate
    14232         imdate = istats.copy()
    14233         while idate < estdate:
    14234             imdate = imdate.copy()
    14235             imdate[1] = imdate[1] + amount
    14236             if imdate[1] > 12:
    14237                 imdate[0] = imdate[0] + 1
    14238                 imdate[1] = 1
    14239             idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
    14240               imdate[5])
    14241             if idate < estdate:
    14242                 if stmatdates.has_key(imdate[1]):
    14243                     vals = stmatdates[imdate[1]]
    14244                     vals = vals + [imdate]
    14245                     stmatdates[imdate[1]] = vals
    14246                 else:
    14247                     stmatdates[imdate[1]] = [imdate]
    14248             else: break
    14249         if stmatdates.has_key(imdate[1]):
    14250             vals = stmatdates[imdate[1]]
    14251             vals = vals + [imdate]
    14252             stmatdates[imdate[1]] = vals
    14253         else:
    14254             stmatdates[imdate[1]] = [imdate]
    14255 
    14256     else:
    14257         print errormsg
    14258         print '  ' + fname + ": period '" + per + "' not ready !!"
    14259         print '    available ones:', availper
    14260         quit(-1)
    14261 
    14262     # Getting the slices
    14263     ##
    14264     slices = []
    14265 
    14266     # Passing to CF-times
    14267     Srefdate0 = tu.split(' ')[2]
    14268     Sreftime0 = tu.split(' ')[3]
    14269     tunits = tu.split(' ')[0]
    14270 
    14271     if Srefdate0.find('/') != -1: Srdate = Srefdate0.replace('/', '')
    14272     elif Srefdate0.find('-') != -1: Srdate = Srefdate0.replace('-','')
    14273     else:
    14274         print errormsg
    14275         print '  ' + fname + ": I do not know how to create a [YYYY][MM][DD] " +     \
    14276           "date with '" + Srefdate0 + "' !!"
    14277         print "    expecting '/' or '-' splitting characters"
    14278         quit(-1)
    14279 
    14280     if Sreftime0.find(':') != -1: Srtime = Sreftime0.replace(':', '')
    14281     else:
    14282         print errormsg
    14283         print '  ' + fname + ": I do not know how to create a [HH][MI][SS] " +       \
    14284           "date with '" + Sreftime0 + "' !!"
    14285         print "    expecting ':' splitting character"
    14286         quit(-1)
    14287 
    14288     Srefdate = Srdate + Srtime
    14289     Nslices = len(stmatdates.keys())
    14290     stmdates = np.zeros((Nslices,6), dtype=np.int)
    14291     #print '  dates slices ________'
    14292     if per[0:3] == 'agg':
    14293         Lper=len(per)
    14294         pern = len[3:Lper]
    14295         print '  ' + fname + ": slices as aggregations of '" + pern + "' !!"
    14296     else:
    14297         for islc in range(Nslices):
    14298             #print islc, '  ', stmatdates[islc]
    14299             stmdates[islc,:] = stmatdates[islc]
    14300 
    14301         # CF-standard format of the dates of the slices
    14302         cfstdates = realdatetime_CFcompilant(stmdates, Srefdate, tunits)
    14303 
    14304         Nslices = len(cfstdates)
    14305         itt = 0
    14306         icfst = 0
    14307         for it in range(dimt):
    14308             timeslice = [itt,itt,1]
    14309             timeslice[1] = it
    14310             icfstup = np.min([icfst+1, Nslices-1])
    14311             if not(tv[it] >= cfstdates[icfst] and tv[it] < cfstdates[icfstup]):
    14312                 slices.append(timeslice)
    14313                 itt = it
    14314                 icfst = icfst + 1
    14315         slices.append(timeslice)
    14316 
    14317     Nslices = len(slices)
    14318     print infmsg
    14319     print '    ' + fname + ' slices _______'
    14320     for isl in range(Nslices-1):
    14321         print '    ', isl, stmatdates[isl], stmatdates[isl+1], ':', slices[isl]
    14322 
    14323     return slices, Nslices
    14324 
    14325 tv = []
    14326 values = []
    14327 totT = 24.*60.
    14328 dT = 7.345
    14329 for it in range(0,10):
    14330    for itt in range(3):
    14331         tv.append(29*24*60.+it*totT+dT*itt)
    14332         #values.append((it*3+itt)*1.)
    14333         values.append(it*1.)
    14334 
    14335 #print 'tv:', tv
    14336 #print 'values:', values
    14337 
    14338 vals = np.array(values)
    14339 itdim = 0
    14340 tu = 'minutes since 1949-12-01 00:00:00'
    14341 per = 'day'
    14342 calend = 'standard'
    14343 amount = 1
    14344 
    14345 tdes = temporal_desc(tv, tu, calend)
    14346 
    14347 tslc, Ntslc = time_slices(tv, tu, calend, per, amount)
    14348 
    14349 #statn = 'mean'
    14350 #tstat = np.ones((Ntslc), dtype=np.float)
    14351 #for islc in range(Ntslc):
    14352 #    timeslcv = tslc[islc]
    14353 #    timeslc = [slice(timeslcv[0],timeslcv[1],timeslcv[2])]
    14354 #    tvals = vals[tuple(timeslc)]
    14355 #    if statn == 'min':
    14356 #        tstat[islc] = np.min(tvals, axis=itdim)
    14357 #    elif statn == 'max':
    14358 #        tstat[islc] = np.max(tvals, axis=itdim)
    14359 #    elif statn == 'mean':
    14360 #        tstat[islc] = np.mean(tvals, axis=itdim)
    14361 #    elif statn == 'mean2':
    14362 #        tstat[islc] = np.mean(tvals**2, axis=itdim)
    14363 #    elif statn == 'std':
    14364 #        tstat[islc] = np.std(tvals, axis=itdim)
    14365 
    14366 #    timeslc = tslc[islc]
    14367 #    itime = tv[timeslc[0]]
    14368 #    etime = tv[timeslc[1]]
    14369 #    print islc, ':', itime, '-', etime, '<>', tvals, '=', tstat[islc]
    14370 #quit()
    1437114051
    1437214052def Latin_Greek(char):
    … …  
    1494914629        else:
    1495014630            tvals = {'DJF':[12,1,2], 'MAM':[3,4,5], 'JJA':[6,7,8], 'SON':[9,10,11]}
     14631            seasn = ['DJF', 'MAM', 'JJA', 'SON']
    1495114632            months12 = np.arange(12)+1
    1495214633            for it in range(dimt):
    14953                 for tvv in tvals:
    14954                     tvvv = tvals[tvv]
     14634                for tvv in range(4):
     14635                    tvvv = tvals[seasn[tvv]]
    1495514636                    iper = tvvv[0]
    1495614637                    eper = tvvv[2]
    14957                     if cyclevar_within(months12, iper, eper, mattimes[it,itsec]):
     14638                    if cyclevar_within(months12, iper, eper, mattimes[it,1]):
    1495814639                        if tvsec.has_key(tvv):
    1495914640                            vvv = tvsec[tvv]
    … …  
    1496214643                            vvv = [it]
    1496314644                        tvsec[tvv] = vvv
    14964                         break           
     14645                        break
     14646                print it, 'Lluis:', mattimes[it,:], tvv, mattimes[it,1]
     14647               
    1496514648        temporal_desc[timesec]=tvsec
    1496614649
    … …  
    1496914652    return temporal_desc
    1497014653
    14971 #quit()
    14972 
    14973 
     14654def time_slices(tv, tu, cal, per, amount, temp_desc):
     14655    """ Function to return temporal slices of a series of times for a given amount
     14656        of periods
     14657      tv: time vallues (CF format)
     14658      tu: CF time-units
     14659      cal: calendar
     14660      per: period
     14661        'year': full year [01/01 00:00:00 - 12/31 23:59:59]
     14662        'seas': full seasons 'DJF':[12,1,2], 'MAM':[3,4,5], 'JJA':[6,7,8],
     14663          'SON':[9,10,11] from [startSEAS/01 00:00:00 - endSEAS/[30,29,27/28] 23:59:59]
     14664        'month': full months [01 00:00:00 - [31,30,28/29] 23:59:59]
     14665        'day': full days [00:00:00 - 23:59:59]
     14666        'hour': full hours [00:00 - 59:59]
     14667        'minute': full minutes [00 - 59]
     14668        'aggmonth': aggregation mulitple months:
     14669          yy1/[1,...,12]/01 00:00:00 - yym/[1,...,12]/01 00:00:00
     14670        'aggday': aggregation mulitple days:
     14671          yy1/mm1/[01,...,31,30,28/29] 00:00:00 - yym/mm1/[01,...,31,30,28/29] 00:00:00
     14672        'agghour': aggregation mulitple hours:
     14673          yy1/mm1/dd1 [00,...,23]:00:00 - yym/mmm/ddm [00,...,23]:59:59
     14674      amount: amount of periods
     14675      temp_desc: disctionary with the description of the time_steps (from
     14676        `temporal_desc' function)
     14677    """
     14678    import datetime as dt
     14679    fname = 'time_slices'
     14680
     14681    seasons = {0:[12,1,2], 1:[3,4,5], 2:[6,7,8], 3:[9,10,11]}
     14682    months12 = np.arange(12)+1
     14683
     14684    availper = ['year', 'season', 'month', 'day', 'hour', 'minute', 'aggseason',     \
     14685      'aggmonth', 'aggday', 'agghour']
     14686
     14687    if not searchInlist(availper, per):
     14688        print errormsg
     14689        print '  ' + fname + ": period '" + per + "' not ready !!"
     14690        print '    available ones:', availper
     14691        quit(-1)
     14692
     14693
     14694    mattimes = CFtimesvar_datetime(tv, tu, cal)
     14695
     14696    if type(tv) == type(np.arange(2)):
     14697        dimt = tv.shape[0]
     14698    elif type(tv) == type(range(2)):
     14699        dimt = len(tv)
     14700    else:
     14701        print errormsg
     14702        print '  ' + fname + ': time values type ', type(tv), 'not ready !!'
     14703        print '    available ones:', type(np.arange(1)), type(range(2))
     14704        quit(-1)
     14705    #print 'mat times _______'   
     14706    #for it in range(dimt):
     14707    #    print '  ',it, mattimes[it,:]
     14708
     14709    # Beginning and Ending of time data
     14710    it = 0
     14711    idate = dt.datetime(mattimes[it,0], mattimes[it,1], mattimes[it,2],             \
     14712      mattimes[it,3], mattimes[it,4], mattimes[it,5])
     14713    it = dimt-1
     14714    edate = dt.datetime(mattimes[it,0], mattimes[it,1], mattimes[it,2],             \
     14715      mattimes[it,3], mattimes[it,4], mattimes[it,5])
     14716
     14717    istats = mattimes[0,:]
     14718    estats = mattimes[dimt-1,:]
     14719    #print 'idate:', idate, 'edate:', edate
     14720
     14721    # Getting the dates
     14722    if per == 'year':
     14723        # Arranging accordingly beginning/ending statistics periods
     14724        istats[1] = 1
     14725        istats[2] = 1
     14726        istats[3:5] = 0
     14727        estats[0] = estats[0] + 1
     14728        estats[1] = 1
     14729        estats[2] = 1
     14730        estats[3:5] = 0
     14731
     14732        istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
     14733          istats[5])
     14734        estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
     14735          estats[5])
     14736
     14737        iper = 0
     14738        # Getting statistics periods
     14739        stmatdates = {iper: [istats]}
     14740        # Single periods per stats
     14741        idate = istdate
     14742        imdate = istats.copy()
     14743        while idate < estdate:
     14744            iper = iper + 1
     14745            imdate = imdate.copy()
     14746            imdate[0] = imdate[0] + amount
     14747            idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
     14748              imdate[5])
     14749            if idate < estdate: stmatdates[iper] = [imdate]
     14750            else: break
     14751        stmatdates[iper] = [imdate]
     14752
     14753    elif per == 'season':
     14754        # Arranging accordingly beginning/ending statistics periods
     14755        istats[2] = 1
     14756        istats[3:5] = 0
     14757        for iseas in range(4):
     14758            vseas = seasons[iseas]
     14759            bseas = vseas[0]
     14760            eseas = vseas[2]
     14761            if cyclevar_within(months12, bseas, eseas, istats[1]):
     14762                istats[1] = bseas
     14763
     14764        mone = estats[1]
     14765        for iseas in range(4):
     14766            vseas = seasons[iseas]
     14767            bseas = vseas[0]
     14768            eseas = vseas[2]
     14769            if cyclevar_within(months12, bseas, eseas, estats[1]):
     14770                estats[1] = bseas
     14771                bereak
     14772        # Checking if the end of season is on the following year...
     14773        if bseas < mone: estats[0] = estats[0] + 1
     14774
     14775        estats[2] = days_month(estats[0],estats[1])
     14776        estats[3] = 23
     14777        estats[4] = 59
     14778        estats[5] = 59
     14779
     14780        istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
     14781          istats[5])
     14782        estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
     14783          estats[5])
     14784
     14785        iper = 0
     14786        stmatdates = {}
     14787        # Getting statistics periods
     14788        for iseas in range(4):
     14789            vseas = seasons[iseas]
     14790            bseas = vseas[0]
     14791            eseas = vseas[2]
     14792            if cyclevar_within(months12, bseas, eseas, istats[1]):
     14793                lval = list(istats)
     14794                stmatdates[iper] = lval
     14795                break           
     14796
     14797        # Single periods per stats
     14798        idate = istdate
     14799        imdate = istats.copy()
     14800        while idate < estdate:
     14801            iper = iper + 1
     14802            imdate[1] = imdate[1] + 3
     14803            if imdate[1] > 12:
     14804                imdate[0] = imdate[0] + 1
     14805                imdate[1] = imdate[1] - 12
     14806            idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
     14807              imdate[5])
     14808            for iseas in range(4):
     14809                vseas = seasons[iseas]
     14810                bseas = vseas[0]
     14811                eseas = vseas[2]
     14812                if cyclevar_within(months12, bseas, eseas, imdate[1]):
     14813                    lval = list(imdate)
     14814                    stmatdates[iper] = lval
     14815                    break
     14816        lval = list(imdate)
     14817        stmatdates[iper] = lval
     14818        printing_dictionary(stmatdates)
     14819
     14820    elif per == 'month':
     14821        # Arranging accordingly beginning/ending statistics periods
     14822        istats[2] = 1
     14823        istats[3:5] = 0
     14824        estats[1] = estats[1]+1
     14825        if estats[1] > 12:
     14826            estats[0] = estats[0] + 1
     14827            estats[1] = 1
     14828        estats[2] = 1
     14829        estats[3:5] = 0
     14830
     14831        istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
     14832          istats[5])
     14833        estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
     14834          estats[5])
     14835
     14836        iper = 0
     14837        # Getting statistics periods
     14838        stmatdates = {iper: [istats]}
     14839        idate = istdate
     14840        imdate = istats.copy()
     14841        while idate < estdate:
     14842            iper = iper + 1
     14843            imdate = imdate.copy()
     14844            imdate[1] = imdate[1] + amount
     14845            if imdate[1] > 12:
     14846                imdate[0] = imdate[0] + 1
     14847                imdate[1] = 1
     14848            idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
     14849              imdate[5])
     14850            if idate < estdate: stmatdates[iper] = [imdate]
     14851            else: break
     14852        stmatdates[iper] = [imdate]
     14853
     14854    elif per == 'day':
     14855        # Arranging accordingly beginning/ending statistics periods
     14856        istats[3:5] = 0
     14857        iper = 0
     14858        estats[2] = estats[2]+1
     14859        dmon = days_month(estats[0], estats[1])
     14860        if estats[2] > dmon:
     14861            estats[1] = estats[1] + 1
     14862            if estats[1] > 12:
     14863                estats[0] = estats[0] + 1
     14864                estats[1] = 1
     14865        estats[3:5] = 0
     14866
     14867        istdate = dt.datetime(istats[0], istats[1], istats[2], istats[3], istats[4], \
     14868          istats[5])
     14869        estdate = dt.datetime(estats[0], estats[1], estats[2], estats[3], estats[4], \
     14870          estats[5])
     14871
     14872        iper = 0
     14873        # Getting statistics periods
     14874        stmatdates = {iper: [istats]}
     14875        idate = istdate
     14876        imdate = istats.copy()
     14877        while idate < estdate:
     14878            iper = iper + 1
     14879            imdate = imdate.copy()
     14880            dmon = days_month(imdate[0], imdate[1])
     14881            imdate[2] = imdate[2] + amount
     14882            if imdate[2] > dmon:
     14883                imdate[2] = imdate[2] - dmon
     14884                imdate[1] = imdate[1] + 1
     14885                if imdate[1] > 12:
     14886                    imdate[0] = imdate[0] + 1
     14887                    imdate[1] = 1
     14888            idate= dt.datetime(imdate[0], imdate[1], imdate[2], imdate[3], imdate[4],\
     14889              imdate[5])
     14890            if idate < estdate: stmatdates[iper] = [imdate]
     14891            else: break
     14892        stmatdates[iper] = [imdate]
     14893
     14894    elif per[0:3] == 'agg':
     14895        Lper = len(per)
     14896        pern = per[3:Lper]
     14897        print "  " + fname + ": temporal '" + pern + "' aggregation statistics !"
     14898
     14899    # Getting the slices
     14900    ##
     14901    slices = []
     14902
     14903    # Passing to CF-times
     14904    Srefdate0 = tu.split(' ')[2]
     14905    Sreftime0 = tu.split(' ')[3]
     14906    tunits = tu.split(' ')[0]
     14907
     14908    if Srefdate0.find('/') != -1: Srdate = Srefdate0.replace('/', '')
     14909    elif Srefdate0.find('-') != -1: Srdate = Srefdate0.replace('-','')
     14910    else:
     14911        print errormsg
     14912        print '  ' + fname + ": I do not know how to create a [YYYY][MM][DD] " +     \
     14913          "date with '" + Srefdate0 + "' !!"
     14914        print "    expecting '/' or '-' splitting characters"
     14915        quit(-1)
     14916
     14917    if Sreftime0.find(':') != -1: Srtime = Sreftime0.replace(':', '')
     14918    else:
     14919        print errormsg
     14920        print '  ' + fname + ": I do not know how to create a [HH][MI][SS] " +       \
     14921          "date with '" + Sreftime0 + "' !!"
     14922        print "    expecting ':' splitting character"
     14923        quit(-1)
     14924
     14925    Srefdate = Srdate + Srtime
     14926    #print '  dates slices ________'
     14927    if per[0:3] == 'agg':
     14928        Lper=len(per)
     14929        pern = per[3:Lper]
     14930        vtemp_desc = temp_desc[pern]
     14931        print vtemp_desc
     14932        vvtdesc = list(vtemp_desc.keys())
     14933        vvtdesc.sort()
     14934
     14935        Nslices = len(vtemp_desc)
     14936        slices = []
     14937        for islc in range(Nslices):
     14938            vvtemp_desc = vtemp_desc[vvtdesc[islc]]
     14939            print 'islc,', islc, 'vtemp:', vvtemp_desc
     14940            Ltt = len(vvtemp_desc)
     14941            slices.append(vvtemp_desc)
     14942       
     14943    else:
     14944        Nslices = len(stmatdates.keys())
     14945        stmdates = np.zeros((Nslices,6), dtype=np.int)
     14946        for islc in range(Nslices):
     14947            lval = stmatdates[islc]
     14948            stmdates[islc,:] = np.array(lval)
     14949
     14950        # CF-standard format of the dates of the slices
     14951        cfstdates = realdatetime_CFcompilant(stmdates, Srefdate, tunits)
     14952
     14953        Nslices = len(cfstdates)
     14954        itt = 0
     14955        icfst = 0
     14956        for it in range(dimt):
     14957            timeslice = [itt,itt,1]
     14958            timeslice[1] = it
     14959            icfstup = np.min([icfst+1, Nslices-1])
     14960            if not(tv[it] >= cfstdates[icfst] and tv[it] < cfstdates[icfstup]):
     14961                slices.append(timeslice)
     14962                itt = it
     14963                icfst = icfst + 1
     14964        slices.append(timeslice)
     14965
     14966        Nslices = len(slices)
     14967        print infmsg
     14968        print '    ' + fname + ':', Nslices, ' slices _______'
     14969        for isl in range(Nslices-1):
     14970            print '    ', isl, stmatdates[isl], stmatdates[isl+1], ':', slices[isl]
     14971        print '    ', isl+1, stmatdates[isl+1], stmatdates[isl+2], ':', slices[isl+1]
     14972
     14973    return slices, Nslices
     14974
     14975tv = []
     14976values = []
     14977totT = 24.
     14978dT = 7.345
     14979iTT = 29.
     14980for it in range(0,10):
     14981   for itt in range(3):
     14982        tv.append(iTT+it*totT+dT*itt)
     14983        #values.append((it*3+itt)*1.)
     14984        values.append(it*1.)
     14985
     14986print 'tv:', tv
     14987print 'values:', values
     14988
     14989vals = np.array(values)
     14990itdim = 0
     14991tu = 'days since 1949-12-01 00:00:00'
     14992per = 'aggseason'
     14993calend = 'standard'
     14994amount = 1
     14995
     14996tdes = temporal_desc(tv, tu, calend)
     14997
     14998tslc, Ntslc = time_slices(tv, tu, calend, per, amount, tdes)
     14999
     15000statn = 'mean'
     15001tstat = np.ones((Ntslc), dtype=np.float)
     15002for islc in range(Ntslc):
     15003    timeslcv = tslc[islc]
     15004    if per[0:3] != 'agg':
     15005        timeslc = [slice(timeslcv[0],timeslcv[1],timeslcv[2])]
     15006        tvals = vals[tuple(timeslc)]
     15007    else:
     15008        newshape = list(vals.shape)
     15009        aggslc = tslc[islc]
     15010        NNtslc = len(aggslc)
     15011        newshape[itdim] = NNtslc
     15012        tvals = np.zeros(tuple(newshape), dtype=vals.dtype)
     15013        origshape = list(vals.shape)
     15014        iishape = list(vals.shape)
     15015        for iislc in range(NNtslc):
     15016            origshape[itdim] = aggslc[iislc]
     15017            iishape[itdim] = iislc
     15018            print iislc, ':', iishape, '<>', origshape
     15019            tvals[tuple(iishape)] = vals[tuple(origshape)]
     15020
     15021    if statn == 'min':
     15022        tstat[islc] = np.min(tvals, axis=itdim)
     15023    elif statn == 'max':
     15024        tstat[islc] = np.max(tvals, axis=itdim)
     15025    elif statn == 'mean':
     15026        tstat[islc] = np.mean(tvals, axis=itdim)
     15027    elif statn == 'mean2':
     15028        tstat[islc] = np.mean(tvals**2, axis=itdim)
     15029    elif statn == 'std':
     15030        tstat[islc] = np.std(tvals, axis=itdim)
     15031
     15032    timeslc = tslc[islc]
     15033    itime = tv[timeslc[0]]
     15034    etime = tv[timeslc[1]]
     15035    print islc, ':', itime, '-', etime, '<>', tvals, '=', tstat[islc]
     15036quit()
     15037
     15038
Note: See TracChangeset for help on using the changeset viewer.