climate-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From jo...@apache.org
Subject [1/2] git commit: removed process.py dependency and fixed Nio toNETCDF logical errors
Date Mon, 27 Oct 2014 19:21:30 GMT
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 <k_whitehall@yahoo.com>
Authored: Fri Oct 24 15:23:29 2014 -0700
Committer: Kim Whitehall <k_whitehall@yahoo.com>
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()


Mime
View raw message