http://git-wip-us.apache.org/repos/asf/climate/blob/9bb21700/mccsearch/code/mccSearch.py
----------------------------------------------------------------------
diff --git a/mccsearch/code/mccSearch.py b/mccsearch/code/mccSearch.py
index 15ce0c6..6f1fa26 100644
--- a/mccsearch/code/mccSearch.py
+++ b/mccsearch/code/mccSearch.py
@@ -56,7 +56,7 @@ TRES = 1                              #temporal resolution in 
hrs
 LAT_DISTANCE = 111.0   #the avg distance in km for 1deg lat for the region 
being considered 
 LON_DISTANCE = 111.0   #the avg distance in km for 1deg lon for the region 
being considered
 STRUCTURING_ELEMENT = [[0,1,0],[1,1,1],[0,1,0]] #the matrix for determining 
the pattern for the contiguous boxes and must
-                                                                               
                #have same rank of the matrix it is being compared against 
+                                                #have same rank of the matrix 
it is being compared against 
 #criteria for determining cloud elements and edges
 T_BB_MAX = 243  #warmest temp to allow (-30C to -55C according to Morel and 
Sensi 2002)
 T_BB_MIN = 218  #cooler temp for the center of the system
@@ -84,2079 +84,2079 @@ PRUNED_GRAPH = nx.DiGraph()
 #************************ Begin Functions *************************
 #******************************************************************
 def readMergData(dirname, filelist = None):
-       '''
-       Purpose::
-               Read MERG data into RCMES format
-       
-       Input::
-               dirname: a string representing the directory to the MERG files 
in NETCDF format
-               filelist (optional): a list of strings representing the 
filenames betweent the start and end dates provided
-       
-       Output::
-               A 3D masked array (t,lat,lon) with only the variables which 
meet the minimum temperature 
-               criteria for each frame
-
-       Assumptions::
-               The MERG data has been converted to NETCDF using LATS4D
-               The data has the same lat/lon format
-
-       TODO:: figure out how to use netCDF4 to do the clipping tmp = 
netCDF4.Dataset(filelist[0])
-
-       '''
-
-       global LAT
-       global LON
-
-       # these strings are specific to the MERG data
-       mergVarName = 'ch4'
-       mergTimeVarName = 'time'
-       mergLatVarName = 'latitude'
-       mergLonVarName = 'longitude'
-       
-       filelistInstructions = dirname + '/*'
-       if filelist == None:
-               filelist = glob.glob(filelistInstructions)
-
-       
-       #sat_img is the array that will contain all the masked frames
-       mergImgs = []
-       #timelist of python time strings
-       timelist = [] 
-       time2store = None
-       tempMaskedValueNp =[]
-       
-
-       filelist.sort()
-       nfiles = len(filelist)
-
-       # Crash nicely if there are no netcdf files
-       if nfiles == 0:
-               print 'Error: no files in this directory! Exiting elegantly'
-               sys.exit()
-       else:
-               # Open the first file in the list to read in lats, lons and 
generate the  grid for comparison
-               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
-               latsraw =[]
-               lonsraw = []
-               nygrd = len(LAT[:, 0]); nxgrd = len(LON[0, :])
-               tmp.close
-       
-       for files in filelist:
-               try:
-                       thisFile = Dataset(files, format='NETCDF4')
-                       #clip the dataset according to user lat,lon coordinates
-                       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
-                       tempMaskedValue = 
ma.zeros((tempRaw.shape)).astype('int16')
-                       for index, value in maenumerate(tempMask): 
-                               time_index, lat_index, lon_index = index        
                
-                               tempMaskedValue[time_index,lat_index, 
lon_index]=value  
-                       
-                       xtimes = thisFile.variables[mergTimeVarName]
-                       #convert this time to a python datastring
-                       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)
-                       mergImgs.extend(tempMaskedValue) 
-                       thisFile.close
-                       thisFile = None
-                       
-               except:
-                       print "bad file! ", files
-
-       mergImgs = ma.array(mergImgs)
-
-       return mergImgs, timelist
+    '''
+    Purpose::
+        Read MERG data into RCMES format
+    
+    Input::
+        dirname: a string representing the directory to the MERG files in 
NETCDF format
+        filelist (optional): a list of strings representing the filenames 
betweent the start and end dates provided
+    
+    Output::
+        A 3D masked array (t,lat,lon) with only the variables which meet the 
minimum temperature 
+        criteria for each frame
+
+    Assumptions::
+        The MERG data has been converted to NETCDF using LATS4D
+        The data has the same lat/lon format
+
+    TODO:: figure out how to use netCDF4 to do the clipping tmp = 
netCDF4.Dataset(filelist[0])
+
+    '''
+
+    global LAT
+    global LON
+
+    # these strings are specific to the MERG data
+    mergVarName = 'ch4'
+    mergTimeVarName = 'time'
+    mergLatVarName = 'latitude'
+    mergLonVarName = 'longitude'
+    
+    filelistInstructions = dirname + '/*'
+    if filelist == None:
+        filelist = glob.glob(filelistInstructions)
+
+    
+    #sat_img is the array that will contain all the masked frames
+    mergImgs = []
+    #timelist of python time strings
+    timelist = [] 
+    time2store = None
+    tempMaskedValueNp =[]
+    
+
+    filelist.sort()
+    nfiles = len(filelist)
+
+    # Crash nicely if there are no netcdf files
+    if nfiles == 0:
+        print 'Error: no files in this directory! Exiting elegantly'
+        sys.exit()
+    else:
+        # Open the first file in the list to read in lats, lons and generate 
the  grid for comparison
+        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
+        latsraw =[]
+        lonsraw = []
+        nygrd = len(LAT[:, 0]); nxgrd = len(LON[0, :])
+        tmp.close
+    
+    for files in filelist:
+        try:
+            thisFile = Dataset(files, format='NETCDF4')
+            #clip the dataset according to user lat,lon coordinates
+            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
+            tempMaskedValue = ma.zeros((tempRaw.shape)).astype('int16')
+            for index, value in maenumerate(tempMask): 
+                time_index, lat_index, lon_index = index                       
+                tempMaskedValue[time_index,lat_index, lon_index]=value 
+            
+            xtimes = thisFile.variables[mergTimeVarName]
+            #convert this time to a python datastring
+            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)
+            mergImgs.extend(tempMaskedValue) 
+            thisFile.close
+            thisFile = None
+            
+        except:
+            print "bad file! ", files
+
+    mergImgs = ma.array(mergImgs)
+
+    return mergImgs, timelist
 #******************************************************************
 def findCloudElements(mergImgs,timelist,TRMMdirName=None):
