climate-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From lewi...@apache.org
Subject [2/3] climate git commit: CLIMATE-922 Optimize traversing in mccSearch
Date Mon, 11 Sep 2017 15:32:18 GMT
http://git-wip-us.apache.org/repos/asf/climate/blob/cf4fb57f/mccsearch/code/mccSearch.py
----------------------------------------------------------------------
diff --cc mccsearch/code/mccSearch.py
index c63b243,3a84dd5..0d3b669
--- a/mccsearch/code/mccSearch.py
+++ b/mccsearch/code/mccSearch.py
@@@ -44,2626 -44,2114 +44,2891 @@@ import ocw.plotter as plotte
  
  #----------------------- GLOBAL VARIABLES --------------------------
  # --------------------- User defined variables ---------------------
 -#FYI the lat lon values are not necessarily inclusive of the points given. These are the limits
 -#the first point closest the the value (for the min) from the MERG data is used, etc.
 -LATMIN = '5.0' #min latitude; -ve values in the SH e.g. 5S = -5
 -LATMAX = '19.0' #max latitude; -ve values in the SH e.g. 5S = -5 20.0
 -LONMIN = '-5.0' #min longitude; -ve values in the WH e.g. 59.8W = -59.8 -30
 -LONMAX = '5.0' #min longitude; -ve values in the WH e.g. 59.8W = -59.8  30
 -XRES = 4.0				#x direction spatial resolution in km
 -YRES = 4.0				#y direction spatial resolution in km
 -TRES = 1 				#temporal resolution in hrs
 -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 
 -#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
 -CONVECTIVE_FRACTION = 0.90 #the min temp/max temp that would be expected in a CE.. this is highly conservative (only a 10K difference)
 -MIN_MCS_DURATION = 3	#minimum time for a MCS to exist
 -AREA_MIN = 2400.0		#minimum area for CE criteria in km^2 according to Vila et al. (2008) is 2400
 -MIN_OVERLAP= 10000.00   #km^2  from Williams and Houze 1987, indir ref in Arnaud et al 1992
 +# FYI the lat lon values are not necessarily inclusive of the points given. These are the limits
 +# the first point closest the the value (for the min) from the MERG data
 +# is used, etc.
 +LATMIN = '5.0'  # min latitude; -ve values in the SH e.g. 5S = -5
 +LATMAX = '19.0'  # max latitude; -ve values in the SH e.g. 5S = -5 20.0
 +LONMIN = '-5.0'  # min longitude; -ve values in the WH e.g. 59.8W = -59.8 -30
 +LONMAX = '5.0'  # min longitude; -ve values in the WH e.g. 59.8W = -59.8  30
 +XRES = 4.0  # x direction spatial resolution in km
 +YRES = 4.0  # y direction spatial resolution in km
 +TRES = 1  # temporal resolution in hrs
 +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
 +# the matrix for determining the pattern for the contiguous boxes and must
 +STRUCTURING_ELEMENT = [[0, 1, 0], [1, 1, 1], [0, 1, 0]]
 +# have same rank of the matrix it is being compared against
 +# criteria for determining cloud elements and edges
 +# warmest temp to allow (-30C to -55C according to Morel and Sensi 2002)
 +T_BB_MAX = 243
 +T_BB_MIN = 218  # cooler temp for the center of the system
 +# the min temp/max temp that would be expected in a CE.. this is highly
 +# conservative (only a 10K difference)
 +CONVECTIVE_FRACTION = 0.90
 +MIN_MCS_DURATION = 3  # minimum time for a MCS to exist
 +# minimum area for CE criteria in km^2 according to Vila et al. (2008) is 2400
 +AREA_MIN = 2400.0
 +# km^2  from Williams and Houze 1987, indir ref in Arnaud et al 1992
 +MIN_OVERLAP = 10000.00
  
  #---the MCC criteria
 -ECCENTRICITY_THRESHOLD_MAX = 1.0  #tending to 1 is a circle e.g. hurricane, 
 -ECCENTRICITY_THRESHOLD_MIN = 0.50 #tending to 0 is a linear e.g. squall line
 -OUTER_CLOUD_SHIELD_AREA = 80000.0 #km^2
 -INNER_CLOUD_SHIELD_AREA = 30000.0 #km^2
 -OUTER_CLOUD_SHIELD_TEMPERATURE = 233 #in K
 -INNER_CLOUD_SHIELD_TEMPERATURE = 213 #in K
 -MINIMUM_DURATION = 6 #min number of frames the MCC must exist for (assuming hrly frames, MCCs is 6hrs)
 -MAXIMUM_DURATION = 24#max number of framce the MCC can last for 
 +ECCENTRICITY_THRESHOLD_MAX = 1.0  # tending to 1 is a circle e.g. hurricane,
 +ECCENTRICITY_THRESHOLD_MIN = 0.50  # tending to 0 is a linear e.g. squall line
 +OUTER_CLOUD_SHIELD_AREA = 80000.0  # km^2
 +INNER_CLOUD_SHIELD_AREA = 30000.0  # km^2
 +OUTER_CLOUD_SHIELD_TEMPERATURE = 233  # in K
 +INNER_CLOUD_SHIELD_TEMPERATURE = 213  # in K
 +# min number of frames the MCC must exist for (assuming hrly frames, MCCs
 +# is 6hrs)
 +MINIMUM_DURATION = 6
 +MAXIMUM_DURATION = 24  # max number of framce the MCC can last for
  #------------------- End user defined Variables -------------------
 -edgeWeight = [1,2,3] #weights for the graph edges
 -#graph object fo the CEs meeting the criteria
 +edgeWeight = [1, 2, 3]  # weights for the graph edges
 +# graph object fo the CEs meeting the criteria
  CLOUD_ELEMENT_GRAPH = nx.DiGraph()
 -#graph meeting the CC criteria
 +# graph meeting the CC criteria
  PRUNED_GRAPH = nx.DiGraph()
  #------------------------ End GLOBAL VARS -------------------------
 +
  #************************ 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
 +
 +
 +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 is 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 BaseException:
 +            print(("bad file! %s" % 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	
 +
 +
 +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 range(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 range(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 as e:
 +                print(("Error is %s" % 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 BaseException:
 +                        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: %s" % CLOUD_ELEMENT_GRAPH.number_of_nodes()))
 +    print(("number of edges are: %s" % 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
 -
 -	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
 -   
 -	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)
 -		
 -		fileHr = fileHr1 if int(fileHr1) % temporalRes == 0 else (int(fileHr1) / temporalRes) * temporalRes
 -		fileHr = '0' + str(fileHr) if fileHr < 10 else str(fileHr)
 -
 -		TRMMfileName = TRMMdirName+"/3B42."+ str(fileDate) + "."+ 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):
 -			if value > 0:
 -				lat_index, lon_index = index
 -				currTimeValue = 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
 -#******************************************************************	
 +    '''
 +    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
 +
 +    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
 +
 +    '''
 +    allCEnodesTRMMdata = []
 +    TRMMdataDict = {}
 +    precipTotal = 0.0
 +
 +    os.chdir((MAINDIRECTORY + '/MERGnetcdfCEs/'))
 +    imgFilename = ''
 +    temporalRes = 3  # 3 hours for TRMM
 +
 +    # sort files
 +    files = list(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)
++        fileHr = fileHr1 if int(fileHr1) % temporalRes == 0 else (int(fileHr1) / temporalRes) * temporalRes
++        fileHr = '0' + str(fileHr) if fileHr < 10 else str(fileHr)
 +
-         TRMMfileName = TRMMdirName + "/3B42." + \
-             str(fileDate) + "." + str(fileHr) + ".7A.nc"
++        TRMMfileName = TRMMdirName + "/3B42." + str(fileDate) + "." + 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 BaseException:
 +            minCEprecipRate = 0.0
 +
 +        try:
 +            maxCEprecipRate = np.max(
 +                finalCETRMMvalues[np.nonzero(finalCETRMMvalues)])
 +        except BaseException:
 +            maxCEprecipRate = 0.0
 +
 +        # add info to CLOUDELEMENTSGRAPH
 +        # TODO try block
 +        for eachdict in CLOUD_ELEMENT_GRAPH.nodes(CEuniqueID):
 +            if eachdict[1]['uniqueID'] == CEuniqueID:
 +                if 'cloudElementPrecipTotal' not in list(eachdict[1].keys()):
 +                    eachdict[1]['cloudElementPrecipTotal'] = precipTotal
 +                if 'cloudElementLatLonTRMM' not in list(eachdict[1].keys()):
 +                    eachdict[1]['cloudElementLatLonTRMM'] = finalCETRMMvalues
 +                if 'TRMMArea' not in list(eachdict[1].keys()):
 +                    eachdict[1]['TRMMArea'] = TRMMArea
 +                if 'CETRMMmin' not in list(eachdict[1].keys()):
 +                    eachdict[1]['CETRMMmin'] = minCEprecipRate
 +                if 'CETRMMmax' not in list(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 range(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: %s", PRUNED_GRAPH.number_of_nodes()))
-     print(("number of edges are: %s", PRUNED_GRAPH.number_of_edges()))
-     print(("*" * 80))
++    print("number of nodes are: %s", PRUNED_GRAPH.number_of_nodes())
++    print("number of edges are: %s", 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
 - 
 -	'''
 -	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],[], set(), [])
 -			#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 == lNod

<TRUNCATED>

Mime
View raw message