![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
- Timestamp:
- Mar 12, 2019, 5:14:12 PM (6 years ago)
- Location:
- trunk/tools
- Files:
-
- 3 edited
-
generic_tools.py (modified) (4 diffs)
-
nc_var.py (modified) (1 diff)
-
nc_var_tools.py (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/generic_tools.py
r2403 r2404 14758 14758 14759 14759 tvsec = {} 14760 if timesec != 'season': 14760 if timesec == 'day': 14761 # This must be julian day without Feb 29 14762 dimt = mattimes.shape[0] 14763 tvals = range(dimt) 14761 14764 for it in range(dimt): 14762 for tvv in tvals: 14763 if mattimes[it,itsec] == tvv: 14765 if mattimes[it,1] == 2 and mattimes[it,2] == 29: 14766 print ' ' + fname + ": for time section '" + timesce + \ 14767 "' no leap date February 29 !!" 14768 else: 14769 date= dt.datetime(mattimes[it,0], mattimes[it,1], mattimes[it,2],\ 14770 mattimes[it,3], mattimes[it,4], mattimes[it,5]) 14771 julday = int(date.strftime("%j")) 14772 for tvv in range(365): 14773 if julday == tvv: 14764 14774 if tvsec.has_key(tvv): 14765 14775 vvv = tvsec[tvv] … … 14769 14779 tvsec[tvv] = vvv 14770 14780 break 14771 else: 14781 14782 elif timesec == 'season': 14772 14783 tvals = {'DJF':[12,1,2], 'MAM':[3,4,5], 'JJA':[6,7,8], 'SON':[9,10,11]} 14773 14784 seasn = ['DJF', 'MAM', 'JJA', 'SON'] … … 14779 14790 eper = tvvv[2] 14780 14791 if cyclevar_within(months12, iper, eper, mattimes[it,1]): 14792 if tvsec.has_key(tvv): 14793 vvv = tvsec[tvv] 14794 vvv.append(it) 14795 else: 14796 vvv = [it] 14797 tvsec[tvv] = vvv 14798 break 14799 else: 14800 for it in range(dimt): 14801 for tvv in tvals: 14802 if mattimes[it,itsec] == tvv: 14781 14803 if tvsec.has_key(tvv): 14782 14804 vvv = tvsec[tvv] … … 15092 15114 Naggv = aggvals.shape[0] 15093 15115 slices = [] 15094 for iagg in range(Naggv): 15095 iaggv = aggvals[iagg] 15096 found = False 15097 for islc in range(Nslices): 15098 if vvtdesc[islc] == iaggv: 15099 vvtemp_desc = vtemp_desc[vvtdesc[islc]] 15100 slices.append(vvtemp_desc) 15101 found = True 15102 break 15103 if not found: slices.append(None) 15116 if pern != 'day': 15117 for iagg in range(Naggv): 15118 iaggv = aggvals[iagg] 15119 found = False 15120 for islc in range(Nslices): 15121 if vvtdesc[islc] == iaggv: 15122 vvtemp_desc = vtemp_desc[vvtdesc[islc]] 15123 slices.append(vvtemp_desc) 15124 found = True 15125 break 15126 if not found: slices.append(None) 15127 else: 15128 # We do not want 29 Feb. 15129 for iagg in range(Naggv): 15130 iaggv = aggvals[iagg] 15131 found = False 15132 for islc in range(Nslices): 15133 if vvtdesc[islc] == iaggv: 15134 vvtemp_desc = vtemp_desc[vvtdesc[islc]] 15135 slices.append(vvtemp_desc) 15136 found = True 15137 break 15138 if not found: slices.append(None) 15104 15139 Nslices = Naggv + 0 15105 15140 else: -
trunk/tools/nc_var.py
r2384 r2404 61 61 ## e.g. # nc_var.py -o curve_section -f /home/lluis/PY/test.nc -S 'gridline,x,y,8.,8.,16.,16.,32' -v all 62 62 ## e.g. # nc_var.py -o merge_files -S 'plev|plev,time|time:merged.nc' -f 'ncobs/AliceSprings/snd_94326_198201-198201.nc,ncobs/AliceSprings/snd_94326_198202-198202.nc' -v 'plev,time' 63 ## e.g. # nc_var.py -o temporal_stats -S 'Time:WRFtime:day@1@min,LTday@-3@1@min :bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -f wrfout_d01_1995-01-01_00:00:00 -v T2,Q263 ## e.g. # nc_var.py -o temporal_stats -S 'Time:WRFtime:day@1@min,LTday@-3@1@min,agghour@1@mean:bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -f wrfout_d01_1995-01-01_00:00:00 -v T2,Q2 64 64 ## e.g. # nc_var.py -o retrieve_stations -f wrfout_d01_1995-01-01_00:00:00 -S 'tmin_percentiles.nc:stname:None:stlon:stlat:None:nearest:west_east:XLONG:south_north:XLAT:HGT:Time:WRFtime' -v T2,QVAPOR 65 65 ## e.g. # nc_var.py -o compute_slice2Dstats -S 'XLAT,-63.,19.,2.,HGT,250.,7000.,500.,Time|Times:west_east|XLONG:south_north|XLAT' -f wrfout_d01_1995-01-01_00:00:00 -v T2,Q2 -
trunk/tools/nc_var_tools.py
r2402 r2404 26132 26132 only that ones with [timedimn]) 26133 26133 """ 26134 import datetime as dt 26134 26135 fname = 'temporal_stats' 26135 26136 … … 26288 26289 statslices[iTst] = [Ntimeslice, timeslice] 26289 26290 26290 print 'Lluis slices period:', period, 'amount:', amount, ' ________'26291 for it in range(Ntimeslice):26292 print it, ':', timeslice[it]26291 #print 'Lluis slices period:', period, 'amount:', amount, ' ________' 26292 #for it in range(Ntimeslice): 26293 # print it, ':', timeslice[it] 26293 26294 26294 26295 # Creation of output file … … 26387 26388 PerN = periodn[3:len(periodn)] 26388 26389 if PerN == 'year': 26390 iit = 0 26389 26391 dtsec = 365 * 24 * 3600 26390 26392 elif PerN == 'month': 26393 iit = 1 26391 26394 dtsec = 365 * 24 * 3600 26392 26395 elif PerN == 'day': 26396 tuv = timeu.split(' ') 26397 if len(tuv) == 4: 26398 refD= tuv[2] + ' ' + tuv[3] 26399 elif len(tuv) == 3: 26400 refD= tuv[2] + ' 00:00:00' 26401 yr = int(refD[0:4]) 26402 mm = int(refD[5:7]) 26403 dd = int(refD[8:10]) 26404 hr = int(refD[11:13]) 26405 mn = int(refD[14:16]) 26406 ss = int(refD[17:19]) 26407 26408 refdateu = dt.datetime(yr, mm, dd, hr, mn, ss) 26409 idate = dt.datetime(yr+1, 1, 1, 0, 0, 0) 26410 26411 diffd = idate - refdateu 26412 if gen.searchInlist(dir(diffd), 'total_seconds'): 26413 iit = diffd.total_seconds()/(24.*3600.) 26414 else: 26415 iit = diffd.days/1. + diffd.seconds/(24.*3600.) 26393 26416 dtsec = 24 * 3600 26394 26417 elif PerN == 'hour': 26418 iit = 0 26395 26419 dtsec = 3600 26396 26420 elif PerN == 'minute': 26421 iit = 0 26397 26422 dtsec = 60 26398 26423 elif PerN == 'second': 26424 iit = 0 26399 26425 dtsec = 1 26400 26426 else: … … 26407 26433 # Getting slices in axis-time units 26408 26434 ut = timeu.split(' ')[0] 26409 if ut == 'year': dtsec / (365 * 24 * 3600.) 26410 elif ut == 'month': dtsec / (365 * 24 * 3600.) 26411 elif ut == 'day': dtsec / (24 * 3600.) 26412 elif ut == 'hour': dtsec / (3600.) 26413 elif ut == 'minute': dtsec / (60.) 26414 elif ut == 'second': dtsec / 1. 26435 if ut == 'days': 26436 ddt = dtsec / (24 * 3600.) 26437 elif ut == 'hours': 26438 ddt = dtsec / (3600.) 26439 elif ut == 'minutes': 26440 ddt = dtsec / (60.) 26441 elif ut == 'seconds': 26442 ddt = dtsec / 1. 26443 else: 26444 print errormsg 26445 print ' ' + fname + ": temporal units '" + ut + "' not ready !!" 26446 print ' available ones:', availcftunits 26447 quit(-1) 26415 26448 26416 26449 timeslcev = tslcs[0] … … 26424 26457 ' for ' + periodn, timeu) 26425 26458 newvarb.setncattr('calendar', timec) 26426 for it in range(Ntslc -1):26459 for it in range(Ntslc): 26427 26460 timeslcev = tslcs[it] 26428 Nvalsslc = len(timeslcev) 26429 print 'timeslcev:', tslcs[it], 'Nvals:', Nvalsslc 26430 islc = newvart[it]+dtsec/2 26431 timeslcev = tslcs[it+1] 26432 eslc = newvart[it]-dtsec/2 26433 print it, ' Lluis islc:', islc, 'eslc:', eslc 26434 newvart[it] = dtsec 26461 if timeslcev is not None: 26462 Nvalsslc = len(timeslcev) 26463 if iit == 1: 26464 newvart[it] = iit*ddt/2.+it*ddt 26465 else: 26466 newvart[it] = (iit+it)*ddt+ddt/2. 26467 26468 islc = newvart[it]-ddt/2 26469 eslc = newvart[it]+ddt/2 26435 26470 newvarb[it,:] = [islc, eslc] 26436 26471 onewnc.sync() … … 26522 26557 for islc in range(Ntslc): 26523 26558 timeslcv = tslcs[islc] 26524 newvarN[islc] = timeslcv[1]-timeslcv[0]+126525 26559 26526 26560 if periodn[0:3] != 'agg': 26561 newvarN[islc] = timeslcv[1]-timeslcv[0]+1 26562 aggslc = 1 26527 26563 timeslc = [] 26528 26564 shapetimeslc = [] … … 26540 26576 tvals = ovar[tuple(timeslc)] 26541 26577 else: 26542 newshape = list(ovar.shape)26578 # Length for aggregation 26543 26579 aggslc = tslcs[islc] 26544 NNtslc = len(aggslc) 26545 newshape[itdim] = NNtslc 26546 tvals = np.zeros(tuple(newshape), dtype=ovar.dtype) 26547 origshape = list(ovar.shape) 26548 iishape = list(ovar.shape) 26549 for iislc in range(NNtslc): 26550 origshape[itdim] = aggslc[iislc] 26551 iishape[itdim] = iislc 26552 tvals[tuple(iishape)] = ovar[tuple(origshape)] 26553 26554 if statn == 'min': 26555 newvarv[tuple(tstatslc)] = np.min(tvals, axis=itdim) 26556 elif statn == 'max': 26557 newvarv[tuple(tstatslc)] = np.max(tvals, axis=itdim) 26558 elif statn == 'mean': 26559 newvarv[tuple(tstatslc)] = np.mean(tvals, axis=itdim) 26560 elif statn == 'mean2': 26561 newvarv[tuple(tstatslc)] = np.mean(tvals**2, axis=itdim) 26562 elif statn == 'std': 26563 newvarv[tuple(tstatslc)] = np.std(tvals, axis=itdim) 26564 elif statn == 'percentiles': 26565 axesS = str(len(tvals.shape)-itdim) 26566 if tvals.shape[itdim] < Npercen: 26567 tstatslc = [slice(0,Npercen+1)] + tstatslc 26568 print warnmsg 26569 print ' ' + fname + ': Not enough temporal data:', \ 26570 tvals.shape[itdim], ' to copmute:',Npercen, 'percentiles !!' 26571 dimspercen = [Npercen+1] 26572 for i in range(len(tvals.shape)): 26573 if i != itdim: 26574 dimspercen.append(tvals.shape[i]) 26575 percens = np.ones(tuple(dimspercen), dtype=np.float)*gen.fillValueF 26576 newvarv[tuple(tstatslc)] = percens[:] 26580 if aggslc is not None: 26581 NNtslc = len(aggslc) 26582 26583 newvarN[islc] = NNtslc 26584 timeslc = [] 26585 shapetimeslc = [] 26586 tstatslc = [] 26587 for dn in vardims: 26588 if dn == timedimn: 26589 timeslc.append(NNtslc) 26590 shapetimeslc.append(NNtslc) 26591 tstatslc.append(islc) 26592 else: 26593 Ldim = len(onc.dimensions[dn]) 26594 timeslc.append(slice(0,Ldim)) 26595 shapetimeslc.append(Ldim) 26596 tstatslc.append(slice(0,Ldim)) 26597 26598 tvals = np.zeros(tuple(shapetimeslc), dtype=ovar.dtype) 26599 # Getting slices 26600 for iagg in range(NNtslc): 26601 timeslc[itdim] = aggslc[iagg] 26602 aggtimeslc = timeslc + [] 26603 aggtimeslc[itdim] = iagg 26604 tvals[tuple(aggtimeslc)] = ovar[tuple(timeslc)] 26577 26605 else: 26578 tvalst = tvals.transpose() 26579 if len(tvals.shape) == 4: 26580 percenst = fsci.module_scientific.percentiles_r_k4d(tvalst, \ 26581 axesS, npercen=Npercen+1, d1=tvals.shape[3], \ 26582 d2=tvals.shape[2], d3=tvals.shape[1], d4=tvals.shape[0]) 26583 elif len(tvals.shape) == 3: 26584 percenst = fsci.module_scientific.percentiles_r_k3d(tvalst, \ 26585 axesS, npercen=Npercen+1, d1=tvals.shape[2], \ 26586 d2=tvals.shape[1], d3=tvals.shape[0]) 26587 elif len(tvals.shape) == 2: 26588 percenst = fsci.module_scientific.percentiles_r_k2d(tvalst, \ 26589 axesS, npercen=Npercen+1, d1=tvals.shape[1], \ 26590 d2=tvals.shape[0]) 26591 elif len(tvals.shape) == 1: 26592 quants = Quantiles(tvals, Npercen) 26593 percenst = quants.quantilesv 26606 newvarN[islc] = 0 26607 timeslc = [] 26608 shapetimeslc = [] 26609 tstatslc = [] 26610 for dn in vardims: 26611 if dn == timedimn: 26612 tstatslc.append(islc) 26613 else: 26614 tstatslc.append(slice(0,Ldim)) 26615 if aggslc is not None: 26616 if statn == 'min': 26617 newvarv[tuple(tstatslc)] = np.min(tvals, axis=itdim) 26618 elif statn == 'max': 26619 newvarv[tuple(tstatslc)] = np.max(tvals, axis=itdim) 26620 elif statn == 'mean': 26621 newvarv[tuple(tstatslc)] = np.mean(tvals, axis=itdim) 26622 elif statn == 'mean2': 26623 newvarv[tuple(tstatslc)] = np.mean(tvals**2, axis=itdim) 26624 elif statn == 'std': 26625 newvarv[tuple(tstatslc)] = np.std(tvals, axis=itdim) 26626 elif statn == 'percentiles': 26627 axesS = str(len(tvals.shape)-itdim) 26628 if tvals.shape[itdim] < Npercen: 26629 tstatslc = [slice(0,Npercen+1)] + tstatslc 26630 print warnmsg 26631 print ' ' + fname + ': Not enough temporal data:', \ 26632 tvals.shape[itdim], ' to copmute:',Npercen, 'percentiles !!' 26633 dimspercen = [Npercen+1] 26634 for i in range(len(tvals.shape)): 26635 if i != itdim: 26636 dimspercen.append(tvals.shape[i]) 26637 percens = np.ones(tuple(dimspercen), dtype=np.float)*gen.fillValueF 26638 newvarv[tuple(tstatslc)] = percens[:] 26594 26639 else: 26595 print errormsg 26596 print ' ' +fname+ ': rank of temporal data:', tvals.shape, \ 26597 'not ready !!' 26598 quit(-1) 26599 percens = percenst.transpose() 26600 26601 shapev = list(percens.shape) 26602 for idd in range(len(shapev)): 26603 if idd != itdim+1: 26604 shapev[idd] = slice(0,shapev[idd]) 26640 tvalst = tvals.transpose() 26641 if len(tvals.shape) == 4: 26642 percenst = fsci.module_scientific.percentiles_r_k4d(tvalst, \ 26643 axesS, npercen=Npercen+1, d1=tvals.shape[3], \ 26644 d2=tvals.shape[2], d3=tvals.shape[1], d4=tvals.shape[0]) 26645 elif len(tvals.shape) == 3: 26646 percenst = fsci.module_scientific.percentiles_r_k3d(tvalst, \ 26647 axesS, npercen=Npercen+1, d1=tvals.shape[2], \ 26648 d2=tvals.shape[1], d3=tvals.shape[0]) 26649 elif len(tvals.shape) == 2: 26650 percenst = fsci.module_scientific.percentiles_r_k2d(tvalst, \ 26651 axesS, npercen=Npercen+1, d1=tvals.shape[1], \ 26652 d2=tvals.shape[0]) 26653 elif len(tvals.shape) == 1: 26654 quants = Quantiles(tvals, Npercen) 26655 percenst = quants.quantilesv 26605 26656 else: 26606 shapev[itdim+1] = 0 26607 26608 tstatslc = [slice(0,Npercen+1)] + tstatslc 26609 newvarv[tuple(tstatslc)] = percens[tuple(shapev)] 26657 print errormsg 26658 print ' ' +fname+ ': rank of temporal data:', tvals.shape, \ 26659 'not ready !!' 26660 quit(-1) 26661 percens = percenst.transpose() 26662 26663 shapev = list(percens.shape) 26664 for idd in range(len(shapev)): 26665 if idd != itdim+1: 26666 shapev[idd] = slice(0,shapev[idd]) 26667 else: 26668 shapev[itdim+1] = 0 26669 26670 tstatslc = [slice(0,Npercen+1)] + tstatslc 26671 newvarv[tuple(tstatslc)] = percens[tuple(shapev)] 26672 else: 26673 print errormsg 26674 print ' ' + fname + ": statisitcs '" + statn + "' not ready !!" 26675 statnor = Lstatn.keys() 26676 statnor.sort() 26677 print ' available ones:', statnor 26678 quit(-1) 26679 onewnc.sync() 26610 26680 else: 26611 print errormsg 26612 print ' ' + fname + ": statisitcs '" + statn + "' not ready !!" 26613 statnor = Lstatn.keys() 26614 statnor.sort() 26615 print ' available ones:', statnor 26616 quit(-1) 26617 onewnc.sync() 26618 26681 newvarv[tuple(tstatslc)] = gen.fillValueF 26619 26682 #newvar[:] = tstat 26620 26683 onewnc.sync()
Note: See TracChangeset
for help on using the changeset viewer.