-       '''
-       Purpose::
-               Determines the contiguous boxes for a given time of the 
satellite images i.e. each frame
-               using scipy ndimage package
-       
-       Input:: 
-               mergImgs: masked numpy array in (time,lat,lon),T_bb 
representing the satellite data. This is masked based on the
-               maximum acceptable temperature, T_BB_MAX
-               timelist: a list of python datatimes
-               TRMMdirName (optional): string representing the path where to 
find the TRMM datafiles
-               
-       Output::
-               CLOUD_ELEMENT_GRAPH: a Networkx directed graph where each node 
contains the information in cloudElementDict
-               The nodes are determined according to the area of contiguous 
squares. The nodes are linked through weighted edges.
-
-               cloudElementDict = {'uniqueID': unique tag for this CE, 
-                                                       'cloudElementTime': 
time of the CE,
-                                                       'cloudElementLatLon': 
(lat,lon,value) of MERG data of CE, 
-                                                       
'cloudElementCenter':list of floating-point [lat,lon] representing the CE's 
center 
-                                                       
'cloudElementArea':floating-point representing the area of the CE, 
-                                                       
'cloudElementEccentricity': floating-point representing the shape of the CE, 
-                                                       
'cloudElementTmax':integer representing the maximum Tb in CE, 
-                                                       'cloudElementTmin': 
integer representing the minimum Tb in CE, 
-                                                       
'cloudElementPrecipTotal':floating-point representing the sum of all rainfall 
in CE if TRMMdirName entered,
-                                                       
'cloudElementLatLonTRMM':(lat,lon,value) of TRMM data in CE if TRMMdirName 
entered, 
-                                                       'TRMMArea': 
floating-point representing the CE if TRMMdirName entered,
-                                                       
'CETRMMmax':floating-point representing the max rate in the CE if TRMMdirName 
entered, 
-                                                       
'CETRMMmin':floating-point representing the min rate in the CE if TRMMdirName 
entered}
-       Assumptions::
-               Assumes we are dealing with MERG data which is 4kmx4km 
resolved, thus the smallest value 
-               required according to Vila et al. (2008) is 2400km^2 
-               therefore, 2400/16 = 150 contiguous squares
-       '''
-
-       frame = ma.empty((1,mergImgs.shape[1],mergImgs.shape[2]))
-       CEcounter = 0
-       frameCEcounter = 0
-       frameNum = 0
-       cloudElementEpsilon = 0.0
-       cloudElementDict = {} 
-       cloudElementCenter = []         #list with two elements [lat,lon] for 
the center of a CE
-       prevFrameCEs = []                       #list for CEs in previous frame
-       currFrameCEs = []                       #list for CEs in current frame
-       cloudElementLat = []            #list for a particular CE's lat values
-       cloudElementLon = []            #list for a particular CE's lon values
-       cloudElementLatLons = []        #list for a particular CE's (lat,lon) 
values
-       
-       prevLatValue = 0.0
-       prevLonValue = 0.0
-       TIR_min = 0.0
-       TIR_max = 0.0
-       temporalRes = 3 # TRMM data is 3 hourly
-       precipTotal = 0.0
-       CETRMMList =[]
-       precip =[]
-       TRMMCloudElementLatLons =[]
-
-       minCELatLimit = 0.0
-       minCELonLimit = 0.0
-       maxCELatLimit = 0.0
-       maxCELonLimit = 0.0
-       
-       nygrd = len(LAT[:, 0]); nxgrd = len(LON[0, :])
-       
-       #openfile for storing ALL cloudElement information 
-       cloudElementsFile = 
open((MAINDIRECTORY+'/textFiles/cloudElements.txt'),'wb')
-       #openfile for storing cloudElement information meeting user criteria 
i.e. MCCs in this case
-       cloudElementsUserFile = 
open((MAINDIRECTORY+'/textFiles/cloudElementsUserFile.txt'),'w')
-       
-       #NB in the TRMM files the info is hours since the time thus 00Z file 
has in 01, 02 and 03 times
-       for t in xrange(mergImgs.shape[0]):
-               #-------------------------------------------------
-               # #textfile name for saving the data for arcgis
-               # thisFileName = MAINDIRECTORY+'/' + 
(str(timelist[t])).replace(" ", "_") + '.txt'
-               # cloudElementsTextFile = open(thisFileName,'w')
-               #-------------------------------------------------
-
-               #determine contiguous locations with temeperature below the 
warmest temp i.e. cloudElements in each frame
-               frame, CEcounter = ndimage.measurements.label(mergImgs[t,:,:], 
structure=STRUCTURING_ELEMENT)
-               frameCEcounter=0
-               frameNum += 1
-
-               #for each of the areas identified, check to determine if it a 
valid CE via an area and T requirement
-               for count in xrange(CEcounter):
-                       #[0] is time dimension. Determine the actual values 
from the data
-                       #loc is a masked array
-                       try:
-                               loc = ndimage.find_objects(frame==(count+1))[0]
-                       except Exception, e:
-                               print "Error is ", e
-                               continue
-
-
-                       cloudElement = mergImgs[t,:,:][loc]
-                       labels, lcounter = ndimage.label(cloudElement)
-                       
-                       #determine the true lats and lons for this particular CE
-                       cloudElementLat = LAT[loc[0],0]
-                       cloudElementLon = LON[0,loc[1]] 
-                       
-                       #determine number of boxes in this cloudelement
-                       numOfBoxes = np.count_nonzero(cloudElement)
-                       cloudElementArea = numOfBoxes*XRES*YRES
-
-                       #If the area is greater than the area required, or if 
the area is smaller than the suggested area, check if it meets a convective 
fraction requirement
-                       #consider as CE
-
-                       if cloudElementArea >= AREA_MIN or (cloudElementArea < 
AREA_MIN and ((ndimage.minimum(cloudElement, 
labels=labels))/float((ndimage.maximum(cloudElement, labels=labels)))) < 
CONVECTIVE_FRACTION ):
-
-                               #get some time information and labeling info
-                               frameTime = str(timelist[t])
-                               frameCEcounter +=1
-                               CEuniqueID = 
'F'+str(frameNum)+'CE'+str(frameCEcounter) 
-
-                               
#-------------------------------------------------
-                               #textfile name for accesing CE data using 
MATLAB code
-                               # thisFileName = MAINDIRECTORY+'/' + 
(str(timelist[t])).replace(" ", "_") + CEuniqueID +'.txt'
-                               # cloudElementsTextFile = open(thisFileName,'w')
-                               
#-------------------------------------------------
-
-                               # ------ NETCDF File stuff for brightness temp 
stuff ------------------------------------
-                               thisFileName = MAINDIRECTORY 
+'/MERGnetcdfCEs/cloudElements'+ (str(timelist[t])).replace(" ", "_") + 
CEuniqueID +'.nc'
-                               currNetCDFCEData = Dataset(thisFileName, 'w', 
format='NETCDF4')
-                               currNetCDFCEData.description = 'Cloud Element 
'+CEuniqueID + ' temperature data'
-                               currNetCDFCEData.calendar = 'standard'
-                               currNetCDFCEData.conventions = 'COARDS'
-                               # dimensions
-                               currNetCDFCEData.createDimension('time', None)
-                               currNetCDFCEData.createDimension('lat', 
len(LAT[:,0]))
-                               currNetCDFCEData.createDimension('lon', 
len(LON[0,:]))
-                               # variables
-                               tempDims = ('time','lat', 'lon',)
-                               times = currNetCDFCEData.createVariable('time', 
'f8', ('time',))
-                               times.units = 'hours since '+ 
str(timelist[t])[:-6]
-                               latitudes = 
currNetCDFCEData.createVariable('latitude', 'f8', ('lat',))
-                               longitudes = 
currNetCDFCEData.createVariable('longitude', 'f8', ('lon',))
-                               brightnesstemp = 
currNetCDFCEData.createVariable('brightnesstemp', 'i16',tempDims )
-                               brightnesstemp.units = 'Kelvin'
-                               # NETCDF data
-                               dates=[timelist[t]+timedelta(hours=0)]
-                               times[:] =  date2num(dates,units=times.units)
-                               longitudes[:] = LON[0,:]
-                               longitudes.units = "degrees_east" 
-                               longitudes.long_name = "Longitude" 
-
-                               latitudes[:] =  LAT[:,0]
-                               latitudes.units = "degrees_north"
-                               latitudes.long_name ="Latitude"
-                               
-                               #generate array of zeros for brightness 
temperature
-                               brightnesstemp1 = ma.zeros((1,len(latitudes), 
len(longitudes))).astype('int16')
-                               #-----------End most of NETCDF file stuff 
------------------------------------
-
-                               #if other dataset (TRMM) assumed to be a 
precipitation dataset was entered
-                               if TRMMdirName:
-                                       #------------------TRMM stuff 
-------------------------------------------------
-                                       fileDate = 
((str(timelist[t])).replace(" ", "")[:-8]).replace("-","")
-                                       fileHr1 = (str(timelist[t])).replace(" 
", "")[-8:-6]
-                                       
-                                       if int(fileHr1) % temporalRes == 0:
-                                               fileHr = fileHr1
-                                       else:
-                                               fileHr = 
(int(fileHr1)/temporalRes) * temporalRes
-                                       if fileHr < 10:
-                                               fileHr = '0'+str(fileHr)
-                                       else:
-                                               str(fileHr)
-
-                                       #open TRMM file for the resolution info 
and to create the appropriate sized grid
-                                       TRMMfileName = TRMMdirName+'/3B42.'+ 
fileDate + "."+str(fileHr)+".7A.nc"
-                                       
-                                       TRMMData = Dataset(TRMMfileName,'r', 
format='NETCDF4')
-                                       precipRate = 
TRMMData.variables['pcp'][:,:,:]
-                                       latsrawTRMMData = 
TRMMData.variables['latitude'][:]
-                                       lonsrawTRMMData = 
TRMMData.variables['longitude'][:]
-                                       lonsrawTRMMData[lonsrawTRMMData > 180] 
= lonsrawTRMMData[lonsrawTRMMData>180] - 360.
-                                       LONTRMM, LATTRMM = 
np.meshgrid(lonsrawTRMMData, latsrawTRMMData)
-
-                                       nygrdTRMM = len(LATTRMM[:,0]); 
nxgrdTRMM = len(LONTRMM[0,:])
-                                       precipRateMasked = 
ma.masked_array(precipRate, mask=(precipRate < 0.0))
-                                       #---------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)
-                                       
#----------------------------------------------------------------------------------
-               
-                                       # #get the lat/lon info from 
cloudElement
-                                       #get the lat/lon info from the file
-                                       latCEStart = LAT[0][0]
-                                       latCEEnd = LAT[-1][0]
-                                       lonCEStart = LON[0][0]
-                                       lonCEEnd = LON[0][-1]
-                                       
-                                       #get the lat/lon info for TRMM data 
(different resolution)
-                                       latStartT = 
find_nearest(latsrawTRMMData, latCEStart)
-                                       latEndT = find_nearest(latsrawTRMMData, 
latCEEnd)
-                                       lonStartT = 
find_nearest(lonsrawTRMMData, lonCEStart)
-                                       lonEndT = find_nearest(lonsrawTRMMData, 
lonCEEnd)
-                                       latStartIndex = 
np.where(latsrawTRMMData == latStartT)
-                                       latEndIndex = np.where(latsrawTRMMData 
== latEndT)
-                                       lonStartIndex = 
np.where(lonsrawTRMMData == lonStartT)
-                                       lonEndIndex = np.where(lonsrawTRMMData 
== lonEndT)
-
-                                       #get the relevant TRMM info 
-                                       CEprecipRate = 
precipRate[:,(latStartIndex[0][0]-1):latEndIndex[0][0],(lonStartIndex[0][0]-1):lonEndIndex[0][0]]
-                                       TRMMData.close()
-                                       
-                                       # ------ NETCDF File info for writing 
TRMM CE rainfall ------------------------------------
-                                       thisFileName = 
MAINDIRECTORY+'/TRMMnetcdfCEs/TRMM' + (str(timelist[t])).replace(" ", "_") + 
CEuniqueID +'.nc'
-                                       currNetCDFTRMMData = 
Dataset(thisFileName, 'w', format='NETCDF4')
-                                       currNetCDFTRMMData.description = 'Cloud 
Element '+CEuniqueID + ' precipitation data'
-                                       currNetCDFTRMMData.calendar = 'standard'
-                                       currNetCDFTRMMData.conventions = 
'COARDS'
-                                       # dimensions
-                                       
currNetCDFTRMMData.createDimension('time', None)
-                                       
currNetCDFTRMMData.createDimension('lat', len(LAT[:,0]))
-                                       
currNetCDFTRMMData.createDimension('lon', len(LON[0,:]))
-                                       
-                                       # variables
-                                       TRMMprecip = ('time','lat', 'lon',)
-                                       times = 
currNetCDFTRMMData.createVariable('time', 'f8', ('time',))
-                                       times.units = 'hours since '+ 
str(timelist[t])[:-6]
-                                       latitude = 
currNetCDFTRMMData.createVariable('latitude', 'f8', ('lat',))
-                                       longitude = 
currNetCDFTRMMData.createVariable('longitude', 'f8', ('lon',))
-                                       rainFallacc = 
currNetCDFTRMMData.createVariable('precipitation_Accumulation', 'f8',TRMMprecip 
)
-                                       rainFallacc.units = 'mm'
-
-                                       longitude[:] = LON[0,:]
-                                       longitude.units = "degrees_east" 
-                                       longitude.long_name = "Longitude" 
-
-                                       latitude[:] =  LAT[:,0]
-                                       latitude.units = "degrees_north"
-                                       latitude.long_name ="Latitude"
-
-                                       finalCETRMMvalues = 
ma.zeros((brightnesstemp.shape))
-                                       #-----------End most of NETCDF file 
stuff ------------------------------------
-
-                               #populate cloudElementLatLons by unpacking the 
original values from loc to get the actual value for lat and lon
-                               #TODO: KDW - too dirty... play with 
itertools.izip or zip and the enumerate with this
-                               #                       as cloudElement is 
masked
-                               for index,value in np.ndenumerate(cloudElement):
-                                       if value != 0 : 
-                                               lat_index,lon_index = index
-                                               lat_lon_tuple = 
(cloudElementLat[lat_index], cloudElementLon[lon_index],value)
-
-                                               #generate the comma separated 
file for GIS
-                                               
cloudElementLatLons.append(lat_lon_tuple)
-
-                                               #temp data for CE NETCDF file
-                                               
brightnesstemp1[0,int(np.where(LAT[:,0]==cloudElementLat[lat_index])[0]),int(np.where(LON[0,:]==cloudElementLon[lon_index])[0])]
 = value
-                                               
-                                               if TRMMdirName:
-                                                       
finalCETRMMvalues[0,int(np.where(LAT[:,0]==cloudElementLat[lat_index])[0]),int(np.where(LON[0,:]==cloudElementLon[lon_index])[0])]
 = 
regriddedTRMM[int(np.where(LAT[:,0]==cloudElementLat[lat_index])[0]),int(np.where(LON[0,:]==cloudElementLon[lon_index])[0])]
-                                                       
CETRMMList.append((cloudElementLat[lat_index], cloudElementLon[lon_index], 
finalCETRMMvalues[0,cloudElementLat[lat_index], cloudElementLon[lon_index]]))
-
-
-                               brightnesstemp[:] = brightnesstemp1[:]
-                               currNetCDFCEData.close()
-
-                               if TRMMdirName:
-
-                                       #calculate the total precip associated 
with the feature
-                                       for index, value in 
np.ndenumerate(finalCETRMMvalues):
-                                               precipTotal += value 
-                                               precip.append(value)
-                       
-                                       rainFallacc[:] = finalCETRMMvalues[:]
-                                       currNetCDFTRMMData.close()
-                                       TRMMnumOfBoxes = 
np.count_nonzero(finalCETRMMvalues)
-                                       TRMMArea = TRMMnumOfBoxes*XRES*YRES
-                                       try:
-                                               maxCEprecipRate = 
np.max(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
-                                               minCEprecipRate = 
np.min(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
-                                       except:
-                                               pass
-
-                               #sort cloudElementLatLons by lats
-                               cloudElementLatLons.sort(key=lambda tup: 
tup[0])        
-
-                               #determine if the cloud element the shape 
-                               cloudElementEpsilon = eccentricity 
(cloudElement)
-                               cloudElementsUserFile.write("\n\nTime is: %s" 
%(str(timelist[t])))
-                               cloudElementsUserFile.write("\nCEuniqueID is: 
%s" %CEuniqueID)
-                               latCenter, lonCenter = 
ndimage.measurements.center_of_mass(cloudElement, labels=labels)
-                               
-                               #latCenter and lonCenter are given according to 
the particular array defining this CE
-                               #so you need to convert this value to the 
overall domain truth
-                               latCenter = cloudElementLat[round(latCenter)]
-                               lonCenter = cloudElementLon[round(lonCenter)]
-                               cloudElementsUserFile.write("\nCenter (lat,lon) 
is: %.2f\t%.2f" %(latCenter, lonCenter))
-                               cloudElementCenter.append(latCenter)
-                               cloudElementCenter.append(lonCenter)
-                               cloudElementsUserFile.write("\nNumber of boxes 
are: %d" %numOfBoxes)
-                               cloudElementsUserFile.write("\nArea is: %.4f 
km^2" %(cloudElementArea))
-                               cloudElementsUserFile.write("\nAverage 
brightness temperature is: %.4f K" %ndimage.mean(cloudElement, labels=labels))
-                               cloudElementsUserFile.write("\nMin brightness 
temperature is: %.4f K" %ndimage.minimum(cloudElement, labels=labels))
-                               cloudElementsUserFile.write("\nMax brightness 
temperature is: %.4f K" %ndimage.maximum(cloudElement, labels=labels))
-                               cloudElementsUserFile.write("\nBrightness 
temperature variance is: %.4f K" %ndimage.variance(cloudElement, labels=labels))
-                               cloudElementsUserFile.write("\nConvective 
fraction is: %.4f " %(((ndimage.minimum(cloudElement, 
labels=labels))/float((ndimage.maximum(cloudElement, labels=labels))))*100.0))
-                               cloudElementsUserFile.write("\nEccentricity is: 
%.4f " %(cloudElementEpsilon))
-                               #populate the dictionary
-                               if TRMMdirName:
-                                       cloudElementDict = {'uniqueID': 
CEuniqueID, 'cloudElementTime': timelist[t],'cloudElementLatLon': 
cloudElementLatLons, 'cloudElementCenter':cloudElementCenter, 
'cloudElementArea':cloudElementArea, 
'cloudElementEccentricity':cloudElementEpsilon, 'cloudElementTmax':TIR_max, 
'cloudElementTmin': TIR_min, 
'cloudElementPrecipTotal':precipTotal,'cloudElementLatLonTRMM':CETRMMList, 
'TRMMArea': TRMMArea,'CETRMMmax':maxCEprecipRate, 'CETRMMmin':minCEprecipRate}
-                               else:
-                                       cloudElementDict = {'uniqueID': 
CEuniqueID, 'cloudElementTime': timelist[t],'cloudElementLatLon': 
cloudElementLatLons, 'cloudElementCenter':cloudElementCenter, 
'cloudElementArea':cloudElementArea, 
'cloudElementEccentricity':cloudElementEpsilon, 'cloudElementTmax':TIR_max, 
'cloudElementTmin': TIR_min,}
-                               
-                               #current frame list of CEs
-                               currFrameCEs.append(cloudElementDict)
-                               
-                               #draw the graph node
-                               CLOUD_ELEMENT_GRAPH.add_node(CEuniqueID, 
cloudElementDict)
-                               
-                               if frameNum != 1:
-                                       for cloudElementDict in prevFrameCEs:
-                                               thisCElen = 
len(cloudElementLatLons)
-                                               percentageOverlap, areaOverlap 
= cloudElementOverlap(cloudElementLatLons, 
cloudElementDict['cloudElementLatLon'])
-                                               
-                                               #change weights to integers 
because the built in shortest path chokes on floating pts according to Networkx 
doc
-                                               #according to Goyens et al, two 
CEs are considered related if there is atleast 95% overlap between them for 
consecutive imgs a max of 2 hrs apart
-                                               if percentageOverlap >= 0.95: 
-                                                       
CLOUD_ELEMENT_GRAPH.add_edge(cloudElementDict['uniqueID'], CEuniqueID, 
weight=edgeWeight[0])
-                                                       
-                                               elif percentageOverlap >= 0.90 
and percentageOverlap < 0.95 :
-                                                       
CLOUD_ELEMENT_GRAPH.add_edge(cloudElementDict['uniqueID'], CEuniqueID, 
weight=edgeWeight[1])
-
-                                               elif areaOverlap >= MIN_OVERLAP:
-                                                       
CLOUD_ELEMENT_GRAPH.add_edge(cloudElementDict['uniqueID'], CEuniqueID, 
weight=edgeWeight[2])
-
-                               else:
-                                       #TODO: remove this else as we only wish 
for the CE details
-                                       #ensure only the non-zero elements are 
considered
-                                       #store intel in allCE file
-                                       labels, _ = ndimage.label(cloudElement)
-                                       
cloudElementsFile.write("\n-----------------------------------------------")
-                                       cloudElementsFile.write("\n\nTime is: 
%s" %(str(timelist[t])))
-                                       # cloudElementLat = LAT[loc[0],0]
-                                       # cloudElementLon = LON[0,loc[1]] 
-                                       
-                                       #populate cloudElementLatLons by 
unpacking the original values from loc
-                                       #TODO: KDW - too dirty... play with 
itertools.izip or zip and the enumerate with this
-                                       #                       as cloudElement 
is masked
-                                       for index,value in 
np.ndenumerate(cloudElement):
-                                               if value != 0 : 
-                                                       lat_index,lon_index = 
index
-                                                       lat_lon_tuple = 
(cloudElementLat[lat_index], cloudElementLon[lon_index])
-                                                       
cloudElementLatLons.append(lat_lon_tuple)
-       
-                                       cloudElementsFile.write("\nLocation of 
rejected CE (lat,lon) points are: %s" %cloudElementLatLons)
-                                       #latCenter and lonCenter are given 
according to the particular array defining this CE
-                                       #so you need to convert this value to 
the overall domain truth
-                                       latCenter, lonCenter = 
ndimage.measurements.center_of_mass(cloudElement, labels=labels)
-                                       latCenter = 
cloudElementLat[round(latCenter)]
-                                       lonCenter = 
cloudElementLon[round(lonCenter)]
-                                       cloudElementsFile.write("\nCenter 
(lat,lon) is: %.2f\t%.2f" %(latCenter, lonCenter))
-                                       cloudElementsFile.write("\nNumber of 
boxes are: %d" %numOfBoxes)
-                                       cloudElementsFile.write("\nArea is: 
%.4f km^2" %(cloudElementArea))
-                                       cloudElementsFile.write("\nAverage 
brightness temperature is: %.4f K" %ndimage.mean(cloudElement, labels=labels))
-                                       cloudElementsFile.write("\nMin 
brightness temperature is: %.4f K" %ndimage.minimum(cloudElement, 
labels=labels))
-                                       cloudElementsFile.write("\nMax 
brightness temperature is: %.4f K" %ndimage.maximum(cloudElement, 
labels=labels))
-                                       cloudElementsFile.write("\nBrightness 
temperature variance is: %.4f K" %ndimage.variance(cloudElement, labels=labels))
-                                       cloudElementsFile.write("\nConvective 
fraction is: %.4f " %(((ndimage.minimum(cloudElement, 
labels=labels))/float((ndimage.maximum(cloudElement, labels=labels))))*100.0))
-                                       cloudElementsFile.write("\nEccentricity 
is: %.4f " %(cloudElementEpsilon))
-                                       
cloudElementsFile.write("\n-----------------------------------------------")
-                                       
-                       #reset list for the next CE
-                       nodeExist = False
-                       cloudElementCenter=[]
-                       cloudElement = []
-                       cloudElementLat=[]
-                       cloudElementLon =[]
-                       cloudElementLatLons =[]
-                       brightnesstemp1 =[]
-                       brightnesstemp =[]
-                       finalCETRMMvalues =[]
-                       CEprecipRate =[]
-                       CETRMMList =[]
-                       precipTotal = 0.0
-                       precip=[]
-                       TRMMCloudElementLatLons=[]
-                       
-               #reset for the next time
-               prevFrameCEs =[]
-               prevFrameCEs = currFrameCEs
-               currFrameCEs =[]
-                                               
-       cloudElementsFile.close
-       cloudElementsUserFile.close
-       #if using ARCGIS data store code, uncomment this file close line
-       #cloudElementsTextFile.close
-
-       #clean up graph - remove parent and childless nodes
-       outAndInDeg = CLOUD_ELEMENT_GRAPH.degree_iter()
-       toRemove = [node[0] for node in outAndInDeg if node[1]<1]
-       CLOUD_ELEMENT_GRAPH.remove_nodes_from(toRemove)
-       
-       print "number of nodes are: ", CLOUD_ELEMENT_GRAPH.number_of_nodes()
-       print "number of edges are: ", CLOUD_ELEMENT_GRAPH.number_of_edges()
-       print ("*"*80)
-
-       #hierachial graph output
-       graphTitle = "Cloud Elements observed over somewhere from 0000Z to 
0000Z" 
-       #drawGraph(CLOUD_ELEMENT_GRAPH, graphTitle, edgeWeight)
-
-       return CLOUD_ELEMENT_GRAPH      
+    '''
+    Purpose::
+        Determines the contiguous boxes for a given time of the satellite 
images i.e. each frame
+        using scipy ndimage package
+    
+    Input::    
+        mergImgs: masked numpy array in (time,lat,lon),T_bb representing the 
satellite data. This is masked based on the
+        maximum acceptable temperature, T_BB_MAX
+        timelist: a list of python datatimes
+        TRMMdirName (optional): string representing the path where to find the 
TRMM datafiles
+        
+    Output::
+        CLOUD_ELEMENT_GRAPH: a Networkx directed graph where each node 
contains the information in cloudElementDict
+        The nodes are determined according to the area of contiguous squares. 
The nodes are linked through weighted edges.
+
+        cloudElementDict = {'uniqueID': unique tag for this CE, 
+                            'cloudElementTime': time of the CE,
+                            'cloudElementLatLon': (lat,lon,value) of MERG data 
of CE, 
+                            'cloudElementCenter':list of floating-point 
[lat,lon] representing the CE's center 
+                            'cloudElementArea':floating-point representing the 
area of the CE, 
+                            'cloudElementEccentricity': floating-point 
representing the shape of the CE, 
+                            'cloudElementTmax':integer representing the 
maximum Tb in CE, 
+                            'cloudElementTmin': integer representing the 
minimum Tb in CE, 
+                            'cloudElementPrecipTotal':floating-point 
representing the sum of all rainfall in CE if TRMMdirName entered,
+                            'cloudElementLatLonTRMM':(lat,lon,value) of TRMM 
data in CE if TRMMdirName entered, 
+                            'TRMMArea': floating-point representing the CE if 
TRMMdirName entered,
+                            'CETRMMmax':floating-point representing the max 
rate in the CE if TRMMdirName entered, 
+                            'CETRMMmin':floating-point representing the min 
rate in the CE if TRMMdirName entered}
+    Assumptions::
+        Assumes we are dealing with MERG data which is 4kmx4km resolved, thus 
the smallest value 
+        required according to Vila et al. (2008) is 2400km^2 
+        therefore, 2400/16 = 150 contiguous squares
+    '''
+
+    frame = ma.empty((1,mergImgs.shape[1],mergImgs.shape[2]))
+    CEcounter = 0
+    frameCEcounter = 0
+    frameNum = 0
+    cloudElementEpsilon = 0.0
+    cloudElementDict = {} 
+    cloudElementCenter = []            #list with two elements [lat,lon] for 
the center of a CE
+    prevFrameCEs = []                  #list for CEs in previous frame
+    currFrameCEs = []                  #list for CEs in current frame
+    cloudElementLat = []               #list for a particular CE's lat values
+    cloudElementLon = []               #list for a particular CE's lon values
+    cloudElementLatLons = []   #list for a particular CE's (lat,lon) values
+    
+    prevLatValue = 0.0
+    prevLonValue = 0.0
+    TIR_min = 0.0
+    TIR_max = 0.0
+    temporalRes = 3 # TRMM data is 3 hourly
+    precipTotal = 0.0
+    CETRMMList =[]
+    precip =[]
+    TRMMCloudElementLatLons =[]
+
+    minCELatLimit = 0.0
+    minCELonLimit = 0.0
+    maxCELatLimit = 0.0
+    maxCELonLimit = 0.0
+    
+    nygrd = len(LAT[:, 0]); nxgrd = len(LON[0, :])
+    
+    #openfile for storing ALL cloudElement information 
+    cloudElementsFile = 
open((MAINDIRECTORY+'/textFiles/cloudElements.txt'),'wb')
+    #openfile for storing cloudElement information meeting user criteria i.e. 
MCCs in this case
+    cloudElementsUserFile = 
open((MAINDIRECTORY+'/textFiles/cloudElementsUserFile.txt'),'w')
+    
+    #NB in the TRMM files the info is hours since the time thus 00Z file has 
in 01, 02 and 03 times
+    for t in xrange(mergImgs.shape[0]):
+        #-------------------------------------------------
+        # #textfile name for saving the data for arcgis
+        # thisFileName = MAINDIRECTORY+'/' + (str(timelist[t])).replace(" ", 
"_") + '.txt'
+        # cloudElementsTextFile = open(thisFileName,'w')
+        #-------------------------------------------------
+
+        #determine contiguous locations with temeperature below the warmest 
temp i.e. cloudElements in each frame
+        frame, CEcounter = ndimage.measurements.label(mergImgs[t,:,:], 
structure=STRUCTURING_ELEMENT)
+        frameCEcounter=0
+        frameNum += 1
+
+        #for each of the areas identified, check to determine if it a valid CE 
via an area and T requirement
+        for count in xrange(CEcounter):
+            #[0] is time dimension. Determine the actual values from the data
+            #loc is a masked array
+            try:
+                loc = ndimage.find_objects(frame==(count+1))[0]
+            except Exception, e:
+                print "Error is ", e
+                continue
+
+
+            cloudElement = mergImgs[t,:,:][loc]
+            labels, lcounter = ndimage.label(cloudElement)
+            
+            #determine the true lats and lons for this particular CE
+            cloudElementLat = LAT[loc[0],0]
+            cloudElementLon = LON[0,loc[1]] 
+            
+            #determine number of boxes in this cloudelement
+            numOfBoxes = np.count_nonzero(cloudElement)
+            cloudElementArea = numOfBoxes*XRES*YRES
+
+            #If the area is greater than the area required, or if the area is 
smaller than the suggested area, check if it meets a convective fraction 
requirement
+            #consider as CE
+
+            if cloudElementArea >= AREA_MIN or (cloudElementArea < AREA_MIN 
and ((ndimage.minimum(cloudElement, 
labels=labels))/float((ndimage.maximum(cloudElement, labels=labels)))) < 
CONVECTIVE_FRACTION ):
+
+                #get some time information and labeling info
+                frameTime = str(timelist[t])
+                frameCEcounter +=1
+                CEuniqueID = 'F'+str(frameNum)+'CE'+str(frameCEcounter) 
+
+                #-------------------------------------------------
+                #textfile name for accesing CE data using MATLAB code
+                # thisFileName = MAINDIRECTORY+'/' + 
(str(timelist[t])).replace(" ", "_") + CEuniqueID +'.txt'
+                # cloudElementsTextFile = open(thisFileName,'w')
+                #-------------------------------------------------
+
+                # ------ NETCDF File stuff for brightness temp stuff 
------------------------------------
+                thisFileName = MAINDIRECTORY +'/MERGnetcdfCEs/cloudElements'+ 
(str(timelist[t])).replace(" ", "_") + CEuniqueID +'.nc'
+                currNetCDFCEData = Dataset(thisFileName, 'w', format='NETCDF4')
+                currNetCDFCEData.description = 'Cloud Element '+CEuniqueID + ' 
temperature data'
+                currNetCDFCEData.calendar = 'standard'
+                currNetCDFCEData.conventions = 'COARDS'
+                # dimensions
+                currNetCDFCEData.createDimension('time', None)
+                currNetCDFCEData.createDimension('lat', len(LAT[:,0]))
+                currNetCDFCEData.createDimension('lon', len(LON[0,:]))
+                # variables
+                tempDims = ('time','lat', 'lon',)
+                times = currNetCDFCEData.createVariable('time', 'f8', 
('time',))
+                times.units = 'hours since '+ str(timelist[t])[:-6]
+                latitudes = currNetCDFCEData.createVariable('latitude', 'f8', 
('lat',))
+                longitudes = currNetCDFCEData.createVariable('longitude', 
'f8', ('lon',))
+                brightnesstemp = 
currNetCDFCEData.createVariable('brightnesstemp', 'i16',tempDims )
+                brightnesstemp.units = 'Kelvin'
+                # NETCDF data
+                dates=[timelist[t]+timedelta(hours=0)]
+                times[:] =  date2num(dates,units=times.units)
+                longitudes[:] = LON[0,:]
+                longitudes.units = "degrees_east" 
+                longitudes.long_name = "Longitude" 
+
+                latitudes[:] =  LAT[:,0]
+                latitudes.units = "degrees_north"
+                latitudes.long_name ="Latitude"
+                
+                #generate array of zeros for brightness temperature
+                brightnesstemp1 = ma.zeros((1,len(latitudes), 
len(longitudes))).astype('int16')
+                #-----------End most of NETCDF file stuff 
------------------------------------
+
+                #if other dataset (TRMM) assumed to be a precipitation dataset 
was entered
+                if TRMMdirName:
+                    #------------------TRMM stuff 
-------------------------------------------------
+                    fileDate = ((str(timelist[t])).replace(" ", 
"")[:-8]).replace("-","")
+                    fileHr1 = (str(timelist[t])).replace(" ", "")[-8:-6]
+                    
+                    if int(fileHr1) % temporalRes == 0:
+                        fileHr = fileHr1
+                    else:
+                        fileHr = (int(fileHr1)/temporalRes) * temporalRes
+                    if fileHr < 10:
+                        fileHr = '0'+str(fileHr)
+                    else:
+                        str(fileHr)
+
+                    #open TRMM file for the resolution info and to create the 
appropriate sized grid
+                    TRMMfileName = TRMMdirName+'/3B42.'+ fileDate + 
"."+str(fileHr)+".7A.nc"
+                    
+                    TRMMData = Dataset(TRMMfileName,'r', format='NETCDF4')
+                    precipRate = TRMMData.variables['pcp'][:,:,:]
+                    latsrawTRMMData = TRMMData.variables['latitude'][:]
+                    lonsrawTRMMData = TRMMData.variables['longitude'][:]
+                    lonsrawTRMMData[lonsrawTRMMData > 180] = 
lonsrawTRMMData[lonsrawTRMMData>180] - 360.
+                    LONTRMM, LATTRMM = np.meshgrid(lonsrawTRMMData, 
latsrawTRMMData)
+
+                    nygrdTRMM = len(LATTRMM[:,0]); nxgrdTRMM = 
len(LONTRMM[0,:])
+                    precipRateMasked = ma.masked_array(precipRate, 
mask=(precipRate < 0.0))
+                    #---------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)
+                    
#----------------------------------------------------------------------------------
+        
+                    # #get the lat/lon info from cloudElement
+                    #get the lat/lon info from the file
+                    latCEStart = LAT[0][0]
+                    latCEEnd = LAT[-1][0]
+                    lonCEStart = LON[0][0]
+                    lonCEEnd = LON[0][-1]
+                    
+                    #get the lat/lon info for TRMM data (different resolution)
+                    latStartT = find_nearest(latsrawTRMMData, latCEStart)
+                    latEndT = find_nearest(latsrawTRMMData, latCEEnd)
+                    lonStartT = find_nearest(lonsrawTRMMData, lonCEStart)
+                    lonEndT = find_nearest(lonsrawTRMMData, lonCEEnd)
+                    latStartIndex = np.where(latsrawTRMMData == latStartT)
+                    latEndIndex = np.where(latsrawTRMMData == latEndT)
+                    lonStartIndex = np.where(lonsrawTRMMData == lonStartT)
+                    lonEndIndex = np.where(lonsrawTRMMData == lonEndT)
+
+                    #get the relevant TRMM info 
+                    CEprecipRate = 
precipRate[:,(latStartIndex[0][0]-1):latEndIndex[0][0],(lonStartIndex[0][0]-1):lonEndIndex[0][0]]
+                    TRMMData.close()
+                    
+                    # ------ NETCDF File info for writing TRMM CE rainfall 
------------------------------------
+                    thisFileName = MAINDIRECTORY+'/TRMMnetcdfCEs/TRMM' + 
(str(timelist[t])).replace(" ", "_") + CEuniqueID +'.nc'
+                    currNetCDFTRMMData = Dataset(thisFileName, 'w', 
format='NETCDF4')
+                    currNetCDFTRMMData.description = 'Cloud Element 
'+CEuniqueID + ' precipitation data'
+                    currNetCDFTRMMData.calendar = 'standard'
+                    currNetCDFTRMMData.conventions = 'COARDS'
+                    # dimensions
+                    currNetCDFTRMMData.createDimension('time', None)
+                    currNetCDFTRMMData.createDimension('lat', len(LAT[:,0]))
+                    currNetCDFTRMMData.createDimension('lon', len(LON[0,:]))
+                    
+                    # variables
+                    TRMMprecip = ('time','lat', 'lon',)
+                    times = currNetCDFTRMMData.createVariable('time', 'f8', 
('time',))
+                    times.units = 'hours since '+ str(timelist[t])[:-6]
+                    latitude = currNetCDFTRMMData.createVariable('latitude', 
'f8', ('lat',))
+                    longitude = currNetCDFTRMMData.createVariable('longitude', 
'f8', ('lon',))
+                    rainFallacc = 
currNetCDFTRMMData.createVariable('precipitation_Accumulation', 'f8',TRMMprecip 
)
+                    rainFallacc.units = 'mm'
+
+                    longitude[:] = LON[0,:]
+                    longitude.units = "degrees_east" 
+                    longitude.long_name = "Longitude" 
+
+                    latitude[:] =  LAT[:,0]
+                    latitude.units = "degrees_north"
+                    latitude.long_name ="Latitude"
+
+                    finalCETRMMvalues = ma.zeros((brightnesstemp.shape))
+                    #-----------End most of NETCDF file stuff 
------------------------------------
+
+                #populate cloudElementLatLons by unpacking the original values 
from loc to get the actual value for lat and lon
+                #TODO: KDW - too dirty... play with itertools.izip or zip and 
the enumerate with this
+                #                      as cloudElement is masked
+                for index,value in np.ndenumerate(cloudElement):
+                    if value != 0 : 
+                        lat_index,lon_index = index
+                        lat_lon_tuple = (cloudElementLat[lat_index], 
cloudElementLon[lon_index],value)
+
+                        #generate the comma separated file for GIS
+                        cloudElementLatLons.append(lat_lon_tuple)
+
+                        #temp data for CE NETCDF file
+                        
brightnesstemp1[0,int(np.where(LAT[:,0]==cloudElementLat[lat_index])[0]),int(np.where(LON[0,:]==cloudElementLon[lon_index])[0])]
 = value
+                        
+                        if TRMMdirName:
+                            
finalCETRMMvalues[0,int(np.where(LAT[:,0]==cloudElementLat[lat_index])[0]),int(np.where(LON[0,:]==cloudElementLon[lon_index])[0])]
 = 
regriddedTRMM[int(np.where(LAT[:,0]==cloudElementLat[lat_index])[0]),int(np.where(LON[0,:]==cloudElementLon[lon_index])[0])]
+                            CETRMMList.append((cloudElementLat[lat_index], 
cloudElementLon[lon_index], finalCETRMMvalues[0,cloudElementLat[lat_index], 
cloudElementLon[lon_index]]))
+
+
+                brightnesstemp[:] = brightnesstemp1[:]
+                currNetCDFCEData.close()
+
+                if TRMMdirName:
+
+                    #calculate the total precip associated with the feature
+                    for index, value in np.ndenumerate(finalCETRMMvalues):
+                        precipTotal += value 
+                        precip.append(value)
+            
+                    rainFallacc[:] = finalCETRMMvalues[:]
+                    currNetCDFTRMMData.close()
+                    TRMMnumOfBoxes = np.count_nonzero(finalCETRMMvalues)
+                    TRMMArea = TRMMnumOfBoxes*XRES*YRES
+                    try:
+                        maxCEprecipRate = 
np.max(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
+                        minCEprecipRate = 
np.min(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
+                    except:
+                        pass
+
+                #sort cloudElementLatLons by lats
+                cloudElementLatLons.sort(key=lambda tup: tup[0])       
+
+                #determine if the cloud element the shape 
+                cloudElementEpsilon = eccentricity (cloudElement)
+                cloudElementsUserFile.write("\n\nTime is: %s" 
%(str(timelist[t])))
+                cloudElementsUserFile.write("\nCEuniqueID is: %s" %CEuniqueID)
+                latCenter, lonCenter = 
ndimage.measurements.center_of_mass(cloudElement, labels=labels)
+                
+                #latCenter and lonCenter are given according to the particular 
array defining this CE
+                #so you need to convert this value to the overall domain truth
+                latCenter = cloudElementLat[round(latCenter)]
+                lonCenter = cloudElementLon[round(lonCenter)]
+                cloudElementsUserFile.write("\nCenter (lat,lon) is: 
%.2f\t%.2f" %(latCenter, lonCenter))
+                cloudElementCenter.append(latCenter)
+                cloudElementCenter.append(lonCenter)
+                cloudElementsUserFile.write("\nNumber of boxes are: %d" 
%numOfBoxes)
+                cloudElementsUserFile.write("\nArea is: %.4f km^2" 
%(cloudElementArea))
+                cloudElementsUserFile.write("\nAverage brightness temperature 
is: %.4f K" %ndimage.mean(cloudElement, labels=labels))
+                cloudElementsUserFile.write("\nMin brightness temperature is: 
%.4f K" %ndimage.minimum(cloudElement, labels=labels))
+                cloudElementsUserFile.write("\nMax brightness temperature is: 
%.4f K" %ndimage.maximum(cloudElement, labels=labels))
+                cloudElementsUserFile.write("\nBrightness temperature variance 
is: %.4f K" %ndimage.variance(cloudElement, labels=labels))
+                cloudElementsUserFile.write("\nConvective fraction is: %.4f " 
%(((ndimage.minimum(cloudElement, 
labels=labels))/float((ndimage.maximum(cloudElement, labels=labels))))*100.0))
+                cloudElementsUserFile.write("\nEccentricity is: %.4f " 
%(cloudElementEpsilon))
+                #populate the dictionary
+                if TRMMdirName:
+                    cloudElementDict = {'uniqueID': CEuniqueID, 
'cloudElementTime': timelist[t],'cloudElementLatLon': cloudElementLatLons, 
'cloudElementCenter':cloudElementCenter, 'cloudElementArea':cloudElementArea, 
'cloudElementEccentricity':cloudElementEpsilon, 'cloudElementTmax':TIR_max, 
'cloudElementTmin': TIR_min, 
'cloudElementPrecipTotal':precipTotal,'cloudElementLatLonTRMM':CETRMMList, 
'TRMMArea': TRMMArea,'CETRMMmax':maxCEprecipRate, 'CETRMMmin':minCEprecipRate}
+                else:
+                    cloudElementDict = {'uniqueID': CEuniqueID, 
'cloudElementTime': timelist[t],'cloudElementLatLon': cloudElementLatLons, 
'cloudElementCenter':cloudElementCenter, 'cloudElementArea':cloudElementArea, 
'cloudElementEccentricity':cloudElementEpsilon, 'cloudElementTmax':TIR_max, 
'cloudElementTmin': TIR_min,}
+                
+                #current frame list of CEs
+                currFrameCEs.append(cloudElementDict)
+                
+                #draw the graph node
+                CLOUD_ELEMENT_GRAPH.add_node(CEuniqueID, cloudElementDict)
+                
+                if frameNum != 1:
+                    for cloudElementDict in prevFrameCEs:
+                        thisCElen = len(cloudElementLatLons)
+                        percentageOverlap, areaOverlap = 
cloudElementOverlap(cloudElementLatLons, cloudElementDict['cloudElementLatLon'])
+                        
+                        #change weights to integers because the built in 
shortest path chokes on floating pts according to Networkx doc
+                        #according to Goyens et al, two CEs are considered 
related if there is atleast 95% overlap between them for consecutive imgs a max 
of 2 hrs apart
+                        if percentageOverlap >= 0.95: 
+                            
CLOUD_ELEMENT_GRAPH.add_edge(cloudElementDict['uniqueID'], CEuniqueID, 
weight=edgeWeight[0])
+                            
+                        elif percentageOverlap >= 0.90 and percentageOverlap < 
0.95 :
+                            
CLOUD_ELEMENT_GRAPH.add_edge(cloudElementDict['uniqueID'], CEuniqueID, 
weight=edgeWeight[1])
+
+                        elif areaOverlap >= MIN_OVERLAP:
+                            
CLOUD_ELEMENT_GRAPH.add_edge(cloudElementDict['uniqueID'], CEuniqueID, 
weight=edgeWeight[2])
+
+                else:
+                    #TODO: remove this else as we only wish for the CE details
+                    #ensure only the non-zero elements are considered
+                    #store intel in allCE file
+                    labels, _ = ndimage.label(cloudElement)
+                    
cloudElementsFile.write("\n-----------------------------------------------")
+                    cloudElementsFile.write("\n\nTime is: %s" 
%(str(timelist[t])))
+                    # cloudElementLat = LAT[loc[0],0]
+                    # cloudElementLon = LON[0,loc[1]] 
+                    
+                    #populate cloudElementLatLons by unpacking the original 
values from loc
+                    #TODO: KDW - too dirty... play with itertools.izip or zip 
and the enumerate with this
+                    #                  as cloudElement is masked
+                    for index,value in np.ndenumerate(cloudElement):
+                        if value != 0 : 
+                            lat_index,lon_index = index
+                            lat_lon_tuple = (cloudElementLat[lat_index], 
cloudElementLon[lon_index])
+                            cloudElementLatLons.append(lat_lon_tuple)
+    
+                    cloudElementsFile.write("\nLocation of rejected CE 
(lat,lon) points are: %s" %cloudElementLatLons)
+                    #latCenter and lonCenter are given according to the 
particular array defining this CE
+                    #so you need to convert this value to the overall domain 
truth
+                    latCenter, lonCenter = 
ndimage.measurements.center_of_mass(cloudElement, labels=labels)
+                    latCenter = cloudElementLat[round(latCenter)]
+                    lonCenter = cloudElementLon[round(lonCenter)]
+                    cloudElementsFile.write("\nCenter (lat,lon) is: 
%.2f\t%.2f" %(latCenter, lonCenter))
+                    cloudElementsFile.write("\nNumber of boxes are: %d" 
%numOfBoxes)
+                    cloudElementsFile.write("\nArea is: %.4f km^2" 
%(cloudElementArea))
+                    cloudElementsFile.write("\nAverage brightness temperature 
is: %.4f K" %ndimage.mean(cloudElement, labels=labels))
+                    cloudElementsFile.write("\nMin brightness temperature is: 
%.4f K" %ndimage.minimum(cloudElement, labels=labels))
+                    cloudElementsFile.write("\nMax brightness temperature is: 
%.4f K" %ndimage.maximum(cloudElement, labels=labels))
+                    cloudElementsFile.write("\nBrightness temperature variance 
is: %.4f K" %ndimage.variance(cloudElement, labels=labels))
+                    cloudElementsFile.write("\nConvective fraction is: %.4f " 
%(((ndimage.minimum(cloudElement, 
labels=labels))/float((ndimage.maximum(cloudElement, labels=labels))))*100.0))
+                    cloudElementsFile.write("\nEccentricity is: %.4f " 
%(cloudElementEpsilon))
+                    
cloudElementsFile.write("\n-----------------------------------------------")
+                    
+            #reset list for the next CE
+            nodeExist = False
+            cloudElementCenter=[]
+            cloudElement = []
+            cloudElementLat=[]
+            cloudElementLon =[]
+            cloudElementLatLons =[]
+            brightnesstemp1 =[]
+            brightnesstemp =[]
+            finalCETRMMvalues =[]
+            CEprecipRate =[]
+            CETRMMList =[]
+            precipTotal = 0.0
+            precip=[]
+            TRMMCloudElementLatLons=[]
+            
+        #reset for the next time
+        prevFrameCEs =[]
+        prevFrameCEs = currFrameCEs
+        currFrameCEs =[]
+                        
+    cloudElementsFile.close
+    cloudElementsUserFile.close
+    #if using ARCGIS data store code, uncomment this file close line
+    #cloudElementsTextFile.close
+
+    #clean up graph - remove parent and childless nodes
+    outAndInDeg = CLOUD_ELEMENT_GRAPH.degree_iter()
+    toRemove = [node[0] for node in outAndInDeg if node[1]<1]
+    CLOUD_ELEMENT_GRAPH.remove_nodes_from(toRemove)
+    
+    print "number of nodes are: ", CLOUD_ELEMENT_GRAPH.number_of_nodes()
+    print "number of edges are: ", CLOUD_ELEMENT_GRAPH.number_of_edges()
+    print ("*"*80)
+
+    #hierachial graph output
+    graphTitle = "Cloud Elements observed over somewhere from 0000Z to 0000Z" 
+    #drawGraph(CLOUD_ELEMENT_GRAPH, graphTitle, edgeWeight)
+
+    return CLOUD_ELEMENT_GRAPH 
 #******************************************************************
 def findPrecipRate(TRMMdirName, timelist):
-       ''' 
-       Purpose:: 
-               Determines the precipitation rates for MCSs found if 
TRMMdirName was not entered in findCloudElements this can be used
+    ''' 
+    Purpose:: 
+        Determines the precipitation rates for MCSs found if TRMMdirName was 
not entered in findCloudElements this can be used
 
-       Input:: 
-               TRMMdirName: a string representing the directory for the 
original TRMM netCDF files
-               timelist: a list of python datatimes
+    Input:: 
+        TRMMdirName: a string representing the directory for the original TRMM 
netCDF files
+        timelist: a list of python datatimes
 
-       Output:: a list of dictionary of the TRMM data 
-               NB: also creates netCDF with TRMM data for each CE (for post 
processing) index
-                       in MAINDIRECTORY/TRMMnetcdfCEs
+    Output:: a list of dictionary of the TRMM data 
+        NB: also creates netCDF with TRMM data for each CE (for post 
processing) index
+            in MAINDIRECTORY/TRMMnetcdfCEs
    
-       Assumptions:: Assumes that findCloudElements was run without the 
TRMMdirName value 
+    Assumptions:: Assumes that findCloudElements was run without the 
TRMMdirName value 
  
-       '''
-       allCEnodesTRMMdata =[]
-       TRMMdataDict={}
-       precipTotal = 0.0
-
-       os.chdir((MAINDIRECTORY+'/MERGnetcdfCEs/'))
-       imgFilename = ''
-       temporalRes = 3 #3 hours for TRMM
-       
-       #sort files
-       files = filter(os.path.isfile, glob.glob("*.nc"))
-       files.sort(key=lambda x: os.path.getmtime(x))
-       
-       for afile in files:
-               fullFname = os.path.splitext(afile)[0]
-               noFrameExtension = (fullFname.replace("_","")).split('F')[0]
-               CEuniqueID = 'F' +(fullFname.replace("_","")).split('F')[1]
-               fileDateTimeChar = 
(noFrameExtension.replace(":","")).split('s')[1]
-               fileDateTime = fileDateTimeChar.replace("-","")
-               fileDate = fileDateTime[:-6]
-               fileHr1=fileDateTime[-6:-4]
-
-               cloudElementData = Dataset(afile,'r', format='NETCDF4')
-               brightnesstemp1 = 
cloudElementData.variables['brightnesstemp'][:,:,:] 
-               latsrawCloudElements = cloudElementData.variables['latitude'][:]
-               lonsrawCloudElements = 
cloudElementData.variables['longitude'][:]
-               
-               brightnesstemp = np.squeeze(brightnesstemp1, axis=0)
-               
-               if int(fileHr1) % temporalRes == 0:
-                       fileHr = fileHr1
-               else:
-                       fileHr = (int(fileHr1)/temporalRes) * temporalRes
-               
-               if fileHr < 10:
-                       fileHr = '0'+str(fileHr)
-               else:
-                       str(fileHr)
-
-               TRMMfileName = TRMMdirName+"/3B42."+ str(fileDate) + 
"."+str(fileHr)+".7A.nc"
-               TRMMData = Dataset(TRMMfileName,'r', format='NETCDF4')
-               precipRate = TRMMData.variables['pcp'][:,:,:]
-               latsrawTRMMData = TRMMData.variables['latitude'][:]
-               lonsrawTRMMData = TRMMData.variables['longitude'][:]
-               lonsrawTRMMData[lonsrawTRMMData > 180] = 
lonsrawTRMMData[lonsrawTRMMData>180] - 360.
-               LONTRMM, LATTRMM = np.meshgrid(lonsrawTRMMData, latsrawTRMMData)
-
-               #nygrdTRMM = len(LATTRMM[:,0]); nxgrd = len(LONTRMM[0,:])
-               nygrd = len(LAT[:, 0]); nxgrd = len(LON[0, :])
-
-               precipRateMasked = ma.masked_array(precipRate, mask=(precipRate 
< 0.0))
-               #---------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 = 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
-               latCEStart = LAT[0][0]
-               latCEEnd = LAT[-1][0]
-               lonCEStart = LON[0][0]
-               lonCEEnd = LON[0][-1]
-
-               #get the lat/lon info for TRMM data (different resolution)
-               latStartT = find_nearest(latsrawTRMMData, latCEStart)
-               latEndT = find_nearest(latsrawTRMMData, latCEEnd)
-               lonStartT = find_nearest(lonsrawTRMMData, lonCEStart)
-               lonEndT = find_nearest(lonsrawTRMMData, lonCEEnd)
-               latStartIndex = np.where(latsrawTRMMData == latStartT)
-               latEndIndex = np.where(latsrawTRMMData == latEndT)
-               lonStartIndex = np.where(lonsrawTRMMData == lonStartT)
-               lonEndIndex = np.where(lonsrawTRMMData == lonEndT)
-
-               #get the relevant TRMM info 
-               CEprecipRate = 
precipRate[:,(latStartIndex[0][0]-1):latEndIndex[0][0],(lonStartIndex[0][0]-1):lonEndIndex[0][0]]
-               TRMMData.close()
-                       
-               
-               # ------ NETCDF File stuff ------------------------------------
-               thisFileName = MAINDIRECTORY+'/TRMMnetcdfCEs/'+ fileDateTime + 
CEuniqueID +'.nc'
-               currNetCDFTRMMData = Dataset(thisFileName, 'w', 
format='NETCDF4')
-               currNetCDFTRMMData.description = 'Cloud Element '+CEuniqueID + 
' rainfall data'
-               currNetCDFTRMMData.calendar = 'standard'
-               currNetCDFTRMMData.conventions = 'COARDS'
-               # dimensions
-               currNetCDFTRMMData.createDimension('time', None)
-               currNetCDFTRMMData.createDimension('lat', len(LAT[:,0]))
-               currNetCDFTRMMData.createDimension('lon', len(LON[0,:]))
-               # variables
-               TRMMprecip = ('time','lat', 'lon',)
-               times = currNetCDFTRMMData.createVariable('time', 'f8', 
('time',))
-               times.units = 'hours since '+ fileDateTime[:-6] 
-               latitude = currNetCDFTRMMData.createVariable('latitude', 'f8', 
('lat',))
-               longitude = currNetCDFTRMMData.createVariable('longitude', 
'f8', ('lon',))
-               rainFallacc = 
currNetCDFTRMMData.createVariable('precipitation_Accumulation', 'f8',TRMMprecip 
)
-               rainFallacc.units = 'mm'
-
-               longitude[:] = LON[0,:]
-               longitude.units = "degrees_east" 
-               longitude.long_name = "Longitude" 
-
-               latitude[:] =  LAT[:,0]
-               latitude.units = "degrees_north"
-               latitude.long_name ="Latitude"
-
-               finalCETRMMvalues = ma.zeros((brightnesstemp1.shape))
-               #-----------End most of NETCDF file stuff 
------------------------------------  
-               for index,value in np.ndenumerate(brightnesstemp):
-                       lat_index, lon_index = index
-                       currTimeValue = 0
-                       if value > 0:
-
-                               finalCETRMMvalues[0,lat_index,lon_index] = 
regriddedTRMM[int(np.where(LAT[:,0]==LAT[lat_index,0])[0]), 
int(np.where(LON[0,:]==LON[0,lon_index])[0])]
-                               
-
-               rainFallacc[:] = finalCETRMMvalues
-               currNetCDFTRMMData.close()
-
-               for index, value in np.ndenumerate(finalCETRMMvalues):
-                       precipTotal += value 
-
-               TRMMnumOfBoxes = np.count_nonzero(finalCETRMMvalues)
-               TRMMArea = TRMMnumOfBoxes*XRES*YRES     
-
-               try:
-                       minCEprecipRate = 
np.min(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
-               except:
-                       minCEprecipRate = 0.0
-
-               try:
-                       maxCEprecipRate = 
np.max(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
-               except:
-                       maxCEprecipRate = 0.0
-
-               #add info to CLOUDELEMENTSGRAPH
-               #TODO try block
-               for eachdict in CLOUD_ELEMENT_GRAPH.nodes(CEuniqueID):
-                       if eachdict[1]['uniqueID'] == CEuniqueID:
-                               if not 'cloudElementPrecipTotal' in 
eachdict[1].keys():
-                                       eachdict[1]['cloudElementPrecipTotal'] 
= precipTotal
-                               if not 'cloudElementLatLonTRMM' in 
eachdict[1].keys():
-                                       eachdict[1]['cloudElementLatLonTRMM'] = 
finalCETRMMvalues
-                               if not 'TRMMArea' in eachdict[1].keys():
-                                       eachdict[1]['TRMMArea'] = TRMMArea
-                               if not 'CETRMMmin' in eachdict[1].keys():
-                                       eachdict[1]['CETRMMmin'] = 
minCEprecipRate
-                               if not 'CETRMMmax' in eachdict[1].keys():
-                                       eachdict[1]['CETRMMmax'] = 
maxCEprecipRate
-
-               #clean up
-               precipTotal = 0.0
-               latsrawTRMMData =[]
-               lonsrawTRMMData = []
-               latsrawCloudElements=[]
-               lonsrawCloudElements=[]
-               finalCETRMMvalues =[]
-               CEprecipRate =[]
-               brightnesstemp =[]
-               TRMMdataDict ={}
-
-       return allCEnodesTRMMdata
+    '''
+    allCEnodesTRMMdata =[]
+    TRMMdataDict={}
+    precipTotal = 0.0
+
+    os.chdir((MAINDIRECTORY+'/MERGnetcdfCEs/'))
+    imgFilename = ''
+    temporalRes = 3 #3 hours for TRMM
+    
+    #sort files
+    files = filter(os.path.isfile, glob.glob("*.nc"))
+    files.sort(key=lambda x: os.path.getmtime(x))
+    
+    for afile in files:
+        fullFname = os.path.splitext(afile)[0]
+        noFrameExtension = (fullFname.replace("_","")).split('F')[0]
+        CEuniqueID = 'F' +(fullFname.replace("_","")).split('F')[1]
+        fileDateTimeChar = (noFrameExtension.replace(":","")).split('s')[1]
+        fileDateTime = fileDateTimeChar.replace("-","")
+        fileDate = fileDateTime[:-6]
+        fileHr1=fileDateTime[-6:-4]
+
+        cloudElementData = Dataset(afile,'r', format='NETCDF4')
+        brightnesstemp1 = cloudElementData.variables['brightnesstemp'][:,:,:] 
+        latsrawCloudElements = cloudElementData.variables['latitude'][:]
+        lonsrawCloudElements = cloudElementData.variables['longitude'][:]
+        
+        brightnesstemp = np.squeeze(brightnesstemp1, axis=0)
+        
+        if int(fileHr1) % temporalRes == 0:
+            fileHr = fileHr1
+        else:
+            fileHr = (int(fileHr1)/temporalRes) * temporalRes
+        
+        if fileHr < 10:
+            fileHr = '0'+str(fileHr)
+        else:
+            str(fileHr)
+
+        TRMMfileName = TRMMdirName+"/3B42."+ str(fileDate) + 
"."+str(fileHr)+".7A.nc"
+        TRMMData = Dataset(TRMMfileName,'r', format='NETCDF4')
+        precipRate = TRMMData.variables['pcp'][:,:,:]
+        latsrawTRMMData = TRMMData.variables['latitude'][:]
+        lonsrawTRMMData = TRMMData.variables['longitude'][:]
+        lonsrawTRMMData[lonsrawTRMMData > 180] = 
lonsrawTRMMData[lonsrawTRMMData>180] - 360.
+        LONTRMM, LATTRMM = np.meshgrid(lonsrawTRMMData, latsrawTRMMData)
+
+        #nygrdTRMM = len(LATTRMM[:,0]); nxgrd = len(LONTRMM[0,:])
+        nygrd = len(LAT[:, 0]); nxgrd = len(LON[0, :])
+
+        precipRateMasked = ma.masked_array(precipRate, mask=(precipRate < 0.0))
+        #---------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 = 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
+        latCEStart = LAT[0][0]
+        latCEEnd = LAT[-1][0]
+        lonCEStart = LON[0][0]
+        lonCEEnd = LON[0][-1]
+
+        #get the lat/lon info for TRMM data (different resolution)
+        latStartT = find_nearest(latsrawTRMMData, latCEStart)
+        latEndT = find_nearest(latsrawTRMMData, latCEEnd)
+        lonStartT = find_nearest(lonsrawTRMMData, lonCEStart)
+        lonEndT = find_nearest(lonsrawTRMMData, lonCEEnd)
+        latStartIndex = np.where(latsrawTRMMData == latStartT)
+        latEndIndex = np.where(latsrawTRMMData == latEndT)
+        lonStartIndex = np.where(lonsrawTRMMData == lonStartT)
+        lonEndIndex = np.where(lonsrawTRMMData == lonEndT)
+
+        #get the relevant TRMM info 
+        CEprecipRate = 
precipRate[:,(latStartIndex[0][0]-1):latEndIndex[0][0],(lonStartIndex[0][0]-1):lonEndIndex[0][0]]
+        TRMMData.close()
+            
+        
+        # ------ NETCDF File stuff ------------------------------------
+        thisFileName = MAINDIRECTORY+'/TRMMnetcdfCEs/'+ fileDateTime + 
CEuniqueID +'.nc'
+        currNetCDFTRMMData = Dataset(thisFileName, 'w', format='NETCDF4')
+        currNetCDFTRMMData.description = 'Cloud Element '+CEuniqueID + ' 
rainfall data'
+        currNetCDFTRMMData.calendar = 'standard'
+        currNetCDFTRMMData.conventions = 'COARDS'
+        # dimensions
+        currNetCDFTRMMData.createDimension('time', None)
+        currNetCDFTRMMData.createDimension('lat', len(LAT[:,0]))
+        currNetCDFTRMMData.createDimension('lon', len(LON[0,:]))
+        # variables
+        TRMMprecip = ('time','lat', 'lon',)
+        times = currNetCDFTRMMData.createVariable('time', 'f8', ('time',))
+        times.units = 'hours since '+ fileDateTime[:-6] 
+        latitude = currNetCDFTRMMData.createVariable('latitude', 'f8', 
('lat',))
+        longitude = currNetCDFTRMMData.createVariable('longitude', 'f8', 
('lon',))
+        rainFallacc = 
currNetCDFTRMMData.createVariable('precipitation_Accumulation', 'f8',TRMMprecip 
)
+        rainFallacc.units = 'mm'
+
+        longitude[:] = LON[0,:]
+        longitude.units = "degrees_east" 
+        longitude.long_name = "Longitude" 
+
+        latitude[:] =  LAT[:,0]
+        latitude.units = "degrees_north"
+        latitude.long_name ="Latitude"
+
+        finalCETRMMvalues = ma.zeros((brightnesstemp1.shape))
+        #-----------End most of NETCDF file stuff 
------------------------------------ 
+        for index,value in np.ndenumerate(brightnesstemp):
+            lat_index, lon_index = index
+            currTimeValue = 0
+            if value > 0:
+
+                finalCETRMMvalues[0,lat_index,lon_index] = 
regriddedTRMM[int(np.where(LAT[:,0]==LAT[lat_index,0])[0]), 
int(np.where(LON[0,:]==LON[0,lon_index])[0])]
+                
+
+        rainFallacc[:] = finalCETRMMvalues
+        currNetCDFTRMMData.close()
+
+        for index, value in np.ndenumerate(finalCETRMMvalues):
+            precipTotal += value 
+
+        TRMMnumOfBoxes = np.count_nonzero(finalCETRMMvalues)
+        TRMMArea = TRMMnumOfBoxes*XRES*YRES    
+
+        try:
+            minCEprecipRate = 
np.min(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
+        except:
+            minCEprecipRate = 0.0
+
+        try:
+            maxCEprecipRate = 
np.max(finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
+        except:
+            maxCEprecipRate = 0.0
+
+        #add info to CLOUDELEMENTSGRAPH
+        #TODO try block
+        for eachdict in CLOUD_ELEMENT_GRAPH.nodes(CEuniqueID):
+            if eachdict[1]['uniqueID'] == CEuniqueID:
+                if not 'cloudElementPrecipTotal' in eachdict[1].keys():
+                    eachdict[1]['cloudElementPrecipTotal'] = precipTotal
+                if not 'cloudElementLatLonTRMM' in eachdict[1].keys():
+                    eachdict[1]['cloudElementLatLonTRMM'] = finalCETRMMvalues
+                if not 'TRMMArea' in eachdict[1].keys():
+                    eachdict[1]['TRMMArea'] = TRMMArea
+                if not 'CETRMMmin' in eachdict[1].keys():
+                    eachdict[1]['CETRMMmin'] = minCEprecipRate
+                if not 'CETRMMmax' in eachdict[1].keys():
+                    eachdict[1]['CETRMMmax'] = maxCEprecipRate
+
+        #clean up
+        precipTotal = 0.0
+        latsrawTRMMData =[]
+        lonsrawTRMMData = []
+        latsrawCloudElements=[]
+        lonsrawCloudElements=[]
+        finalCETRMMvalues =[]
+        CEprecipRate =[]
+        brightnesstemp =[]
+        TRMMdataDict ={}
+
+    return allCEnodesTRMMdata
 #******************************************************************    
 def findCloudClusters(CEGraph):
-       '''
-       Purpose:: 
-               Determines the cloud clusters properties from the subgraphs in 
-               the graph i.e. prunes the graph according to the minimum depth
-
-       Input:: 
-               CEGraph: a Networkx directed graph of the CEs with weighted 
edges
-               according the area overlap between nodes (CEs) of consectuive 
frames
-       
-       Output:: 
-               PRUNED_GRAPH: a Networkx directed graph of with CCs/ MCSs
-
-       '''
-
-       seenNode = []
-       allMCSLists =[]
-       pathDictList =[]
-       pathList=[]
-
-       cloudClustersFile = 
open((MAINDIRECTORY+'/textFiles/cloudClusters.txt'),'wb')
-       
-       for eachNode in CEGraph:
-               #check if the node has been seen before
-               if eachNode not in dict(enumerate(zip(*seenNode))):
-                       #look for all trees associated with node as the root
-                       thisPathDistanceAndLength = 
nx.single_source_dijkstra(CEGraph, eachNode)
-                       #determine the actual shortestPath and minimum 
depth/length
-                       maxDepthAndMinPath = 
findMaxDepthAndMinPath(thisPathDistanceAndLength)
-                       if maxDepthAndMinPath:
-                               maxPathLength = maxDepthAndMinPath[0] 
-                               shortestPath = maxDepthAndMinPath[1]
-                               
-                               #add nodes and paths to PRUNED_GRAPH
-                               for i in xrange(len(shortestPath)):
-                                       if 
PRUNED_GRAPH.has_node(shortestPath[i]) is False:
-                                               
PRUNED_GRAPH.add_node(shortestPath[i])
-                                               
-                                       #add edge if necessary
-                                       if i < (len(shortestPath)-1) and 
PRUNED_GRAPH.has_edge(shortestPath[i], shortestPath[i+1]) is False:
-                                               prunedGraphEdgeweight = 
CEGraph.get_edge_data(shortestPath[i], shortestPath[i+1])['weight']
-                                               
PRUNED_GRAPH.add_edge(shortestPath[i], shortestPath[i+1], 
weight=prunedGraphEdgeweight)
-
-                               #note information in a file for consideration 
later i.e. checking to see if it works
-                               cloudClustersFile.write("\nSubtree pathlength 
is %d and path is %s" %(maxPathLength, shortestPath))
-                               #update seenNode info
-                               seenNode.append(shortestPath)   
-
-       print "pruned graph"
-       print "number of nodes are: ", PRUNED_GRAPH.number_of_nodes()
-       print "number of edges are: ", PRUNED_GRAPH.number_of_edges()
-       print ("*"*80)          
-                                       
-       graphTitle = "Cloud Clusters observed over somewhere during sometime"
-       #drawGraph(PRUNED_GRAPH, graphTitle, edgeWeight)
-       cloudClustersFile.close
-       
-       return PRUNED_GRAPH  
+    '''
+    Purpose:: 
+        Determines the cloud clusters properties from the subgraphs in 
+        the graph i.e. prunes the graph according to the minimum depth
+
+    Input:: 
+        CEGraph: a Networkx directed graph of the CEs with weighted edges
+        according the area overlap between nodes (CEs) of consectuive frames
+    
+    Output:: 
+        PRUNED_GRAPH: a Networkx directed graph of with CCs/ MCSs
+
+    '''
+
+    seenNode = []
+    allMCSLists =[]
+    pathDictList =[]
+    pathList=[]
+
+    cloudClustersFile = 
open((MAINDIRECTORY+'/textFiles/cloudClusters.txt'),'wb')
+    
+    for eachNode in CEGraph:
+        #check if the node has been seen before
+        if eachNode not in dict(enumerate(zip(*seenNode))):
+            #look for all trees associated with node as the root
+            thisPathDistanceAndLength = nx.single_source_dijkstra(CEGraph, 
eachNode)
+            #determine the actual shortestPath and minimum depth/length
+            maxDepthAndMinPath = 
findMaxDepthAndMinPath(thisPathDistanceAndLength)
+            if maxDepthAndMinPath:
+                maxPathLength = maxDepthAndMinPath[0] 
+                shortestPath = maxDepthAndMinPath[1]
+                
+                #add nodes and paths to PRUNED_GRAPH
+                for i in xrange(len(shortestPath)):
+                    if PRUNED_GRAPH.has_node(shortestPath[i]) is False:
+                        PRUNED_GRAPH.add_node(shortestPath[i])
+                        
+                    #add edge if necessary
+                    if i < (len(shortestPath)-1) and 
PRUNED_GRAPH.has_edge(shortestPath[i], shortestPath[i+1]) is False:
+                        prunedGraphEdgeweight = 
CEGraph.get_edge_data(shortestPath[i], shortestPath[i+1])['weight']
+                        PRUNED_GRAPH.add_edge(shortestPath[i], 
shortestPath[i+1], weight=prunedGraphEdgeweight)
+
+                #note information in a file for consideration later i.e. 
checking to see if it works
+                cloudClustersFile.write("\nSubtree pathlength is %d and path 
is %s" %(maxPathLength, shortestPath))
+                #update seenNode info
+                seenNode.append(shortestPath)  
+
+    print "pruned graph"
+    print "number of nodes are: ", PRUNED_GRAPH.number_of_nodes()
+    print "number of edges are: ", PRUNED_GRAPH.number_of_edges()
+    print ("*"*80)             
+                    
+    graphTitle = "Cloud Clusters observed over somewhere during sometime"
+    #drawGraph(PRUNED_GRAPH, graphTitle, edgeWeight)
+    cloudClustersFile.close
+    
+    return PRUNED_GRAPH  
 #******************************************************************
 def findMCC (prunedGraph):
-       '''
-       Purpose:: 
-               Determines if subtree is a MCC according to Laurent et al 1998 
criteria
-
-       Input:: 
-               prunedGraph: a Networkx Graph representing the CCs 
-
-       Output:: 
-               finalMCCList: a list of list of tuples representing a MCC
-                        
-       Assumptions: 
-               frames are ordered and are equally distributed in time e.g. 
hrly satellite images
+    '''
+    Purpose:: 
+        Determines if subtree is a MCC according to Laurent et al 1998 criteria
+
+    Input:: 
+        prunedGraph: a Networkx Graph representing the CCs 
+
+    Output:: 
+        finalMCCList: a list of list of tuples representing a MCC
+             
+    Assumptions: 
+        frames are ordered and are equally distributed in time e.g. hrly 
satellite images
  
-       '''
-       MCCList = []
-       MCSList = []
-       definiteMCC = []
-       definiteMCS = []
-       eachList =[]
-       eachMCCList =[]
-       maturing = False
-       decaying = False
-       fNode = ''
-       lNode = ''
-       removeList =[]
-       imgCount = 0
-       imgTitle =''
-       
-       maxShieldNode = ''
-       orderedPath =[]
-       treeTraversalList =[]
-       definiteMCCFlag = False
-       unDirGraph = nx.Graph()
-       aSubGraph = nx.DiGraph()
-       definiteMCSFlag = False
-
-       
-       #connected_components is not available for DiGraph, so generate graph 
as undirected 
-       unDirGraph = PRUNED_GRAPH.to_undirected()
-       subGraph = nx.connected_component_subgraphs(unDirGraph)
-
-       #for each path in the subgraphs determined
-       for path in subGraph:
-               #definite is a subTree provided the duration is longer than 3 
hours
-
-               if len(path.nodes()) > MIN_MCS_DURATION:
-                       orderedPath = path.nodes()
-                       orderedPath.sort(key=lambda 
item:(len(item.split('C')[0]), item.split('C')[0]))
-                       #definiteMCS.append(orderedPath)
-
-                       #build back DiGraph for checking purposes/paper purposes
-                       aSubGraph.add_nodes_from(path.nodes())  
-                       for eachNode in path.nodes():
-                               if prunedGraph.predecessors(eachNode):
-                                       for node in 
prunedGraph.predecessors(eachNode):
-                                               
aSubGraph.add_edge(node,eachNode,weight=edgeWeight[0])
-
-                               if prunedGraph.successors(eachNode):
-                                       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:
-                       imgCount +=1
-                       #----------end build back 
---------------------------------------------
-
-                       mergeList, splitList = hasMergesOrSplits(path)  
-                       #add node behavior regarding neutral, merge, split or 
both
-                       for node in path:
-                               if node in mergeList and node in splitList:
-                                       addNodeBehaviorIdentifier(node,'B')
-                               elif node in mergeList and not node in 
splitList:
-                                       addNodeBehaviorIdentifier(node,'M')
-                               elif node in splitList and not node in 
mergeList:
-                                       addNodeBehaviorIdentifier(node,'S')
-                               else:
-                                       addNodeBehaviorIdentifier(node,'N')
-                       
-
-                       #Do the first part of checking for the MCC feature
-                       #find the path
-                       treeTraversalList = traverseTree(aSubGraph, 
orderedPath[0],[],[])
-                       #print "treeTraversalList is ", treeTraversalList
-                       #check the nodes to determine if a MCC on just the area 
criteria (consecutive nodes meeting the area and temp requirements)
-                       MCCList = checkedNodesMCC(prunedGraph, 
treeTraversalList)
-                       for aDict in MCCList:
-                               for eachNode in aDict["fullMCSMCC"]:
-                                       
addNodeMCSIdentifier(eachNode[0],eachNode[1])
-                               
-                       #do check for if MCCs overlap
-                       if MCCList:
-                               if len(MCCList) > 1:
-                                       for count in range(len(MCCList)): #for 
eachDict in MCCList:
-                                               #if there are more than two 
lists
-                                               if count >= 1:
-                                                       #and the first node in 
this list
-                                                       eachList = list(x[0] 
for x in MCCList[count]["possMCCList"])
-                                                       
eachList.sort(key=lambda nodeID:(len(nodeID.split('C')[0]), 
nodeID.split('C')[0]))
-                                                       if eachList:
-                                                               fNode = 
eachList[0]
-                                                               #get the 
lastNode in the previous possMCC list
-                                                               eachList = 
list(x[0] for x in MCCList[(count-1)]["possMCCList"])
-                                                               
eachList.sort(key=lambda nodeID:(len(nodeID.split('C')[0]), 
nodeID.split('C')[0]))
-                                                               if eachList:
-                                                                       lNode = 
eachList[-1]
-                                                                       if 
lNode in CLOUD_ELEMENT_GRAPH.predecessors(fNode):
-                                                                               
for aNode in CLOUD_ELEMENT_GRAPH.predecessors(fNode):
-                                                                               
        if aNode in eachList and aNode == lNode:
-                                                                               
                #if edge_data is equal or less than to the exisitng edge in the 
tree append one to the other
-                                                                               
                if CLOUD_ELEMENT_GRAPH.get_edge_data(aNode,fNode)['weight'] <= 
CLOUD_ELEMENT_GRAPH.get_edge_data(lNode,fNode)['weight']:
-                                                                               
                        
MCCList[count-1]["possMCCList"].extend(MCCList[count]["possMCCList"]) 
-                                                                               
                        
MCCList[count-1]["fullMCSMCC"].extend(MCCList[count]["fullMCSMCC"])
-                                                                               
                        MCCList[count-1]["durationAandB"] +=  
MCCList[count]["durationAandB"]
-                                                                               
                        MCCList[count-1]["CounterCriteriaA"] += 
MCCList[count]["CounterCriteriaA"]
-                                                                               
                        MCCList[count-1]["highestMCCnode"] = 
MCCList[count]["highestMCCnode"]
-                                                                               
                        MCCList[count-1]["frameNum"] = 
MCCList[count]["frameNum"] 
-                                                                               
                        removeList.append(count)
-                               #update the MCCList
-                               if removeList:
-                                       for i in removeList:
-                                               if (len(MCCList)-1) > i:
-                                                       del MCCList[i]
-                                                       removeList =[]
-                               
-                       #check if the nodes also meet the duration criteria and 
the shape crieria
-                       for eachDict in MCCList:
-                               #order the fullMCSMCC list, then run maximum 
extent and eccentricity criteria 
-                               if (eachDict["durationAandB"] * TRES) >= 
MINIMUM_DURATION and (eachDict["durationAandB"] * TRES) <= MAXIMUM_DURATION:
-                                       eachList = list(x[0] for x in 
eachDict["fullMCSMCC"])
-                                       eachList.sort(key=lambda 
nodeID:(len(nodeID.split('C')[0]), nodeID.split('C')[0]))
-                                       eachMCCList = list(x[0] for x in 
eachDict["possMCCList"])
-                                       eachMCCList.sort(key=lambda 
nodeID:(len(nodeID.split('C')[0]), nodeID.split('C')[0]))
-                                       
-                                       #update the nodemcsidentifer behavior
-                                       #find the first element eachMCCList in 
eachList, and ensure everything ahead of it is indicated as 'I', 
-                                       #find last element in eachMCCList in 
eachList and ensure everything after it is indicated as 'D'
-                                       #ensure that everything between is 
listed as 'M'
-                                       for eachNode in 
eachList[:(eachList.index(eachMCCList[0]))]: 
-                                               
addNodeMCSIdentifier(eachNode,'I')
-
-                                       addNodeMCSIdentifier(eachMCCList[0],'M')
-
-                                       for eachNode in 
eachList[(eachList.index(eachMCCList[-1])+1):]:
-                                               addNodeMCSIdentifier(eachNode, 
'D')
-
-                                       #update definiteMCS list
-                                       for eachNode in 
orderedPath[(orderedPath.index(eachMCCList[-1])+1):]:
-                                               addNodeMCSIdentifier(eachNode, 
'D')
-
-                                       #run maximum extent and eccentricity 
criteria
-                                       maxExtentNode, definiteMCCFlag = 
maxExtentAndEccentricity(eachList)
-                                       #print "maxExtentNode, definiteMCCFlag 
", maxExtentNode, definiteMCCFlag
-                                       if definiteMCCFlag == True:
-                                               definiteMCC.append(eachList)
-
-
-                       definiteMCS.append(orderedPath)
-                       
-                       #reset for next subGraph        
-                       aSubGraph.clear()
-                       orderedPath=[]
-                       MCCList =[]
-                       MCSList =[]
-                       definiteMCSFlag = False
-               
-       return definiteMCC, definiteMCS
+    '''
+    MCCList = []
+    MCSList = []
+    definiteMCC = []
+    definiteMCS = []
+    eachList =[]
+    eachMCCList =[]
+    maturing = False
+    decaying = False
+    fNode = ''
+    lNode = ''
+    removeList =[]
+    imgCount = 0
+    imgTitle =''
+    
+    maxShieldNode = ''
+    orderedPath =[]
+    treeTraversalList =[]
+    definiteMCCFlag = False
+    unDirGraph = nx.Graph()
+    aSubGraph = nx.DiGraph()
+    definiteMCSFlag = False
+
+    
+    #connected_components is not available for DiGraph, so generate graph as 
undirected 
+    unDirGraph = PRUNED_GRAPH.to_undirected()
+    subGraph = nx.connected_component_subgraphs(unDirGraph)
+
+    #for each path in the subgraphs determined
+    for path in subGraph:
+        #definite is a subTree provided the duration is longer than 3 hours
+
+        if len(path.nodes()) > MIN_MCS_DURATION:
+            orderedPath = path.nodes()
+            orderedPath.sort(key=lambda item:(len(item.split('C')[0]), 
item.split('C')[0]))
+            #definiteMCS.append(orderedPath)
+
+            #build back DiGraph for checking purposes/paper purposes
+            aSubGraph.add_nodes_from(path.nodes())     
+            for eachNode in path.nodes():
+                if prunedGraph.predecessors(eachNode):
+                    for node in prunedGraph.predecessors(eachNode):
+                        aSubGraph.add_edge(node,eachNode,weight=edgeWeight[0])
+
+                if prunedGraph.successors(eachNode):
+                    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:
+            imgCount +=1
+            #----------end build back 
---------------------------------------------
+
+            mergeList, splitList = hasMergesOrSplits(path)     
+            #add node behavior regarding neutral, merge, split or both
+            for node in path:
+                if node in mergeList and node in splitList:
+                    addNodeBehaviorIdentifier(node,'B')
+                elif node in mergeList and not node in splitList:
+                    addNodeBehaviorIdentifier(node,'M')
+                elif node in splitList and not node in mergeList:
+                    addNodeBehaviorIdentifier(node,'S')
+                else:
+                    addNodeBehaviorIdentifier(node,'N')
+            
+
+            #Do the first part of checking for the MCC feature
+            #find the path
+            treeTraversalList = traverseTree(aSubGraph, orderedPath[0],[],[])
+            #print "treeTraversalList is ", treeTraversalList
+            #check the nodes to determine if a MCC on just the area criteria 
(consecutive nodes meeting the area and temp requirements)
+            MCCList = checkedNodesMCC(prunedGraph, treeTraversalList)
+            for aDict in MCCList:
+                for eachNode in aDict["fullMCSMCC"]:
+                    addNodeMCSIdentifier(eachNode[0],eachNode[1])
+                
+            #do check for if MCCs overlap
+            if MCCList:
+                if len(MCCList) > 1:
+                    for count in range(len(MCCList)): #for eachDict in MCCList:
+                        #if there are more than two lists
+                        if count >= 1:
+                            #and the first node in this list
+                            eachList = list(x[0] for x in 
MCCList[count]["possMCCList"])
+                            eachList.sort(key=lambda 
nodeID:(len(nodeID.split('C')[0]), nodeID.split('C')[0]))
+                            if eachList:
+                                fNode = eachList[0]
+                                #get the lastNode in the previous possMCC list
+                                eachList = list(x[0] for x in 
MCCList[(count-1)]["possMCCList"])
+                                eachList.sort(key=lambda 
nodeID:(len(nodeID.split('C')[0]), nodeID.split('C')[0]))
+                                if eachList:
+                                    lNode = eachList[-1]
+                                    if lNode in 
CLOUD_ELEMENT_GRAPH.predecessors(fNode):
+                                        for aNode in 
CLOUD_ELEMENT_GRAPH.predecessors(fNode):
+                                            if aNode in eachList and aNode == 
lNode:
+                                                #if edge_data is equal or less 
than to the exisitng edge in the tree append one to the other
+                                                if 
CLOUD_ELEMENT_GRAPH.get_edge_data(aNode,fNode)['weight'] <= 
CLOUD_ELEMENT_GRAPH.get_edge_data(lNode,fNode)['weight']:
+                                                    
MCCList[count-1]["possMCCList"].extend(MCCList[count]["possMCCList"]) 
+                                                    
MCCList[count-1]["fullMCSMCC"].extend(MCCList[count]["fullMCSMCC"])
+                                                    
MCCList[count-1]["durationAandB"] +=  MCCList[count]["durationAandB"]
+                                                    
MCCList[count-1]["CounterCriteriaA"] += MCCList[count]["CounterCriteriaA"]
+                                                    
MCCList[count-1]["highestMCCnode"] = MCCList[count]["highestMCCnode"]
+                                                    
MCCList[count-1]["frameNum"] = MCCList[count]["frameNum"] 
+                                                    removeList.append(count)
+                #update the MCCList
+                if removeList:
+                    for i in removeList:
+                        if (len(MCCList)-1) > i:
+                            del MCCList[i]
+                            removeList =[]
+                
+            #check if the nodes also meet the duration criteria and the shape 
crieria
+            for eachDict in MCCList:
+                #order the fullMCSMCC list, then run maximum extent and 
eccentricity criteria 
+                if (eachDict["durationAandB"] * TRES) >= MINIMUM_DURATION and 
(eachDict["durationAandB"] * TRES) <= MAXIMUM_DURATION:
+                    eachList = list(x[0] for x in eachDict["fullMCSMCC"])
+                    eachList.sort(key=lambda 
nodeID:(len(nodeID.split('C')[0]), nodeID.split('C')[0]))
+                    eachMCCList = list(x[0] for x in eachDict["possMCCList"])
+                    eachMCCList.sort(key=lambda 
nodeID:(len(nodeID.split('C')[0]), nodeID.split('C')[0]))
+                    
+                    #update the nodemcsidentifer behavior
+                    #find the first element eachMCCList in eachList, and 
ensure everything ahead of it is indicated as 'I', 
+                    #find last element in eachMCCList in eachList and ensure 
everything after it is indicated as 'D'
+                    #ensure that everything between is listed as 'M'
+                    for eachNode in 
eachList[:(eachList.index(eachMCCList[0]))]: 
+                        addNodeMCSIdentifier(eachNode,'I')
+
+                    addNodeMCSIdentifier(eachMCCList[0],'M')
+
+                    for eachNode in 
eachList[(eachList.index(eachMCCList[-1])+1):]:
+                        addNodeMCSIdentifier(eachNode, 'D')
+
+                    #update definiteMCS list
+                    for eachNode in 
orderedPath[(orderedPath.index(eachMCCList[-1])+1):]:
+                        addNodeMCSIdentifier(eachNode, 'D')
+
+                    #run maximum extent and eccentricity criteria
+                    maxExtentNode, definiteMCCFlag = 
maxExtentAndEccentricity(eachList)
+                    #print "maxExtentNode, definiteMCCFlag ", maxExtentNode, 
definiteMCCFlag
+                    if def

<TRUNCATED>

Reply via email to