Repository: climate Updated Branches: refs/heads/master 15a93c833 -> 9183f5fea
removed process.py dependency and fixed Nio toNETCDF logical errors Project: http://git-wip-us.apache.org/repos/asf/climate/repo Commit: http://git-wip-us.apache.org/repos/asf/climate/commit/48352bc0 Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/48352bc0 Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/48352bc0 Branch: refs/heads/master Commit: 48352bc04c4fe4a1e8a03a254ecbc659424ad147 Parents: 6d47a57 Author: Kim Whitehall <[email protected]> Authored: Fri Oct 24 15:23:29 2014 -0700 Committer: Kim Whitehall <[email protected]> Committed: Fri Oct 24 15:23:29 2014 -0700 ---------------------------------------------------------------------- mccsearch/code/mccSearch.py | 246 +++++++++++++++++++++++++++++++++---- mccsearch/code/mccSearchUI.py | 2 +- 2 files changed, 222 insertions(+), 26 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/climate/blob/48352bc0/mccsearch/code/mccSearch.py ---------------------------------------------------------------------- diff --git a/mccsearch/code/mccSearch.py b/mccsearch/code/mccSearch.py index e691ed5..7bfaed6 100644 --- a/mccsearch/code/mccSearch.py +++ b/mccsearch/code/mccSearch.py @@ -44,8 +44,9 @@ import matplotlib.colors as mcolors from matplotlib.ticker import FuncFormatter, FormatStrFormatter #existing modules in services -import files -import process +import Nio +#import files +#import process #----------------------- GLOBAL VARIABLES -------------------------- # --------------------- User defined variables --------------------- #FYI the lat lon values are not necessarily inclusive of the points given. These are the limits @@ -53,7 +54,7 @@ import process LATMIN = '5.0' #min latitude; -ve values in the SH e.g. 5S = -5 LATMAX = '19.0' #max latitude; -ve values in the SH e.g. 5S = -5 20.0 LONMIN = '-5.0' #min longitude; -ve values in the WH e.g. 59.8W = -59.8 -30 -LONMAX = '9.0' #min longitude; -ve values in the WH e.g. 59.8W = -59.8 30 +LONMAX = '5.0' #min longitude; -ve values in the WH e.g. 59.8W = -59.8 30 XRES = 4.0 #x direction spatial resolution in km YRES = 4.0 #y direction spatial resolution in km TRES = 1 #temporal resolution in hrs @@ -139,14 +140,25 @@ def readMergData(dirname, filelist = None): sys.exit() else: # Open the first file in the list to read in lats, lons and generate the grid for comparison - #tmp = Nio.open_file(filelist[0], format='nc') - tmp = Dataset(filelist[0], format='NETCDF4') - - #clip the lat/lon grid according to user input - #http://www.pyngl.ucar.edu/NioExtendedSelection.shtml - latsraw = tmp.variables[mergLatVarName][mergLatVarName+"|"+LATMIN+":"+LATMAX].astype('f2') - lonsraw = tmp.variables[mergLonVarName][mergLonVarName+"|"+LONMIN+":"+LONMAX].astype('f2') - lonsraw[lonsraw > 180] = lonsraw[lonsraw > 180] - 360. # convert to -180,180 if necessary + tmp = Dataset(filelist[0], 'r+',format='NETCDF4') + + alllatsraw = tmp.variables[mergLatVarName][:] + alllonsraw = tmp.variables[mergLonVarName][:] + alllonsraw[alllonsraw > 180] = alllonsraw[alllonsraw > 180] - 360. # convert to -180,180 if necessary + + #get the lat/lon info data (different resolution) + latminNETCDF = find_nearest(alllatsraw, float(LATMIN)) + latmaxNETCDF = find_nearest(alllatsraw, float(LATMAX)) + lonminNETCDF = find_nearest(alllonsraw, float(LONMIN)) + lonmaxNETCDF = find_nearest(alllonsraw, float(LONMAX)) + latminIndex = (np.where(alllatsraw == latminNETCDF))[0][0] + latmaxIndex = (np.where(alllatsraw == latmaxNETCDF))[0][0] + lonminIndex = (np.where(alllonsraw == lonminNETCDF))[0][0] + lonmaxIndex = (np.where(alllonsraw == lonmaxNETCDF))[0][0] + + #subsetting the data + latsraw = alllatsraw[latminIndex: latmaxIndex] + lonsraw = alllonsraw[lonminIndex:lonmaxIndex] LON, LAT = np.meshgrid(lonsraw, latsraw) #clean up @@ -157,13 +169,9 @@ def readMergData(dirname, filelist = None): for files in filelist: try: - #thisFile = Nio.open_file(files, format='nc') thisFile = Dataset(files, format='NETCDF4') #clip the dataset according to user lat,lon coordinates - #mask the data and fill with zeros for later - tempRaw = thisFile.variables[mergVarName][mergLatVarName+"|"+LATMIN+":"+LATMAX \ - +" " +mergLonVarName+"|"+LONMIN+":"+LONMAX ].astype('int16') - + tempRaw = thisFile.variables[mergVarName][:,latminIndex:latmaxIndex,lonminIndex:lonmaxIndex].astype('int16') tempMask = ma.masked_array(tempRaw, mask=(tempRaw > T_BB_MAX), fill_value=0) #get the actual values that the mask returned @@ -171,10 +179,10 @@ def readMergData(dirname, filelist = None): for index, value in maenumerate(tempMask): time_index, lat_index, lon_index = index tempMaskedValue[time_index,lat_index, lon_index]=value - - timesRaw = thisFile.variables[mergTimeVarName] + + xtimes = thisFile.variables[mergTimeVarName] #convert this time to a python datastring - time2store, _ = process.getModelTimes(files, mergTimeVarName) + time2store, _ = getModelTimes(xtimes, mergTimeVarName) #extend instead of append because getModelTimes returns a list already and we don't #want a list of list timelist.extend(time2store) @@ -183,7 +191,7 @@ def readMergData(dirname, filelist = None): thisFile = None except: - print "bad file! ", file + print "bad file! ", files mergImgs = ma.array(mergImgs) @@ -373,7 +381,8 @@ def findCloudElements(mergImgs,timelist,TRMMdirName=None): #---------regrid the TRMM data to the MERG dataset ---------------------------------- #regrid using the do_regrid stuff from the Apache OCW regriddedTRMM = ma.zeros((0, nygrd, nxgrd)) - regriddedTRMM = process.do_regrid(precipRateMasked[0,:,:], LATTRMM, LONTRMM, LAT, LON, order=1, mdi= -999999999) + #regriddedTRMM = process.do_regrid(precipRateMasked[0,:,:], LATTRMM, LONTRMM, LAT, LON, order=1, mdi= -999999999) + regriddedTRMM = do_regrid(precipRateMasked[0,:,:], LATTRMM, LONTRMM, LAT, LON, order=1, mdi= -999999999) #---------------------------------------------------------------------------------- # #get the lat/lon info from cloudElement @@ -592,7 +601,7 @@ def findCloudElements(mergImgs,timelist,TRMMdirName=None): #hierachial graph output graphTitle = "Cloud Elements observed over somewhere from 0000Z to 0000Z" - drawGraph(CLOUD_ELEMENT_GRAPH, graphTitle, edgeWeight) + #drawGraph(CLOUD_ELEMENT_GRAPH, graphTitle, edgeWeight) return CLOUD_ELEMENT_GRAPH #****************************************************************** @@ -665,7 +674,8 @@ def findPrecipRate(TRMMdirName, timelist): #---------regrid the TRMM data to the MERG dataset ---------------------------------- #regrid using the do_regrid stuff from the Apache OCW regriddedTRMM = ma.zeros((0, nygrd, nxgrd)) - regriddedTRMM = process.do_regrid(precipRateMasked[0,:,:], LATTRMM, LONTRMM, LAT, LON, order=1, mdi= -999999999) + regriddedTRMM = do_regrid(precipRateMasked[0,:,:], LATTRMM, LONTRMM, LAT, LON, order=1, mdi= -999999999) + #regriddedTRMM = process.do_regrid(precipRateMasked[0,:,:], LATTRMM, LONTRMM, LAT, LON, order=1, mdi= -999999999) #---------------------------------------------------------------------------------- # #get the lat/lon info from @@ -827,7 +837,7 @@ def findCloudClusters(CEGraph): print ("*"*80) graphTitle = "Cloud Clusters observed over somewhere during sometime" - drawGraph(PRUNED_GRAPH, graphTitle, edgeWeight) + #drawGraph(PRUNED_GRAPH, graphTitle, edgeWeight) cloudClustersFile.close return PRUNED_GRAPH @@ -894,7 +904,7 @@ def findMCC (prunedGraph): for node in prunedGraph.successors(eachNode): aSubGraph.add_edge(eachNode,node,weight=edgeWeight[0]) imgTitle = 'CC'+str(imgCount+1) - drawGraph(aSubGraph, imgTitle, edgeWeight) #for eachNode in path: + #drawGraph(aSubGraph, imgTitle, edgeWeight) #for eachNode in path: imgCount +=1 #----------end build back --------------------------------------------- @@ -2456,6 +2466,192 @@ def tempMaskedImages(imgFilename): subprocess.call(gradscmd, shell=True) return #****************************************************************** +def getModelTimes(xtimes, timeVarName): + ''' + Taken from process.py, removed the file opening at the beginning + TODO: Do a better job handling dates here + Routine to convert from model times ('hours since 1900...', 'days since ...') + into a python datetime structure + + Input:: + modelFile - path to the model tile you want to extract the times list and modelTimeStep from + timeVarName - name of the time variable in the model file + + Output:: + times - list of python datetime objects describing model data times + modelTimeStep - 'hourly','daily','monthly','annual' + ''' + + timeFormat = xtimes.units + # search to check if 'since' appears in units + try: + sinceLoc = re.search('since', timeFormat).end() + + except AttributeError: + print 'Error decoding model times: time variable attributes do not contain "since"' + raise + + units = None + TIME_UNITS = ('minutes', 'hours', 'days', 'months', 'years') + # search for 'seconds','minutes','hours', 'days', 'months', 'years' so know units + for unit in TIME_UNITS: + if re.search(unit, timeFormat): + units = unit + break + + # cut out base time (the bit following 'since') + base_time_string = string.lstrip(timeFormat[sinceLoc:]) + # decode base time + base_time = decodeTimeFromString(base_time_string) + + times = [] + + for xtime in xtimes[:]: + # Cast time as an int + #TODO: KDW this may cause problems for data that is hourly with more than one timestep in it + xtime = int(xtime) + + if int(xtime) == 0: + xtime = 1 + + if units == 'minutes': + dt = timedelta(minutes=xtime) + new_time = base_time + dt + elif units == 'hours': + dt = timedelta(hours=int(xtime)) + new_time = base_time + dt# timedelta(hours=int(xtime)) + elif units == 'days': + dt = timedelta(days=xtime) + new_time = base_time + dt + elif units == 'months': + # NB. adding months in python is complicated as month length varies and hence ambiguous. + # Perform date arithmatic manually + # Assumption: the base_date will usually be the first of the month + # NB. this method will fail if the base time is on the 29th or higher day of month + # -as can't have, e.g. Feb 31st. + new_month = int(base_time.month + xtime % 12) + new_year = int(math.floor(base_time.year + xtime / 12.)) + new_time = datetime.datetime(new_year, new_month, base_time.day, base_time.hour, base_time.second, 0) + elif units == 'years': + dt = datetime.timedelta(years=xtime) + new_time = base_time + dt + + times.append(new_time) + + try: + if len(xtimes) == 1: + timeStepLength = 0 + else: + timeStepLength = int(xtimes[1] - xtimes[0] + 1.e-12) + + modelTimeStep = getModelTimeStep(units, timeStepLength) + + #if timeStepLength is zero do not normalize times as this would create an empty list for MERG (hourly) data + if timeStepLength != 0: + times = normalizeDatetimes(times, modelTimeStep) + except: + raise + + return times, modelTimeStep +#****************************************************************** +def getModelTimeStep(units, stepSize): + # Time units are now determined. Determine the time intervals of input data (mdlTimeStep) + ''' + Taken from process.py + ''' + if units == 'minutes': + if stepSize == 60: + modelTimeStep = 'hourly' + elif stepSize == 1440: + modelTimeStep = 'daily' + # 28 days through 31 days + elif 40320 <= stepSize <= 44640: + modelTimeStep = 'monthly' + # 365 days through 366 days + elif 525600 <= stepSize <= 527040: + modelTimeStep = 'annual' + else: + raise Exception('model data time step interval exceeds the max time interval (annual)', units, stepSize) + + elif units == 'hours': + #need a check for fractional hrs and only one hr i.e. stepSize=0 e.g. with MERG data + if stepSize == 0 or stepSize == 1: + modelTimeStep = 'hourly' + elif stepSize == 24: + modelTimeStep = 'daily' + elif 672 <= stepSize <= 744: + modelTimeStep = 'monthly' + elif 8760 <= stepSize <= 8784: + modelTimeStep = 'annual' + else: + raise Exception('model data time step interval exceeds the max time interval (annual)', units, stepSize) + + elif units == 'days': + if stepSize == 1: + modelTimeStep = 'daily' + elif 28 <= stepSize <= 31: + modelTimeStep = 'monthly' + elif 365 <= stepSize <= 366: + modelTimeStep = 'annual' + else: + raise Exception('model data time step interval exceeds the max time interval (annual)', units, stepSize) + + elif units == 'months': + if stepSize == 1: + modelTimeStep = 'monthly' + elif stepSize == 12: + modelTimeStep = 'annual' + else: + raise Exception('model data time step interval exceeds the max time interval (annual)', units, stepSize) + + elif units == 'years': + if stepSize == 1: + modelTimeStep = 'annual' + else: + raise Exception('model data time step interval exceeds the max time interval (annual)', units, stepSize) + + else: + errorMessage = 'the time unit ', units, ' is not currently handled in this version.' + raise Exception(errorMessage) + + return modelTimeStep +#****************************************************************** +def decodeTimeFromString(time_string): + ''' + Taken from process.py + Decodes string into a python datetime object + *Method:* tries a bunch of different time format possibilities and hopefully one of them will hit. + :: + + **Input:** time_string - a string that represents a date/time + + **Output:** mytime - a python datetime object + ''' + # This will deal with times that use decimal seconds + if '.' in time_string: + time_string = time_string.split('.')[0] + '0' + else: + pass + + try: + mytime = datetime.strptime(time_string,'%Y-%m-%d %H') + return mytime + + except ValueError: + pass + + print 'Error decoding time string: string does not match a predefined time format' + return 0 +#****************************************************************** +def do_regrid(q, lat, lon, lat2, lon2, order=1, mdi=-999999999): + """ + This function has been moved to the ocw/dataset_processor module + """ + from ocw import dataset_processor + q2 = dataset_processor._rcmes_spatial_regrid(q, lat, lon, lat2, lon2, order=1) + + return q2 +#****************************************************************** # # METRICS FUNCTIONS FOR MERG.PY # TODO: rewrite these metrics so that they use the data from the http://git-wip-us.apache.org/repos/asf/climate/blob/48352bc0/mccsearch/code/mccSearchUI.py ---------------------------------------------------------------------- diff --git a/mccsearch/code/mccSearchUI.py b/mccsearch/code/mccSearchUI.py index 885ff43..91670e4 100644 --- a/mccsearch/code/mccSearchUI.py +++ b/mccsearch/code/mccSearchUI.py @@ -28,7 +28,7 @@ import subprocess #mccSearch modules import mccSearch -import files +#import files def main(): CEGraph = nx.DiGraph()
