#
#  Licensed to the Apache Software Foundation (ASF) under one or more
#  contributor license agreements.  See the NOTICE file distributed with
#  this work for additional information regarding copyright ownership.
#  The ASF licenses this file to You under the Apache License, Version 2.0
#  (the "License"); you may not use this file except in compliance with
#  the License.  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#

from datetime import timedelta, datetime
import glob
import itertools
from netCDF4 import Dataset, date2num
import numpy as np
import numpy.ma as ma
import os
import pickle
import re
from scipy import ndimage
import string
import subprocess
import sys
import time

import networkx as nx
import matplotlib.dates as mdates
from matplotlib.dates import DateFormatter, HourLocator
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import cm as cmbm
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors
from matplotlib.ticker import FuncFormatter, FormatStrFormatter

import ocw.plotter as plotter

#----------------------- 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
# 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
# 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
CLOUD_ELEMENT_GRAPH = nx.DiGraph()
# 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 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 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 = 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)

        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):
            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 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)

    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 node not in splitList:
                    addNodeBehaviorIdentifier(node, 'M')
                elif node in splitList and node not 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 %s" % 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 eachDict in MCCList:
                    for count in range(len(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)
                    if definiteMCCFlag:
                        definiteMCC.append(eachList)

            definiteMCS.append(orderedPath)

            # reset for next subGraph
            aSubGraph.clear()
            orderedPath = []
            MCCList = []
            MCSList = []
            definiteMCSFlag = False

    return definiteMCC, definiteMCS
#******************************************************************


def traverseTree(subGraph, node, stack, checkedNodes=None):
    '''
    Purpose::
            To traverse a tree using a modified depth-first iterative deepening (DFID) search algorithm

    Input::
            subGraph: a Networkx DiGraph representing a CC
                    lengthOfsubGraph: an integer representing the length of the subgraph
                    node: a string representing the node currently being checked
                    stack: a list of strings representing a list of nodes in a stack functionality
                                    i.e. Last-In-First-Out (LIFO) for sorting the information from each visited node
                    checkedNodes: a list of strings representing the list of the nodes in the traversal

    Output::
            checkedNodes: a list of strings representing the list of the nodes in the traversal

    Assumptions:
            frames are ordered and are equally distributed in time e.g. hrly satellite images
    '''
    if len(checkedNodes) == len(subGraph):
        return checkedNodes

    if not checkedNodes:
        stack = []
        checkedNodes.append(node)

    # check one level infront first...if something does exisit, stick it at
    # the front of the stack
    upOneLevel = subGraph.predecessors(node)
    downOneLevel = subGraph.successors(node)
    for parent in upOneLevel:
        if parent not in checkedNodes and parent not in stack:
            for child in downOneLevel:
                if child not in checkedNodes and child not in stack:
                    stack.insert(0, child)

            stack.insert(0, parent)
    for child in downOneLevel:
        if child not in checkedNodes and child not in stack:
            if len(subGraph.predecessors(child)) > 1 or node in checkedNodes:
                stack.insert(0, child)
            else:
                stack.append(child)

    for eachNode in stack:
        if eachNode not in checkedNodes:
            checkedNodes.append(eachNode)
            return traverseTree(subGraph, eachNode, stack, checkedNodes)

    return checkedNodes


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 node not in splitList:
                    addNodeBehaviorIdentifier(node, 'M')
                elif node in splitList and node not 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 eachDict in MCCList:
                    for count in range(len(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:
                        definiteMCC.append(eachList)

            definiteMCS.append(orderedPath)

            # reset for next subGraph
            aSubGraph.clear()
            orderedPath = []
            MCCList = []
            MCSList = []
            definiteMCSFlag = False

    return definiteMCC, definiteMCS
#******************************************************************


def traverseTree(subGraph, node, stack, bookkeeper_stack, checkedNodes=None):
    '''
    Purpose::
            To traverse a tree using a modified depth-first iterative deepening (DFID) search algorithm

    Input::
            subGraph: a Networkx DiGraph representing a CC
                    lengthOfsubGraph: an integer representing the length of the subgraph
                    node: a string representing the node currently being checked
                    stack: a list of strings representing a list of nodes in a stack functionality
                                    i.e. Last-In-First-Out (LIFO) for sorting the information from each visited node
                    checkedNodes: a list of strings representing the list of the nodes in the traversal

    Output::
            checkedNodes: a list of strings representing the list of the nodes in the traversal

    Assumptions:
            frames are ordered and are equally distributed in time e.g. hrly satellite images

    '''
    if len(checkedNodes) == len(subGraph):
        return checkedNodes

    if not checkedNodes:
        stack = []
        bookkeeper_stack = set()
        checkedNodes.append(node)

    # check one level infront first...if something does exisit, stick it at
    # the front of the stack
    upOneLevel = subGraph.predecessors(node)
    downOneLevel = subGraph.successors(node)
    for parent in upOneLevel:
        if parent not in checkedNodes and parent not in bookkeeper_stack:
            for child in downOneLevel:
                if child not in checkedNodes and child not in bookkeeper_stack:
                    stack.insert(0, child)
                    bookkeeper_stack.add(child)
            stack.insert(0, parent)
            bookkeeper_stack.add(parent)

    for child in downOneLevel:
        if child not in checkedNodes and child not in bookkeeper_stack:
            stack.insert(0, child)
            bookkeeper_stack.add(child)

    for eachNode in stack:
        if eachNode not in checkedNodes:
            checkedNodes.append(eachNode)
            return traverseTree(
                subGraph,
                eachNode,
                stack,
                bookkeeper_stack,
                checkedNodes)

    return checkedNodes
#******************************************************************


def checkedNodesMCC(prunedGraph, nodeList):
    '''
    Purpose ::
            Determine if this path is (or is part of) a MCC and provides
            preliminary information regarding the stages of the feature

    Input::
            prunedGraph: a Networkx Graph representing all the cloud clusters
            nodeList: list of strings (CE ID) from the traversal

    Output::
            potentialMCCList: list of dictionaries representing all possible MCC within the path
                    dictionary = {"possMCCList":[(node,'I')], "fullMCSMCC":[(node,'I')], "CounterCriteriaA": CounterCriteriaA, "durationAandB": durationAandB}
    '''

    CounterCriteriaAFlag = False
    CounterCriteriaBFlag = False
    INITIATIONFLAG = False
    MATURITYFLAG = False
    DECAYFLAG = False
    thisdict = {}  # will have the same items as the cloudElementDict
    cloudElementAreaB = 0.0
    cloudElementAreaA = 0.0
    epsilon = 0.0
    frameNum = 0
    oldNode = ''
    potentialMCCList = []
    durationAandB = 0

    # check for if the list contains only one string/node
    if isinstance(nodeList, str):
        oldNode = nodeList
        nodeList = []
        nodeList.append(oldNode)

    for node in nodeList:
        thisdict = thisDict(node)
        CounterCriteriaAFlag = False
        CounterCriteriaBFlag = False
        existingFrameFlag = False

        if thisdict['cloudElementArea'] >= OUTER_CLOUD_SHIELD_AREA:
            CounterCriteriaAFlag = True
            INITIATIONFLAG = True
            MATURITYFLAG = False

            # check if criteriaA is met
            cloudElementAreaA, criteriaA = checkCriteria(
                thisdict['cloudElementLatLon'], OUTER_CLOUD_SHIELD_TEMPERATURE)
            # TODO: calcuate the eccentricity at this point and read over????or
            # create a new field in the dict

            if cloudElementAreaA >= OUTER_CLOUD_SHIELD_AREA:
                # check if criteriaB is met
                cloudElementAreaB, criteriaB = checkCriteria(
                    thisdict['cloudElementLatLon'], INNER_CLOUD_SHIELD_TEMPERATURE)

                # if Criteria A and B have been met, then the MCC is initiated,
                # i.e. store node as potentialMCC
                if cloudElementAreaB >= INNER_CLOUD_SHIELD_AREA:
                    # TODO: add another field to the dictionary for the
                    # OUTER_AREA_SHIELD area
                    CounterCriteriaBFlag = True
                    # append this information on to the dictionary
                    addInfothisDict(node, cloudElementAreaB, criteriaB)
                    INITIATIONFLAG = False
                    MATURITYFLAG = True
                    stage = 'M'
                    potentialMCCList = updateMCCList(
                        prunedGraph,
                        potentialMCCList,
                        node,
                        stage,
                        CounterCriteriaAFlag,
                        CounterCriteriaBFlag)
                else:
                    # criteria B failed
                    CounterCriteriaBFlag = False
                    if INITIATIONFLAG:
                        stage = 'I'
                        potentialMCCList = updateMCCList(
                            prunedGraph,
                            potentialMCCList,
                            node,
                            stage,
                            CounterCriteriaAFlag,
                            CounterCriteriaBFlag)

                    elif (INITIATIONFLAG == False and MATURITYFLAG == True) or DECAYFLAG == True:
                        DECAYFLAG = True
                        MATURITYFLAG = False
                        stage = 'D'
                        potentialMCCList = updateMCCList(
                            prunedGraph,
                            potentialMCCList,
                            node,
                            stage,
                            CounterCriteriaAFlag,
                            CounterCriteriaBFlag)
            else:
                # criteria A failed
                CounterCriteriaAFlag = False
                CounterCriteriaBFlag = False
                # add as a CE before or after the main feature
                if INITIATIONFLAG or (
                        INITIATIONFLAG == False and MATURITYFLAG == True):
                    stage = "I"
                    potentialMCCList = updateMCCList(
                        prunedGraph,
                        potentialMCCList,
                        node,
                        stage,
                        CounterCriteriaAFlag,
                        CounterCriteriaBFlag)
                elif (INITIATIONFLAG == False and MATURITYFLAG == False) or DECAYFLAG == True:
                    stage = "D"
                    DECAYFLAG = True
                    potentialMCCList = updateMCCList(
                        prunedGraph,
                        potentialMCCList,
                        node,
                        stage,
                        CounterCriteriaAFlag,
                        CounterCriteriaBFlag)
                elif (INITIATIONFLAG == False and MATURITYFLAG == False and DECAYFLAG == False):
                    stage = "I"
                    potentialMCCList = updateMCCList(
                        prunedGraph,
                        potentialMCCList,
                        node,
                        stage,
                        CounterCriteriaAFlag,
                        CounterCriteriaBFlag)

        else:
            # criteria A failed
            CounterCriteriaAFlag = False
            CounterCriteriaBFlag = False
            # add as a CE before or after the main feature
            if INITIATIONFLAG or (
                    INITIATIONFLAG == False and MATURITYFLAG == True):
                stage = "I"
                potentialMCCList = updateMCCList(
                    prunedGraph,
                    potentialMCCList,
                    node,
                    stage,
                    CounterCriteriaAFlag,
                    CounterCriteriaBFlag)
            elif (INITIATIONFLAG == False and MATURITYFLAG == False) or DECAYFLAG == True:
                stage = "D"
                DECAYFLAG = True
                potentialMCCList = updateMCCList(
                    prunedGraph,
                    potentialMCCList,
                    node,
                    stage,
                    CounterCriteriaAFlag,
                    CounterCriteriaBFlag)
            elif (INITIATIONFLAG == False and MATURITYFLAG == False and DECAYFLAG == False):
                stage = "I"
                potentialMCCList = updateMCCList(
                    prunedGraph,
                    potentialMCCList,
                    node,
                    stage,
                    CounterCriteriaAFlag,
                    CounterCriteriaBFlag)

    return potentialMCCList
#******************************************************************


def updateMCCList(
        prunedGraph,
        potentialMCCList,
        node,
        stage,
        CounterCriteriaAFlag,
        CounterCriteriaBFlag):
    '''
    Purpose::
            Utility function to determine if a path is (or is part of) a MCC and provides
                       preliminary information regarding the stages of the feature

    Input::
            prunedGraph: a Networkx Graph representing all the cloud clusters
            potentialMCCList: a list of dictionaries representing the possible MCCs within a path
            node: a string representing the cloud element currently being assessed
            CounterCriteriaAFlag: a boolean value indicating whether the node meets the MCC criteria A according to Laurent et al
            CounterCriteriaBFlag: a boolean value indicating whether the node meets the MCC criteria B according to Laurent et al

    Output::
            potentialMCCList: list of dictionaries representing all possible MCC within the path
                     dictionary = {"possMCCList":[(node,'I')], "fullMCSMCC":[(node,'I')], "CounterCriteriaA": CounterCriteriaA, "durationAandB": durationAandB}

    '''
    existingFrameFlag = False
    existingMCSFrameFlag = False
    predecessorsFlag = False
    predecessorsMCSFlag = False
    successorsFlag = False
    successorsMCSFlag = False
    frameNum = 0

    frameNum = int((node.split('CE')[0]).split('F')[1])
    if potentialMCCList == []:
        # list empty
        stage = 'I'
        if CounterCriteriaAFlag and CounterCriteriaBFlag:
            potentialMCCList.append(
                {
                    "possMCCList": [
                        (node,
                         stage)],
                    "fullMCSMCC": [
                        (node,
                         stage)],
                    "CounterCriteriaA": 1,
                    "durationAandB": 1,
                    "highestMCCnode": node,
                    "frameNum": frameNum})
        elif CounterCriteriaAFlag and CounterCriteriaBFlag == False:
            potentialMCCList.append({"possMCCList": [],
                                     "fullMCSMCC": [(node,
                                                     stage)],
                                     "CounterCriteriaA": 1,
                                     "durationAandB": 0,
                                     "highestMCCnode": "",
                                     "frameNum": 0})
        elif CounterCriteriaAFlag == False and CounterCriteriaBFlag == False:
            potentialMCCList.append({"possMCCList": [],
                                     "fullMCSMCC": [(node,
                                                     stage)],
                                     "CounterCriteriaA": 0,
                                     "durationAandB": 0,
                                     "highestMCCnode": "",
                                     "frameNum": 0})

    else:
        # list not empty
        predecessorsFlag, index = isThereALink(
            prunedGraph, 1, node, potentialMCCList, 1)

        if predecessorsFlag:

            for eachNode in potentialMCCList[index]["possMCCList"]:
                if int((eachNode[0].split('CE')[0]).split('F')[1]) == frameNum:
                    existingFrameFlag = True

            # this MUST come after the check for the existing frame
            if CounterCriteriaAFlag and CounterCriteriaBFlag:
                stage = 'M'
                potentialMCCList[index]["possMCCList"].append((node, stage))
                potentialMCCList[index]["fullMCSMCC"].append((node, stage))

            if not existingFrameFlag:
                if CounterCriteriaAFlag and CounterCriteriaBFlag:
                    stage = 'M'
                    potentialMCCList[index]["CounterCriteriaA"] += 1
                    potentialMCCList[index]["durationAandB"] += 1
                    if frameNum > potentialMCCList[index]["frameNum"]:
                        potentialMCCList[index]["frameNum"] = frameNum
                        potentialMCCList[index]["highestMCCnode"] = node
                    return potentialMCCList

                # if this frameNum doesn't exist and this frameNum is less than
                # the MCC node max frame Num (including 0), then append to
                # fullMCSMCC list
                if frameNum > potentialMCCList[index]["frameNum"] or potentialMCCList[index]["frameNum"] == 0:
                    stage = 'I'
                    if CounterCriteriaAFlag and CounterCriteriaBFlag == False:
                        potentialMCCList.append(
                            {
                                "possMCCList": [],
                                "fullMCSMCC": [
                                    (node,
                                     stage)],
                                "CounterCriteriaA": 1,
                                "durationAandB": 0,
                                "highestMCCnode": "",
                                "frameNum": 0})
                        return potentialMCCList
                    elif CounterCriteriaAFlag == False and CounterCriteriaBFlag == False:
                        potentialMCCList.append(
                            {
                                "possMCCList": [],
                                "fullMCSMCC": [
                                    (node,
                                     stage)],
                                "CounterCriteriaA": 0,
                                "durationAandB": 0,
                                "highestMCCnode": "",
                                "frameNum": 0})
                        return potentialMCCList

            # if predecessor and this frame number already exist in the MCC
            # list, add the current node to the fullMCSMCC list
            if existingFrameFlag:
                if CounterCriteriaAFlag and CounterCriteriaBFlag == False:
                    potentialMCCList[index]["fullMCSMCC"].append((node, stage))
                    potentialMCCList[index]["CounterCriteriaA"] += 1
                    return potentialMCCList
                if not CounterCriteriaAFlag:
                    potentialMCCList[index]["fullMCSMCC"].append((node, stage))
                    return potentialMCCList

        if not predecessorsFlag:
            successorsFlag, index = isThereALink(
                prunedGraph, 2, node, potentialMCCList, 2)

            if successorsFlag:
                for eachNode in potentialMCCList[index]["possMCCList"]:
                    if int(
                            (eachNode[0].split('CE')[0]).split('F')[1]) == frameNum:
                        existingFrameFlag = True

                if CounterCriteriaAFlag and CounterCriteriaBFlag:
                    stage = 'M'
                    potentialMCCList[index]["possMCCList"].append(
                        (node, stage))
                    potentialMCCList[index]["fullMCSMCC"].append((node, stage))

                    if frameNum > potentialMCCList[index]["frameNum"] or potentialMCCList[index]["frameNum"] == 0:
                        potentialMCCList[index]["frameNum"] = frameNum
                        potentialMCCList[index]["highestMCCnode"] = node
                    return potentialMCCList

                if not existingFrameFlag:
                    if stage == 'M':
                        stage = 'D'
                    if CounterCriteriaAFlag and CounterCriteriaBFlag:
                        potentialMCCList[index]["CounterCriteriaA"] += 1
                        potentialMCCList[index]["durationAandB"] += 1
                    elif CounterCriteriaAFlag:
                        potentialMCCList[index]["CounterCriteriaA"] += 1
                    elif CounterCriteriaAFlag == False:
                        potentialMCCList[index]["fullMCSMCC"].append(
                            (node, stage))
                        return potentialMCCList
                        # if predecessor and this frame number already exist in
                        # the MCC list, add the current node to the fullMCSMCC
                        # list
                else:
                    if CounterCriteriaAFlag and CounterCriteriaBFlag == False:
                        potentialMCCList[index]["fullMCSMCC"].append(
                            (node, stage))
                        potentialMCCList[index]["CounterCriteriaA"] += 1
                        return potentialMCCList
                    if not CounterCriteriaAFlag:
                        potentialMCCList[index]["fullMCSMCC"].append(
                            (node, stage))
                        return potentialMCCList

        # if this node isn't connected to exisiting MCCs check if it is
        # connected to exisiting MCSs ...
        if predecessorsFlag == False and successorsFlag == False:
            stage = 'I'
            predecessorsMCSFlag, index = isThereALink(
                prunedGraph, 1, node, potentialMCCList, 2)
            if predecessorsMCSFlag:
                if CounterCriteriaAFlag and CounterCriteriaBFlag:
                    potentialMCCList[index]["possMCCList"].append((node, 'M'))
                    potentialMCCList[index]["fullMCSMCC"].append((node, 'M'))
                    potentialMCCList[index]["durationAandB"] += 1
                    if frameNum > potentialMCCList[index]["frameNum"]:
                        potentialMCCList[index]["frameNum"] = frameNum
                        potentialMCCList[index]["highestMCCnode"] = node
                    return potentialMCCList

                if potentialMCCList[index]["frameNum"] == 0 or frameNum <= potentialMCCList[index]["frameNum"]:
                    if CounterCriteriaAFlag and CounterCriteriaBFlag == False:
                        potentialMCCList[index]["fullMCSMCC"].append(
                            (node, stage))
                        potentialMCCList[index]["CounterCriteriaA"] += 1
                        return potentialMCCList
                    elif CounterCriteriaAFlag == False:
                        potentialMCCList[index]["fullMCSMCC"].append(
                            (node, stage))
                        return potentialMCCList
            else:
                successorsMCSFlag, index = isThereALink(
                    prunedGraph, 2, node, potentialMCCList, 2)
                if successorsMCSFlag:
                    if CounterCriteriaAFlag and CounterCriteriaBFlag:
                        potentialMCCList[index]["possMCCList"].append(
                            (node, 'M'))
                        potentialMCCList[index]["fullMCSMCC"].append(
                            (node, 'M'))

                        potentialMCCList[index]["durationAandB"] += 1
                        if frameNum > potentialMCCList[index]["frameNum"]:
                            potentialMCCList[index]["frameNum"] = frameNum
                            potentialMCCList[index]["highestMCCnode"] = node
                        return potentialMCCList

                    if potentialMCCList[index]["frameNum"] == 0 or frameNum <= potentialMCCList[index]["frameNum"]:
                        if CounterCriteriaAFlag and CounterCriteriaBFlag == False:
                            potentialMCCList[index]["fullMCSMCC"].append(
                                (node, stage))
                            potentialMCCList[index]["CounterCriteriaA"] += 1
                            return potentialMCCList
                        elif CounterCriteriaAFlag == False:
                            potentialMCCList[index]["fullMCSMCC"].append(
                                (node, stage))
                            return potentialMCCList

            # if this node isn't connected to existing MCCs or MCSs, create a
            # new one ...
            if predecessorsFlag == False and predecessorsMCSFlag == False and successorsFlag == False and successorsMCSFlag == False:
                if CounterCriteriaAFlag and CounterCriteriaBFlag:
                    potentialMCCList.append(
                        {
                            "possMCCList": [
                                (node,
                                 stage)],
                            "fullMCSMCC": [
                                (node,
                                 stage)],
                            "CounterCriteriaA": 1,
                            "durationAandB": 1,
                            "highestMCCnode": node,
                            "frameNum": frameNum})
                elif CounterCriteriaAFlag and CounterCriteriaBFlag == False:
                    potentialMCCList.append(
                        {
                            "possMCCList": [],
                            "fullMCSMCC": [
                                (node,
                                 stage)],
                            "CounterCriteriaA": 1,
                            "durationAandB": 0,
                            "highestMCCnode": "",
                            "frameNum": 0})
                elif CounterCriteriaAFlag == False and CounterCriteriaBFlag == False:
                    potentialMCCList.append(
                        {
                            "possMCCList": [],
                            "fullMCSMCC": [
                                (node,
                                 stage)],
                            "CounterCriteriaA": 0,
                            "durationAandB": 0,
                            "highestMCCnode": "",
                            "frameNum": 0})

    return potentialMCCList
#******************************************************************


def isThereALink(prunedGraph, upOrDown, node, potentialMCCList, whichList):
    '''
    Purpose::
            Utility script for updateMCCList mostly because there is no Pythonic way to break out of nested loops

    Input::
            prunedGraph:a Networkx Graph representing all the cloud clusters
            upOrDown: an integer representing 1- to do predecesor check and 2 - to do successor checkedNodesMCC
            node: a string representing the cloud element currently being assessed
            potentialMCCList: a list of dictionaries representing the possible MCCs within a path
            whichList: an integer representing which list ot check in the dictionary; 1- possMCCList, 2- fullMCSMCC

    Output::
            thisFlag: a boolean representing whether the list passed has in the parent or child of the node
            index: an integer representing the location in the potentialMCCList where thisFlag occurs

    '''
    thisFlag = False
    index = -1
    checkList = ""
    if whichList == 1:
        checkList = "possMCCList"
    elif whichList == 2:
        checkList = "fullMCSMCC"

    # check parents
    if upOrDown == 1:
        for aNode in prunedGraph.predecessors(node):
            # reset the index counter for this node search through
            # potentialMCCList
            index = -1
            for MCCDict in potentialMCCList:
                index += 1
                if aNode in list(x[0] for x in MCCDict[checkList]):
                    thisFlag = True
                    # get out of looping so as to avoid the flag being written
                    # over when another node in the predecesor list is checked
                    return thisFlag, index

    # check children
    if upOrDown == 2:
        for aNode in prunedGraph.successors(node):
            # reset the index counter for this node search through
            # potentialMCCList
            index = -1
            for MCCDict in potentialMCCList:
                index += 1

                if aNode in list(x[0] for x in MCCDict[checkList]):
                    thisFlag = True
                    return thisFlag, index

    return thisFlag, index
#******************************************************************


def maxExtentAndEccentricity(eachList):
    '''
    Purpose::
            Perform the final check for MCC based on maximum extent and eccentricity criteria

    Input::
            eachList: a list of strings  representing the node of the possible MCCs within a path

    Output::
            maxShieldNode: a string representing the node with the maximum maxShieldNode
            definiteMCCFlag: a boolean indicating that the MCC has met all requirements

    '''
    maxShieldNode = ''
    maxShieldArea = 0.0
    maxShieldEccentricity = 0.0
    definiteMCCFlag = False

    if eachList:
        for eachNode in eachList:
            if (thisDict(eachNode)['nodeMCSIdentifier'] == 'M' or thisDict(eachNode)[
                    'nodeMCSIdentifier'] == 'D') and thisDict(eachNode)['cloudElementArea'] > maxShieldArea:
                maxShieldNode = eachNode
                maxShieldArea = thisDict(eachNode)['cloudElementArea']

        maxShieldEccentricity = thisDict(maxShieldNode)[
            'cloudElementEccentricity']
        if thisDict(maxShieldNode)['cloudElementEccentricity'] >= ECCENTRICITY_THRESHOLD_MIN and thisDict(
                maxShieldNode)['cloudElementEccentricity'] <= ECCENTRICITY_THRESHOLD_MAX:
            # criteria met
            definiteMCCFlag = True

    return maxShieldNode, definiteMCCFlag
#******************************************************************


def findMaxDepthAndMinPath(thisPathDistanceAndLength):
    '''
    Purpose::
            To determine the maximum depth and min path for the headnode

    Input::
            tuple of dictionaries representing the shortest distance and paths for a node in the tree as returned by nx.single_source_dijkstra
            thisPathDistanceAndLength({distance}, {path})
                    {distance} = nodeAsString, valueAsInt, {path} = nodeAsString, pathAsList

    Output::
            tuple of the max pathLength and min pathDistance as a tuple (like what was input)
                    minDistanceAndMaxPath = ({distance},{path})
    '''
    maxPathLength = 0
    minPath = 0

    # maxPathLength for the node in question
    maxPathLength = max(
        len(values) for values in list(
            thisPathDistanceAndLength[1].values()))

    # if the duration is shorter then the min MCS length, then don't store!
    if maxPathLength < MIN_MCS_DURATION:  # MINIMUM_DURATION :
        minDistanceAndMaxPath = ()

    # else find the min path and max depth
    else:
        # max path distance for the node in question
        minPath = max(
            values for values in list(
                thisPathDistanceAndLength[0].values()))

        # check to determine the shortest path from the longest paths returned
        for pathDistance, path in zip(
            list(
                thisPathDistanceAndLength[0].values()), list(
                thisPathDistanceAndLength[1].values())):
            pathLength = len(path)
            # if pathLength is the same as the maxPathLength, then look the
            # pathDistance to determine if the min
            if pathLength == maxPathLength:
                if pathDistance <= minPath:
                    minPath = pathLength
                    # store details if absolute minPath and deepest
                    minDistanceAndMaxPath = (pathDistance, path)
    return minDistanceAndMaxPath
#******************************************************************


def thisDict(thisNode):
    '''
    Purpose::
            Return dictionary from graph if node exist in tree

    Input::
            thisNode: a string representing the CE to get the information for

    Output ::
            eachdict[1]: a dictionary representing the info associated with thisNode from the graph

    '''
    for eachdict in CLOUD_ELEMENT_GRAPH.nodes(thisNode):
        if eachdict[1]['uniqueID'] == thisNode:
            return eachdict[1]
#******************************************************************


def checkCriteria(thisCloudElementLatLon, aTemperature):
    '''
    Purpose::
            Determine if criteria B is met for a CEGraph

    Input::
            thisCloudElementLatLon: 2D array of (lat,lon) variable from the node dictionary being currently considered
            aTemperature:a integer representing the temperature maximum for masking

    Output ::
            cloudElementArea: a floating-point number representing the area in the array that meet the criteria - criteriaB

    '''
    cloudElementCriteriaBLatLon = []

    frame, CEcounter = ndimage.measurements.label(
        thisCloudElementLatLon, structure=STRUCTURING_ELEMENT)
    frameCEcounter = 0
    # determine min and max values in lat and lon, then use this to generate
    # teh array from LAT,LON meshgrid

    minLat = min(x[0] for x in thisCloudElementLatLon)
    maxLat = max(x[0]for x in thisCloudElementLatLon)
    minLon = min(x[1]for x in thisCloudElementLatLon)
    maxLon = max(x[1]for x in thisCloudElementLatLon)

    minLatIndex = np.argmax(LAT[:, 0] == minLat)
    maxLatIndex = np.argmax(LAT[:, 0] == maxLat)
    minLonIndex = np.argmax(LON[0, :] == minLon)
    maxLonIndex = np.argmax(LON[0, :] == maxLon)

    criteriaBframe = ma.zeros(
        ((abs(maxLatIndex - minLatIndex) + 1), (abs(maxLonIndex - minLonIndex) + 1)))

    for x in thisCloudElementLatLon:
        # to store the values of the subset in the new array, remove the minLatIndex and minLonindex from the
        # index given in the original array to get the indices for the new
        # array
        criteriaBframe[(np.argmax(LAT[:, 0] == x[0]) - minLatIndex),
                       (np.argmax(LON[0, :] == x[1]) - minLonIndex)] = x[2]

    # keep only those values < aTemperature
    tempMask = ma.masked_array(criteriaBframe, mask=(
        criteriaBframe >= aTemperature), fill_value=0)

    # get the actual values that the mask returned
    criteriaB = ma.zeros((criteriaBframe.shape)).astype('int16')

    for index, value in maenumerate(tempMask):
        lat_index, lon_index = index
        criteriaB[lat_index, lon_index] = value

    for count in range(CEcounter):
        #[0] is time dimension. Determine the actual values from the data
        # loc is a masked array
        #***** returns elements down then across thus (6,4) is 6 arrays deep of size 4
        try:

            loc = ndimage.find_objects(criteriaB)[0]
        except BaseException:            # this would mean that no objects were found meeting criteria B
            print("no objects at this temperature!")
            cloudElementArea = 0.0
            return cloudElementArea, cloudElementCriteriaBLatLon

        try:
            cloudElementCriteriaB = ma.zeros((criteriaB.shape))
            cloudElementCriteriaB = criteriaB[loc]
        except BaseException:
            print("YIKESS")
            print("CEcounter %s %s" % (CEcounter, criteriaB.shape))
            print("criteriaB %s" % criteriaB)

        for index, value in np.ndenumerate(cloudElementCriteriaB):
            if value != 0:
                t, lat, lon = index
                # add back on the minLatIndex and minLonIndex to find the true
                # lat, lon values
                lat_lon_tuple = (LAT[(lat), 0], LON[0, (lon)], value)
                cloudElementCriteriaBLatLon.append(lat_lon_tuple)

        cloudElementArea = np.count_nonzero(
            cloudElementCriteriaB) * XRES * YRES
        # do some cleaning up
        tempMask = []
        criteriaB = []
        cloudElementCriteriaB = []

        return cloudElementArea, cloudElementCriteriaBLatLon
#******************************************************************


def hasMergesOrSplits(nodeList):
    '''
    Purpose::
            Determine if nodes within a path defined from shortest_path splittingNodeDict
    Input::
            nodeList: list of strings representing the nodes from a path
    Output::
            splitList: a list of strings representing all the nodes in the path that split
            mergeList: a list of strings representing all the nodes in the path that merged
    '''
    mergeList = []
    splitList = []

    for node, numParents in list(PRUNED_GRAPH.in_degree(nodeList).items()):
        if numParents > 1:
            mergeList.append(node)

    for node, numChildren in list(PRUNED_GRAPH.out_degree(nodeList).items()):
        if numChildren > 1:
            splitList.append(node)
    # sort
    splitList.sort(key=lambda item: (
        len(item.split('C')[0]), item.split('C')[0]))
    mergeList.sort(key=lambda item: (
        len(item.split('C')[0]), item.split('C')[0]))

    return mergeList, splitList
#******************************************************************


def allAncestors(path, aNode):
    '''
    Purpose::
            Utility script to provide the path leading up to a nodeList

    Input::
            path: a list of strings representing the nodes in the path
            aNode: a string representing a node to be checked for parents

    Output::
            path: a list of strings representing the list of the nodes connected to aNode through its parents
            numOfChildren: an integer representing the number of parents of the node passed
    '''

    numOfParents = PRUNED_GRAPH.in_degree(aNode)
    try:
        if PRUNED_GRAPH.predecessors(aNode) and numOfParents <= 1:
            path = path + PRUNED_GRAPH.predecessors(aNode)
            thisNode = PRUNED_GRAPH.predecessors(aNode)[0]
            return allAncestors(path, thisNode)
        else:
            path = path + aNode
            return path, numOfParents
    except BaseException:
        return path, numOfParents
#******************************************************************


def allDescendants(path, aNode):
    '''
    Purpose::
            Utility script to provide the path leading up to a nodeList

    Input::
            path: a list of strings representing the nodes in the path
            aNode: a string representing a node to be checked for children

    Output::
            path: a list of strings representing the list of the nodes connected to aNode through its children
            numOfChildren: an integer representing the number of children of the node passed
    '''

    numOfChildren = PRUNED_GRAPH.out_degree(aNode)
    try:
        if PRUNED_GRAPH.successors(aNode) and numOfChildren <= 1:
            path = path + PRUNED_GRAPH.successors(aNode)
            thisNode = PRUNED_GRAPH.successors(aNode)[0]
            return allDescendants(path, thisNode)
        else:
            path = path + aNode
            #i.e. PRUNED_GRAPH.predecessors(aNode) is empty
            return path, numOfChildren
    # i.e. PRUNED_GRAPH.predecessors(aNode) threw an exception
    except BaseException:
        return path, numOfChildren
#******************************************************************


def addInfothisDict(thisNode, cloudElementArea, criteriaB):
    '''
    Purpose::
            Update original dictionary node with information

    Input::
            thisNode: a string representing the unique ID of a node
            cloudElementArea: a floating-point number representing the area of the cloud element
            criteriaB: a masked array of floating-point numbers representing the lat,lons meeting the criteria

    Output:: None
    '''
    for eachdict in CLOUD_ELEMENT_GRAPH.nodes(thisNode):
        if eachdict[1]['uniqueID'] == thisNode:
            eachdict[1]['CriteriaBArea'] = cloudElementArea
            eachdict[1]['CriteriaBLatLon'] = criteriaB
    return
#******************************************************************


def addNodeBehaviorIdentifier(thisNode, nodeBehaviorIdentifier):
    '''
    Purpose:: add an identifier to the node dictionary to indicate splitting, merging or neither node

    Input::
            thisNode: a string representing the unique ID of a node
            nodeBehaviorIdentifier: a string representing the behavior S- split, M- merge, B- both split and merge, N- neither split or merge

    Output :: None

    '''
    for eachdict in CLOUD_ELEMENT_GRAPH.nodes(thisNode):
        if eachdict[1]['uniqueID'] == thisNode:
            if 'nodeBehaviorIdentifier' not in list(eachdict[1].keys()):
                eachdict[1]['nodeBehaviorIdentifier'] = nodeBehaviorIdentifier
    return
#******************************************************************


def addNodeMCSIdentifier(thisNode, nodeMCSIdentifier):
    '''
    Purpose::
            Add an identifier to the node dictionary to indicate splitting, merging or neither node

    Input::
            thisNode: a string representing the unique ID of a node
            nodeMCSIdentifier: a string representing the stage of the MCS lifecyle  'I' for Initiation, 'M' for Maturity, 'D' for Decay

    Output :: None

    '''
    for eachdict in CLOUD_ELEMENT_GRAPH.nodes(thisNode):
        if eachdict[1]['uniqueID'] == thisNode:
            if 'nodeMCSIdentifier' not in list(eachdict[1].keys()):
                eachdict[1]['nodeMCSIdentifier'] = nodeMCSIdentifier
    return
#******************************************************************


def updateNodeMCSIdentifier(thisNode, nodeMCSIdentifier):
    '''
    Purpose::
            Update an identifier to the node dictionary to indicate splitting, merging or neither node

    Input::
            thisNode: thisNode: a string representing the unique ID of a node
            nodeMCSIdentifier: a string representing the stage of the MCS lifecyle  'I' for Initiation, 'M' for Maturity, 'D' for Decay

    Output :: None

    '''
    for eachdict in CLOUD_ELEMENT_GRAPH.nodes(thisNode):
        if eachdict[1]['uniqueID'] == thisNode:
            eachdict[1]['nodeMCSIdentifier'] = nodeBehaviorIdentifier

    return
#******************************************************************


def eccentricity(cloudElementLatLon):
    '''
    Purpose::
            Determines the eccentricity (shape) of contiguous boxes
            Values tending to 1 are more circular by definition, whereas
            values tending to 0 are more linear

    Input::
            cloudElementLatLon: 2D array in (lat,lon) representing T_bb contiguous squares

    Output::
            epsilon: a floating-point representing the eccentricity of the matrix passed

    '''

    epsilon = 0.0

    # loop over all lons and determine longest (non-zero) col
    # loop over all lats and determine longest (non-zero) row
    for latLon in cloudElementLatLon:
        # assign a matrix to determine the legit values

        nonEmptyLons = sum(sum(cloudElementLatLon) > 0)
        nonEmptyLats = sum(sum(cloudElementLatLon.transpose()) > 0)

        lonEigenvalues = 1.0 * nonEmptyLats / \
            (nonEmptyLons + 0.001)  # for long oval on y axis
        latEigenvalues = 1.0 * nonEmptyLons / \
            (nonEmptyLats + 0.001)  # for long oval on x-axs
        epsilon = min(latEigenvalues, lonEigenvalues)

    return epsilon
#******************************************************************


def cloudElementOverlap(currentCELatLons, previousCELatLons):
    '''
    Purpose::
            Determines the percentage overlap between two list of lat-lons passed

    Input::
            currentCELatLons: a list of tuples for the current CE
            previousCELatLons: a list of tuples for the other CE being considered

    Output::
            percentageOverlap: a floating-point representing the number of overlapping lat_lon tuples
            areaOverlap: a floating-point number representing the area overlapping

    '''

    latlonprev = []
    latloncurr = []
    count = 0
    percentageOverlap = 0.0
    areaOverlap = 0.0

    # remove the temperature from the tuples for currentCELatLons and
    # previousCELatLons then check for overlap
    latlonprev = [(x[0], x[1]) for x in previousCELatLons]
    latloncurr = [(x[0], x[1]) for x in currentCELatLons]

    # find overlap
    count = len(list(set(latloncurr) & set(latlonprev)))

    # find area overlap
    areaOverlap = count * XRES * YRES

    # find percentage
    percentageOverlap = max(((count * 1.0) / (len(latloncurr) * 1.0)),
                            ((count * 1.0) / (len(latlonprev) * 1.0)))

    return percentageOverlap, areaOverlap
#******************************************************************


def findCESpeed(node, MCSList):
    '''
    Purpose::
            To determine the speed of the CEs uses vector displacement delta_lat/delta_lon (y/x)

    Input::
            node: a string representing the CE
            MCSList: a list of strings representing the feature

    Output::
            CEspeed: a floating-point number representing the speed of the CE

    '''

    delta_lon = 0.0
    delta_lat = 0.0
    CEspeed = []
    theSpeed = 0.0

    theList = CLOUD_ELEMENT_GRAPH.successors(node)
    nodeLatLon = thisDict(node)['cloudElementCenter']

    for aNode in theList:
        if aNode in MCSList:
            # if aNode is part of the MCSList then determine distance
            aNodeLatLon = thisDict(aNode)['cloudElementCenter']
            # calculate CE speed
            # checking the lats
            # nodeLatLon[0] += 90.0
            # aNodeLatLon[0] += 90.0
            # delta_lat = (nodeLatLon[0] - aNodeLatLon[0])
            delta_lat = (
                (thisDict(node)['cloudElementCenter'][0] + 90.0) - (
                    thisDict(aNode)['cloudElementCenter'][0] + 90.0))
            # nodeLatLon[1] += 360.0
            # aNodeLatLon[1] += 360.0
            # delta_lon = (nodeLatLon[1] - aNodeLatLon[1])
            delta_lon = (
                (thisDict(node)['cloudElementCenter'][1] + 360.0) - (
                    thisDict(aNode)['cloudElementCenter'][1] + 360.0))

            try:
                # convert to s --> m/s
                theSpeed = abs(
                    (((delta_lat / delta_lon) * LAT_DISTANCE * 1000) / (TRES * 3600)))
            except BaseException:
                theSpeed = 0.0

            CEspeed.append(theSpeed)

            # print "~~~ ", thisDict(aNode)['uniqueID']
            # print "*** ", nodeLatLon, thisDict(node)['cloudElementCenter']
            # print "*** ", aNodeLatLon, thisDict(aNode)['cloudElementCenter']

    if not CEspeed:
        return 0.0
    else:
        return min(CEspeed)
#******************************************************************
#
#			UTILITY SCRIPTS FOR MCCSEARCH.PY
#
#******************************************************************


def maenumerate(mArray):
    '''
    Purpose::
            Utility script for returning the actual values from the masked array
            Taken from: http://stackoverflow.com/questions/8620798/numpy-ndenumerate-for-masked-arrays

    Input::
            mArray: the masked array returned from the ma.array() command


    Output::
            maskedValues: 3D (t,lat,lon), value of only masked values

    '''

    mask = ~mArray.mask.ravel()
    # beware yield fast, but generates a type called "generate" that does not
    # allow for array methods
    for index, maskedValue in zip(np.ndenumerate(mArray), mask):
        if maskedValue:
            yield index
#******************************************************************


def createMainDirectory(mainDirStr):
    '''
    Purpose::
            To create the main directory for storing information and
            the subdirectories for storing information
    Input::
            mainDir: a directory for where all information generated from
                    the program are to be stored
    Output:: None

    '''
    global MAINDIRECTORY

    MAINDIRECTORY = mainDirStr
    # if directory doesnt exist, creat it
    if not os.path.exists(MAINDIRECTORY):
        os.makedirs(MAINDIRECTORY)

    os.chdir((MAINDIRECTORY))
    # create the subdirectories
    try:
        os.makedirs('images')
        os.makedirs('textFiles')
        os.makedirs('MERGnetcdfCEs')
        os.makedirs('TRMMnetcdfCEs')
    except BaseException:
        print("Directory exists already!!!")
        # TODO: some nice way of prompting if it is ok to continue...or just
        # leave

    return
#******************************************************************


def checkForFiles(startTime, endTime, thisDir, fileType):
    '''
    Purpose:: To ensure all the files between the starttime and endTime
                      exist in the directory supplied

    Input::
                    startTime: a string yyyymmmddhh representing the starttime
                    endTime: a string yyyymmmddhh representing the endTime
                    thisDir: a string representing the directory path where to
                            look for the file
                    fileType: an integer representing the type of file in the directory
                            1 - MERG original files, 2 - TRMM original files

    Output::
                    status: a boolean representing whether all files exists

    '''
    filelist = []
    startFilename = ''
    endFilename = ''
    currFilename = ''
    status = False
    startyr = int(startTime[:4])
    startmm = int(startTime[4:6])
    startdd = int(startTime[6:8])
    starthr = int(startTime[-2:])
    endyr = int(endTime[:4])
    endmm = int(endTime[4:6])
    enddd = int(endTime[6:8])
    endhh = int(endTime[-2:])
    curryr = startyr
    currmm = startmm
    currdd = startdd
    currhr = starthr
    currmmStr = ''
    currddStr = ''
    currhrStr = ''
    endmmStr = ''
    endddStr = ''
    endhhStr = ''

    # check that the startTime is before the endTime
    if fileType == 1:
        # print "fileType is 1"
        startFilename = "merg_" + startTime + "_4km-pixel.nc"
        endFilename = thisDir + "/merg_" + endTime + "_4km-pixel.nc"

    if fileType == 2:
        # TODO:: determine closest time for TRMM files for end
        # http://disc.sci.gsfc.nasa.gov/additional/faq/precipitation_faq.shtml#convert
        # How do I extract time information from the TRMM 3B42 file name? section
        # startFilename = "3B42."+startTime[:8]+"."+currhr+".7A.nc"
        # endFilename = "3B42."+endTime[:8]+"."+endTime[-2:]+".7A.nc"
        if starthr % 3 == 2:
            currhr += 1
        elif starthr % 3 == 1:
            currhr -= 1
        else:
            currhr = starthr

        curryr, currmmStr, currddStr, currhrStr, _, _, _ = findTime(
            curryr, currmm, currdd, currhr)

        startFilename = "3B42." + \
            str(curryr) + currmmStr + currddStr + "." + currhrStr + ".7A.nc"
        if endhh % 3 == 2:
            endhh += 1
        elif endhh % 3 == 1:
            endhh -= 1

        endyr, endmmStr, endddStr, endhhStr, _, _, _ = findTime(
            endyr, endmm, enddd, endhh)

        endFilename = thisDir + "/3B42." + \
            str(endyr) + endmmStr + endddStr + "." + endhhStr + ".7A.nc"

    # check for files between startTime and endTime
    currFilename = thisDir + "/" + startFilename

    while currFilename is not endFilename:

        if not os.path.isfile(currFilename):
            print("file is missing! Filename: ", currFilename)
            status = False
            return status, filelist
        else:
            # create filelist
            filelist.append(currFilename)

        status = True
        if currFilename == endFilename:
            break

        # generate new currFilename
        if fileType == 1:
            currhr += 1
        elif fileType == 2:
            currhr += 3

        curryr, currmmStr, currddStr, currhrStr, currmm, currdd, currhr = findTime(
            curryr, currmm, currdd, currhr)

        if fileType == 1:
            currFilename = thisDir + "/" + "merg_" + \
                str(curryr) + currmmStr + currddStr + currhrStr + "_4km-pixel.nc"
        if fileType == 2:
            currFilename = thisDir + "/" + "3B42." + \
                str(curryr) + currmmStr + currddStr + "." + currhrStr + ".7A.nc"

    return status, filelist
#******************************************************************


def findTime(curryr, currmm, currdd, currhr):
    '''
    Purpose:: To determine the new yr, mm, dd, hr

    Input:: curryr, an integer representing the year
                    currmm, an integer representing the month
                    currdd, an integer representing the day
                    currhr, an integer representing the hour

    Output::curryr, an integer representing the year
                    currmm, an integer representing the month
                    currdd, an integer representing the day
                    currhr, an integer representing the hour
    '''
    if currhr > 23:
        currhr = 0
        currdd += 1
        if currdd > 30 and (currmm == 4 or currmm ==
                            6 or currmm == 9 or currmm == 11):
            currmm += 1
        elif currdd > 31 and (currmm == 1 or currmm == 3 or currmm == 5 or currmm == 7 or currmm == 8 or currmm == 10):
            currmm += 1
            currdd = 1
        elif currdd > 31 and currmm == 12:
            currmm = 1
            currdd = 1
            curryr += 1
        elif currdd > 28 and currmm == 2 and (curryr % 4) != 0:
            currmm = 3
            currdd = 1
        elif (curryr % 4) == 0 and currmm == 2 and currdd > 29:
            currmm = 3
            currdd = 1

    if currmm < 10:
        currmmStr = "0" + str(currmm)
    else:
        currmmStr = str(currmm)

    if currdd < 10:
        currddStr = "0" + str(currdd)
    else:
        currddStr = str(currdd)

    if currhr < 10:
        currhrStr = "0" + str(currhr)
    else:
        currhrStr = str(currhr)

    return curryr, currmmStr, currddStr, currhrStr, currmm, currdd, currhr
#******************************************************************


def find_nearest(thisArray, value):
    '''
    Purpose :: to determine the value within an array closes to
                    another value

    Input ::
    Output::
    '''
    idx = (np.abs(thisArray - value)).argmin()
    return thisArray[idx]
#******************************************************************

#******************************************************************


def postProcessingNetCDF(dataset, dirName=None):
    '''
    Purpose::
            Utility script displaying the data in NETCDF4 files

    Input::
            dataset: integer representing original MERG (1) or post-processed MERG data (2) or post-processed TRMM(3)
            string: Directory to the location of the raw (MERG) files, preferably zipped

    Output::
       Generates 2D plots in location as specfied in the code

    '''
    coreDir = os.path.dirname(os.path.abspath(__file__))
    imgFilename = ''

    if dataset == 1:
        var = 'ch4'
        plotTitle = 'Original MERG data '
    elif dataset == 2:
        var = 'brightnesstemp'
        plotTitle = 'MERG CE data'
    elif dataset == 3:
        var = 'precipitation_Accumulation'
        plotTitle = 'TRMM CE data'

    # sort files
    os.chdir((dirName + '/'))
    files = list(filter(os.path.isfile, glob.glob("*.nc")))
    files.sort(key=lambda x: os.path.getmtime(x))

    for eachfile in files:
        fullFname = os.path.splitext(eachfile)[0]
        fnameNoExtension = fullFname.split('.nc')[0]

        fname = dirName + '/' + fnameNoExtension + '.nc'

        if os.path.isfile(fname):
            fileData = Dataset(fname, 'r', format='NETCDF4')
            file_variable = fileData.variables[var][:]
            lats = fileData.variables['latitude'][:]
            lons = fileData.variables['longitude'][:]
            LONDATA, LATDATA = np.meshgrid(lons, lats)
            nygrd = len(LATDATA[:, 0])
            nxgrd = len(LONDATA[0, :])
            fileData.close()

        imgFilename = MAINDIRECTORY + '/images/' + fnameNoExtension + '.gif'

        if dataset == 3:
            createPrecipPlot(np.squeeze(file_variable, axis=0),
                             LATDATA[:, 0], LONDATA[0, :], plotTitle, imgFilename)
        else:
            plotter.draw_contour_map(
                file_variable, LATDATA[:, 0], LONDATA[0, :], imgFilename, ptitle=plotTitle)

    return
#******************************************************************


def drawGraph(thisGraph, graphTitle, edgeWeight=None):
    '''
    Purpose::
            Utility function to draw graph in the hierachial format

    Input::
            thisGraph: a Networkx directed graph
            graphTitle: a string representing the graph title
            edgeWeight: (optional) a list of integers representing the edge weights in the graph

    Output:: None

    '''

    imgFilename = MAINDIRECTORY + '/images/' + graphTitle + ".gif"
    fig = plt.figure(facecolor='white', figsize=(16, 12))

    edge95 = [
        (u, v) for (
            u, v, d) in thisGraph.edges(
            data=True) if d['weight'] == edgeWeight[0]]
    edge90 = [
        (u, v) for (
            u, v, d) in thisGraph.edges(
            data=True) if d['weight'] == edgeWeight[1]]
    edegeOverlap = [
        (u, v) for (
            u, v, d) in thisGraph.edges(
            data=True) if d['weight'] == edgeWeight[2]]

    nx.write_dot(thisGraph, 'test.dot')
    plt.title(graphTitle)
    pos = nx.graphviz_layout(thisGraph, prog='dot')
    # draw graph in parts
    # nodes
    nx.draw_networkx_nodes(thisGraph, pos, with_labels=True, arrows=False)
    # edges
    nx.draw_networkx_edges(
        thisGraph,
        pos,
        edgelist=edge95,
        alpha=0.5,
        arrows=False)
    nx.draw_networkx_edges(
        thisGraph,
        pos,
        edgelist=edge90,
        edge_color='b',
        style='dashed',
        arrows=False)
    nx.draw_networkx_edges(
        thisGraph,
        pos,
        edgelist=edegeOverlap,
        edge_color='y',
        style='dashed',
        arrows=False)
    # labels
    nx.draw_networkx_labels(thisGraph, pos, arrows=False)
    plt.axis('off')
    plt.savefig(imgFilename, facecolor=fig.get_facecolor(), transparent=True)
    # do some clean up...and ensuring that we are in the right dir
    os.chdir((MAINDIRECTORY + '/'))
    subprocess.call('rm test.dot', shell=True)
#******************************************************************


def getModelTimes(xtimes, timeVarName):
    '''
    Taken from process.py, removed the file opening at the beginning
    TODO:  Do a better job handling dates here
    Routine to convert from model times ('hours since 1900...', 'days since ...')
    into a python datetime structure

    Input::
        modelFile - path to the model tile you want to extract the times list and modelTimeStep from
        timeVarName - name of the time variable in the model file

    Output::
        times  - list of python datetime objects describing model data times
        modelTimeStep - 'hourly','daily','monthly','annual'
    '''

    from ocw import utils

    timeFormat = xtimes.units
    # search to check if 'since' appears in units
    try:
        sinceLoc = re.search('since', timeFormat).end()

    except AttributeError:
        print("Error decoding model times: time variable attributes do not contain 'since'")
        raise

    units = None
    TIME_UNITS = ('minutes', 'hours', 'days', 'months', 'years')
    # search for 'seconds','minutes','hours', 'days', 'months', 'years' so
    # know units
    for unit in TIME_UNITS:
        if re.search(unit, timeFormat):
            units = unit
            break

    # cut out base time (the bit following 'since')
    base_time_string = string.lstrip(timeFormat[sinceLoc:])
    # decode base time
    base_time = decodeTimeFromString(base_time_string)

    times = []

    for xtime in xtimes[:]:
        # Cast time as an int
        # TODO: KDW this may cause problems for data that is hourly with more
        # than one timestep in it
        xtime = int(xtime)

        if int(xtime) == 0:
            xtime = 1

        if units == 'minutes':
            dt = timedelta(minutes=xtime)
            new_time = base_time + dt
        elif units == 'hours':
            dt = timedelta(hours=int(xtime))
            new_time = base_time + dt  # timedelta(hours=int(xtime))
        elif units == 'days':
            dt = timedelta(days=xtime)
            new_time = base_time + dt
        elif units == 'months':
            # NB. adding months in python is complicated as month length varies and hence ambiguous.
            # Perform date arithmatic manually
            #  Assumption: the base_date will usually be the first of the month
            #              NB. this method will fail if the base time is on the 29th or higher day of month
            #                      -as can't have, e.g. Feb 31st.
            new_month = int(base_time.month + xtime % 12)
            new_year = int(math.floor(base_time.year + xtime / 12.))
            new_time = datetime.datetime(
                new_year,
                new_month,
                base_time.day,
                base_time.hour,
                base_time.second,
                0)
        elif units == 'years':
            dt = datetime.timedelta(years=xtime)
            new_time = base_time + dt

        times.append(new_time)

    try:
        if len(xtimes) == 1:
            timeStepLength = 0
        else:
            timeStepLength = int(xtimes[1] - xtimes[0] + 1.e-12)

        modelTimeStep = getModelTimeStep(units, timeStepLength)

        # if timeStepLength is zero do not normalize times as this would create
        # an empty list for MERG (hourly) data
        if timeStepLength != 0:
            times = normalizeDatetimes(times, modelTimeStep)
    except BaseException:
        raise

    return times, modelTimeStep
#******************************************************************


def getModelTimeStep(units, stepSize):
    # Time units are now determined. Determine the time intervals of input
    # data (mdlTimeStep)
    '''
    Taken from process.py
    '''
    if units == 'minutes':
        if stepSize == 60:
            modelTimeStep = 'hourly'
        elif stepSize == 1440:
            modelTimeStep = 'daily'
        # 28 days through 31 days
        elif 40320 <= stepSize <= 44640:
            modelTimeStep = 'monthly'
        # 365 days through 366 days
        elif 525600 <= stepSize <= 527040:
            modelTimeStep = 'annual'
        else:
            raise Exception(
                'model data time step interval exceeds the max time interval (annual)',
                units,
                stepSize)

    elif units == 'hours':
        # need a check for fractional hrs and only one hr i.e. stepSize=0 e.g.
        # with MERG data
        if stepSize == 0 or stepSize == 1:
            modelTimeStep = 'hourly'
        elif stepSize == 24:
            modelTimeStep = 'daily'
        elif 672 <= stepSize <= 744:
            modelTimeStep = 'monthly'
        elif 8760 <= stepSize <= 8784:
            modelTimeStep = 'annual'
        else:
            raise Exception(
                'model data time step interval exceeds the max time interval (annual)',
                units,
                stepSize)

    elif units == 'days':
        if stepSize == 1:
            modelTimeStep = 'daily'
        elif 28 <= stepSize <= 31:
            modelTimeStep = 'monthly'
        elif 365 <= stepSize <= 366:
            modelTimeStep = 'annual'
        else:
            raise Exception(
                'model data time step interval exceeds the max time interval (annual)',
                units,
                stepSize)

    elif units == 'months':
        if stepSize == 1:
            modelTimeStep = 'monthly'
        elif stepSize == 12:
            modelTimeStep = 'annual'
        else:
            raise Exception(
                'model data time step interval exceeds the max time interval (annual)',
                units,
                stepSize)

    elif units == 'years':
        if stepSize == 1:
            modelTimeStep = 'annual'
        else:
            raise Exception(
                'model data time step interval exceeds the max time interval (annual)',
                units,
                stepSize)

    else:
        errorMessage = 'the time unit ', units, ' is not currently handled in this version.'
        raise Exception(errorMessage)

    return modelTimeStep
#******************************************************************


def decodeTimeFromString(time_string):
    '''
    Taken from process.py
     Decodes string into a python datetime object
     *Method:* tries a bunch of different time format possibilities and hopefully one of them will hit.
     ::

       **Input:**  time_string - a string that represents a date/time

       **Output:** mytime - a python datetime object
    '''
    # This will deal with times that use decimal seconds
    if '.' in time_string:
        time_string = time_string.split('.')[0] + '0'
    else:
        pass

    try:
        mytime = datetime.strptime(time_string, '%Y-%m-%d %H')
        return mytime

    except ValueError:
        pass

    print("Error decoding time string: string does not match a predefined time format'")
    return 0
#******************************************************************


def do_regrid(q, lat, lon, lat2, lon2, order=1, mdi=-999999999):
    """
    This function has been moved to the ocw/dataset_processor module
    """
    from ocw import dataset_processor
    q2 = dataset_processor._rcmes_spatial_regrid(
        q, lat, lon, lat2, lon2, order=1)

    return q2
#******************************************************************
#
#			 METRICS FUNCTIONS FOR MERG.PY
# TODO: rewrite these metrics so that they use the data from the
#	file instead of the graph as this reduce mem resources needed
#
#
#******************************************************************


def numberOfFeatures(finalMCCList):
    '''
    Purpose::
            To count the number of MCCs found for the period

    Input::
            finalMCCList: a list of list of strings representing a list of list of nodes representing a MCC

    Output::
            an integer representing the number of MCCs found

    '''
    return len(finalMCCList)
#******************************************************************


def temporalAndAreaInfoMetric(finalMCCList):
    '''
    Purpose::
            To provide information regarding the temporal properties of the MCCs found

    Input::
            finalMCCList: a list of dictionaries representing a list of nodes representing a MCC

    Output::
            allMCCtimes: a list of dictionaries {MCCtimes, starttime, endtime, duration, area} representing a list of dictionaries
                    of MCC temporal details for each MCC in the period considered

    Assumptions::
            the final time hour --> the event lasted throughout that hr, therefore +1 to endtime
    '''
    # TODO: in real data edit this to use datetime
    #starttime =0
    #endtime =0
    #duration = 0
    MCCtimes = []
    allMCCtimes = []
    MCSArea = []

    if finalMCCList:
        for eachMCC in finalMCCList:
            # get the info from the node
            for eachNode in eachMCC:
                MCCtimes.append(thisDict(eachNode)['cloudElementTime'])
                MCSArea.append(thisDict(eachNode)['cloudElementArea'])

            # sort and remove duplicates
            MCCtimes = sorted(set(MCCtimes))
            tdelta = MCCtimes[1] - MCCtimes[0]
            starttime = MCCtimes[0]
            endtime = MCCtimes[-1]
            duration = (endtime - starttime) + tdelta
            print(
                "starttime: %s endtime: %s tdelta: %s duration: %s MCSAreas %s" %
                (starttime, endtime, tdelta, duration, MCSArea))
            allMCCtimes.append({'MCCtimes': MCCtimes,
                                'starttime': starttime,
                                'endtime': endtime,
                                'duration': duration,
                                'MCSArea': MCSArea})
            MCCtimes = []
            MCSArea = []
    else:
        allMCCtimes = []
        tdelta = 0

    return allMCCtimes, tdelta
#******************************************************************


def longestDuration(allMCCtimes):
    '''
    Purpose::
            To determine the longest MCC for the period

    Input::
            allMCCtimes: a list of dictionaries {MCCtimes, starttime, endtime, duration, area} representing a list of dictionaries
                    of MCC temporal details for each MCC in the period considered

    Output::
            an integer - lenMCC: representing the duration of the longest MCC found
               a list of strings - longestMCC: representing the nodes of longest MCC

    Assumptions::

    '''

    # MCCList = []
    # lenMCC = 0
    # longestMCC =[]

    # #remove duplicates
    # MCCList = list(set(finalMCCList))

    # longestMCC = max(MCCList, key = lambda tup:len(tup))
    # lenMCC = len(longestMCC)

    # return lenMCC, longestMCC

    return max([MCC['duration'] for MCC in allMCCtimes])
#******************************************************************


def shortestDuration(allMCCtimes):
    '''
    Purpose:: To determine the shortest MCC for the period

    Input:: list of dictionaries - allMCCtimes {MCCtimes, starttime, endtime, duration): a list of dictionaries
                      of MCC temporal details for each MCC in the period considered

    Output::an integer - lenMCC: representing the duration of the shortest MCC found
                    a list of strings - longestMCC: representing the nodes of shortest MCC

    Assumptions::

    '''
    # lenMCC = 0
    # shortestMCC =[]
    # MCCList =[]

    # #remove duplicates
    # MCCList = list(set(finalMCCList))

    # shortestMCC = min(MCCList, key = lambda tup:len(tup))
    # lenMCC = len(shortestMCC)

    # return lenMCC, shortestMCC
    return min([MCC['duration'] for MCC in allMCCtimes])
#******************************************************************


def averageDuration(allMCCtimes):
    '''
    Purpose:: To determine the average MCC length for the period

    Input:: list of dictionaries - allMCCtimes {MCCtimes, starttime, endtime, duration): a list of dictionaries
                      of MCC temporal details for each MCC in the period considered

    Output::a floating-point representing the average duration of a MCC in the period

    Assumptions::

    '''

    return sum([MCC['duration'] for MCC in allMCCtimes],
               timedelta(seconds=0)) / len(allMCCtimes)
#******************************************************************


def averageTime(allTimes):
    '''
    Purpose::
            To determine the average time in a list of datetimes
            e.g. of use is finding avg starttime,
    Input::
            allTimes: a list of datetimes representing all of a given event e.g. start time

    Output::
            a floating-point number representing the average of the times given

    '''
    avgTime = 0

    for aTime in allTimes:
        avgTime += aTime.second + 60 * aTime.minute + 3600 * aTime.hour

    if len(allTimes) > 1:
        avgTime /= len(allTimes)

    rez = str(avgTime / 3600) + ' ' + str((avgTime %
                                           3600) / 60) + ' ' + str(avgTime % 60)
    return datetime.strptime(rez, "%H %M %S")
#******************************************************************


def averageFeatureSize(finalMCCList):
    '''
    Purpose:: To determine the average MCC size for the period

    Input:: a list of list of strings - finalMCCList: a list of list of nodes representing a MCC

    Output::a floating-point representing the average area of a MCC in the period

    Assumptions::

    '''
    thisMCC = 0.0
    thisMCCAvg = 0.0

    # for each node in the list, get the are information from the dictionary
    # in the graph and calculate the area
    for eachPath in finalMCCList:
        for eachNode in eachPath:
            thisMCC += thisDict(eachNode)['cloudElementArea']

        thisMCCAvg += (thisMCC / len(eachPath))
        thisMCC = 0.0

    # calcuate final average
    return thisMCCAvg / (len(finalMCCList))
#******************************************************************


def commonFeatureSize(finalMCCList):
    '''
    Purpose::
            To determine the common (mode) MCC size for the period

    Input::
            finalMCCList: a list of list of strings representing the list of nodes representing a MCC

    Output::
            a floating-point representing the average area of a MCC in the period

    Assumptions::

    '''
    thisMCC = 0.0
    thisMCCAvg = []

    # for each node in the list, get the area information from the dictionary
    # in the graph and calculate the area
    for eachPath in finalMCCList:
        for eachNode in eachPath:
            thisMCC += eachNode['cloudElementArea']

        thisMCCAvg.append(thisMCC / len(eachPath))
        thisMCC = 0.0

    # calcuate
    hist, bin_edges = np.histogram(thisMCCAvg)
    return hist, bin_edges
#******************************************************************


def precipTotals(finalMCCList):
    '''
    Purpose::
            Precipitation totals associated with a cloud element

    Input::
            finalMCCList: a list of dictionaries representing a list of nodes representing a MCC

    Output::
            precipTotal: a floating-point number representing the total amount of precipitation associated
                    with the feature
    '''
    precipTotal = 0.0
    CEprecip = 0.0
    MCSPrecip = []
    allMCSPrecip = []
    count = 0

    if finalMCCList:
        # print "len finalMCCList is: ", len(finalMCCList)
        for eachMCC in finalMCCList:
            # get the info from the node
            for node in eachMCC:
                eachNode = thisDict(node)
                count += 1
                if count == 1:
                    prevHr = int(
                        str(eachNode['cloudElementTime']).replace(" ", "")[-8:-6])

                currHr = int(
                    str(eachNode['cloudElementTime']).replace(" ", "")[-8:-6])
                if prevHr == currHr:
                    CEprecip += eachNode['cloudElementPrecipTotal']
                else:
                    MCSPrecip.append((prevHr, CEprecip))
                    CEprecip = eachNode['cloudElementPrecipTotal']
                # last value in for loop
                if count == len(eachMCC):
                    MCSPrecip.append((currHr, CEprecip))

                precipTotal += eachNode['cloudElementPrecipTotal']
                prevHr = currHr

            MCSPrecip.append(('0', precipTotal))

            allMCSPrecip.append(MCSPrecip)
            precipTotal = 0.0
            CEprecip = 0.0
            MCSPrecip = []
            count = 0

        print("allMCSPrecip %s" % allMCSPrecip)

    return allMCSPrecip
#******************************************************************


def precipMaxMin(finalMCCList):
    '''
    TODO: this doesnt work the np.min/max function seems to be not working with the nonzero option..possibly a problem upstream with cloudElementLatLonTRMM
    Purpose::
            Precipitation maximum and min rates associated with each CE in MCS
    Input::
            finalMCCList: a list of dictionaries representing a list of nodes representing a MCC

    Output::
            MCSPrecip: a list indicating max and min rate for each CE identified

    '''
    maxCEprecip = 0.0
    minCEprecip = 0.0
    MCSPrecip = []
    allMCSPrecip = []

    if finalMCCList:
        if isinstance(finalMCCList[0], str):  # len(finalMCCList) == 1:
            for node in finalMCCList:
                eachNode = thisDict(node)
                CETRMM = eachNode['cloudElementLatLonTRMM']

                print("all %s" % np.min(CETRMM[np.nonzero(CETRMM)]))
                print("minCEprecip %s" %
                      np.min(eachNode['cloudElementLatLonTRMM']))

                print("maxCEprecip %s" % np.max(eachNode['cloudElementLatLonTRMM'][np.nonzero(
                    eachNode['cloudElementLatLonTRMM'])]))
                sys.exit()
                maxCEprecip = np.max(eachNode['cloudElementLatLonTRMM'][np.nonzero(
                    eachNode['cloudElementLatLonTRMM'])])
                minCEprecip = np.min(eachNode['cloudElementLatLonTRMM'][np.nonzero(
                    eachNode['cloudElementLatLonTRMM'])])
                MCSPrecip.append(
                    (eachNode['uniqueID'], minCEprecip, maxCEprecip))

        else:
            for eachMCC in finalMCCList:
                # get the info from the node
                for node in eachMCC:
                    eachNode = thisDict(node)
                    # find min and max precip
                    maxCEprecip = np.max(eachNode['cloudElementLatLonTRMM'][np.nonzero(
                        eachNode['cloudElementLatLonTRMM'])])
                    minCEprecip = np.min(eachNode['cloudElementLatLonTRMM'][np.nonzero(
                        eachNode['cloudElementLatLonTRMM'])])
                    MCSPrecip.append(
                        (eachNode['uniqueID'], minCEprecip, maxCEprecip))
                allMCSPrecip.append(MCSPrecip)
                MCSPrecip = []

    return MCSPrecip
#******************************************************************
#
#							PLOTS
#
#******************************************************************


def displaySize(finalMCCList):
    '''
    Purpose::
            To create a figure showing the area verse time for each MCS

    Input::
            finalMCCList: a list of list of strings representing the list of nodes representing a MCC

    Output::
            None

    '''
    timeList = []
    count = 1
    imgFilename = ''
    minArea = 10000.0
    maxArea = 0.0
    eachNode = {}

    # for each node in the list, get the area information from the dictionary
    # in the graph and calculate the area

    if finalMCCList:
        for eachMCC in finalMCCList:
            # get the info from the node
            for node in eachMCC:
                eachNode = thisDict(node)
                timeList.append(eachNode['cloudElementTime'])

                if eachNode['cloudElementArea'] < minArea:
                    minArea = eachNode['cloudElementArea']
                if eachNode['cloudElementArea'] > maxArea:
                    maxArea = eachNode['cloudElementArea']

            # sort and remove duplicates
            timeList = sorted(set(timeList))
            tdelta = timeList[1] - timeList[0]
            starttime = timeList[0] - tdelta
            endtime = timeList[-1] + tdelta
            timeList.insert(0, starttime)
            timeList.append(endtime)

            # plot info
            plt.close('all')
            title = 'Area distribution of the MCC over somewhere'
            # figsize=(10,8))#figsize=(16,12))
            fig = plt.figure(facecolor='white', figsize=(18, 10))
            fig, ax = plt.subplots(1, facecolor='white', figsize=(10, 10))

            # the data
            for node in eachMCC:  # for eachNode in eachMCC:
                eachNode = thisDict(node)
                if eachNode['cloudElementArea'] < 80000:  # 2400.00:
                    ax.plot(
                        eachNode['cloudElementTime'],
                        eachNode['cloudElementArea'],
                        'bo',
                        markersize=10)
                elif eachNode['cloudElementArea'] >= 80000.00 and eachNode['cloudElementArea'] < 160000.00:
                    ax.plot(
                        eachNode['cloudElementTime'],
                        eachNode['cloudElementArea'],
                        'yo',
                        markersize=20)
                else:
                    ax.plot(
                        eachNode['cloudElementTime'],
                        eachNode['cloudElementArea'],
                        'ro',
                        markersize=30)

            #axes and labels
            maxArea += 1000.00
            ax.set_xlim(starttime, endtime)
            ax.set_ylim(minArea, maxArea)
            ax.set_ylabel('Area in km^2', fontsize=12)
            ax.set_title(title)
            ax.fmt_xdata = mdates.DateFormatter('%Y-%m-%d%H:%M:%S')
            fig.autofmt_xdate()
            plt.subplots_adjust(bottom=0.2)

            imgFilename = MAINDIRECTORY + '/images/' + str(count) + 'MCS.gif'
            plt.savefig(
                imgFilename,
                facecolor=fig.get_facecolor(),
                transparent=True)

            # if time in not already in the time list, append it
            timeList = []
            count += 1
    return
#******************************************************************


def displayPrecip(finalMCCList):
    '''
    Purpose::
            To create a figure showing the precip rate verse time for each MCS

    Input::
            finalMCCList: a list of dictionaries representing a list of nodes representing a MCC

    Output:: None

    '''
    timeList = []
    oriTimeList = []
    colorBarTime = []
    count = 1
    imgFilename = ''
    TRMMprecipDis = []
    percentagePrecipitating = []  # 0.0
    CEArea = []
    nodes = []
    xy = []
    x = []
    y = []
    precip = []
    partialArea = []
    totalSize = 0.0

    firstTime = True
    xStart = 0.0
    yStart = 0.0

    num_bins = 5

    # for each node in the list, get the area information from the dictionary
    # in the graph and calculate the area

    if finalMCCList:
        for eachMCC in finalMCCList:
            # get the info from the node
            for node in eachMCC:
                eachNode = thisDict(node)
                if firstTime:
                    xStart = eachNode['cloudElementCenter'][1]  # lon
                    yStart = eachNode['cloudElementCenter'][0]  # lat
                timeList.append(eachNode['cloudElementTime'])
                percentagePrecipitating.append(
                    (eachNode['TRMMArea'] / eachNode['cloudElementArea']) * 100.0)
                CEArea.append(eachNode['cloudElementArea'])
                nodes.append(eachNode['uniqueID'])
                # print eachNode['uniqueID'],
                # eachNode['cloudElementCenter'][1],
                # eachNode['cloudElementCenter'][0]
                x.append(eachNode['cloudElementCenter'][1])  # -xStart)
                y.append(eachNode['cloudElementCenter'][0])  # -yStart)

                firstTime = False

            # convert the timeList[] to list of floats
            for i in range(len(timeList)):  # oriTimeList:
                colorBarTime.append(time.mktime(timeList[i].timetuple()))

            totalSize = sum(CEArea)
            partialArea = [(a / totalSize) * 30000 for a in CEArea]

            # print "x ", x
            # print "y ", y

            # plot info
            plt.close('all')

            title = 'Precipitation distribution of the MCS '
            fig, ax = plt.subplots(1, facecolor='white', figsize=(20, 7))

            cmap = plt.jet
            ax.scatter(
                x,
                y,
                s=partialArea,
                c=colorBarTime,
                edgecolors='none',
                marker='o',
                cmap=cmap)
            colorBarTime = []
            colorBarTime = sorted(set(timeList))
            cb = colorbar_index(
                ncolors=len(colorBarTime),
                nlabels=colorBarTime,
                cmap=cmap)

            #axes and labels
            ax.set_xlabel('Degrees Longtude', fontsize=12)
            ax.set_ylabel('Degrees Latitude', fontsize=12)
            ax.set_title(title)
            ax.grid(True)
            plt.subplots_adjust(bottom=0.2)

            for i, txt in enumerate(nodes):
                if CEArea[i] >= 2400.00:
                    ax.annotate('%d' %
                                percentagePrecipitating[i] +
                                '%', (x[i], y[i]))
                precip = []

            imgFilename = MAINDIRECTORY + \
                '/images/MCSprecip' + str(count) + '.gif'
            plt.savefig(
                imgFilename,
                facecolor=fig.get_facecolor(),
                transparent=True)

            # reset for next image
            timeList = []
            percentagePrecipitating = []
            CEArea = []
            x = []
            y = []
            colorBarTime = []
            nodes = []
            precip = []
            count += 1
            firstTime = True
    return
#******************************************************************


def plotPrecipHistograms(finalMCCList, num_bins=5):
    '''
    Purpose::
            To create plots (histograms) of the each TRMMnetcdfCEs files

    Input::
            finalMCCList: a list of dictionaries representing a list of nodes representing a MCC
            num_bins: an integer representing the number of bins
    Output::
            plots
    '''

    precip = []
    imgFilename = " "
    lastTime = " "
    firstTime = True
    MCScount = 0
    MSClen = 0
    thisCount = 0
    totalPrecip = np.zeros((1, 137, 440))

    # TODO: use try except block instead
    if finalMCCList:

        for eachMCC in finalMCCList:
            firstTime = True
            MCScount += 1
            # totalPrecip=np.zeros((1,137,440))
            totalPrecip = np.zeros((1, 413, 412))

            # get the info from the node
            for node in eachMCC:
                eachNode = thisDict(node)
                thisTime = eachNode['cloudElementTime']
                MCSlen = len(eachMCC)
                thisCount += 1

                # this is the precipitation distribution plot from
                # displayPrecip

                if eachNode['cloudElementArea'] >= 2400.0:
                    if (str(thisTime) != lastTime and lastTime !=
                            " ") or thisCount == MCSlen:
                        # plt.close('all')
                        # title = 'TRMM precipitation distribution for '+ str(thisTime)

                        # fig,ax = plt.subplots(1, facecolor='white', figsize=(7,5))

                        # n,binsdg = np.histogram(precip, num_bins)
                        # wid = binsdg[1:] - binsdg[:-1]
                        # plt.bar(binsdg[:-1], n/float(len(precip)), width=wid)

                        # #make percentage plot
                        # formatter = FuncFormatter(to_percent)
                        # plt.xlim(min(binsdg), max(binsdg))
                        # ax.set_xticks(binsdg)
                        # ax.set_xlabel('Precipitation [mm]', fontsize=12)
                        # ax.set_ylabel('Area', fontsize=12)
                        # ax.set_title(title)
                        # # Set the formatter
                        # plt.gca().yaxis.set_major_formatter(formatter)
                        # plt.gca().xaxis.set_major_formatter(FormatStrFormatter('%0.0f'))
                        # imgFilename = MAINDIRECTORY+'/images/'+str(thisTime)+eachNode['uniqueID']+'TRMMMCS.gif'

                        # plt.savefig(imgFilename, transparent=True)
                        data_names = ['Precipitation [mm]', 'Area']
                        plotter.draw_histogram(
                            precip, data_names, imgFilename, num_bins)
                        precip = []

                    # ------ NETCDF File get info -----------------------------
                    thisFileName = MAINDIRECTORY + '/TRMMnetcdfCEs/TRMM' + \
                        str(thisTime).replace(" ", "_") + eachNode['uniqueID'] + '.nc'
                    TRMMData = Dataset(thisFileName, 'r', format='NETCDF4')
                    precipRate = TRMMData.variables['precipitation_Accumulation'][:, :, :]
                    CEprecipRate = precipRate[0, :, :]
                    TRMMData.close()
                    if firstTime:
                        totalPrecip = np.zeros((CEprecipRate.shape))

                    totalPrecip = np.add(totalPrecip, precipRate)
                    # ------ End NETCDF File ----------------------------------
                    for index, value in np.ndenumerate(CEprecipRate):
                        if value != 0.0:
                            precip.append(value)

                    lastTime = str(thisTime)
                    firstTime = False
                else:
                    lastTime = str(thisTime)
                    firstTime = False
    return
#******************************************************************


def plotAccTRMM(finalMCCList):
    '''
    Purpose::
            (1) generate a file with the accumulated precipiation for the MCS
            (2) generate the appropriate image
            TODO: NB: as the domain changes, will need to change XDEF and YDEF by hand to accomodate the new domain
            TODO: look into getting the info from the NETCDF file

    Input::
            finalMCCList: a list of dictionaries representing a list of nodes representing a MCC

    Output::
            a netcdf file containing the accumulated precip
            a 2D matlablibplot
    '''
    os.chdir((MAINDIRECTORY + '/TRMMnetcdfCEs'))
    fname = ''
    imgFilename = ''
    firstPartName = ''
    firstTime = True
    replaceExpXDef = ''

    # generate the file name using MCCTimes
    # if the file name exists, add it to the accTRMM file
    for path in finalMCCList:
        for eachNode in path:
            thisNode = thisDict(eachNode)
            fname = 'TRMM' + \
                str(thisNode['cloudElementTime']).replace(" ", "_") + thisNode['uniqueID'] + '.nc'

            if os.path.isfile(fname):
                # open NetCDF file add info to the accu
                # print "opening TRMM file ", fname
                TRMMCEData = Dataset(fname, 'r', format='NETCDF4')
                precipRate = TRMMCEData.variables['precipitation_Accumulation'][:]
                lats = TRMMCEData.variables['latitude'][:]
                lons = TRMMCEData.variables['longitude'][:]
                LONTRMM, LATTRMM = np.meshgrid(lons, lats)
                nygrdTRMM = len(LATTRMM[:, 0])
                nxgrdTRMM = len(LONTRMM[0, :])
                precipRate = ma.masked_array(
                    precipRate, mask=(precipRate < 0.0))
                TRMMCEData.close()

                if firstTime:
                    firstPartName = str(
                        thisNode['cloudElementTime']).replace(
                        " ", "_") + '-'
                    accuPrecipRate = ma.zeros((precipRate.shape))
                    firstTime = False

                accuPrecipRate += precipRate

        imgFilename = MAINDIRECTORY + '/images/MCSaccu' + firstPartName + \
            str(thisNode['cloudElementTime']).replace(" ", "_") + '.gif'

    # create new netCDF file
    accuTRMMFile = MAINDIRECTORY + '/TRMMnetcdfCEs/accu' + firstPartName + \
        str(thisNode['cloudElementTime']).replace(" ", "_") + '.nc'
    # write the file
    accuTRMMData = Dataset(accuTRMMFile, 'w', format='NETCDF4')
    accuTRMMData.description = 'Accumulated precipitation data'
    accuTRMMData.calendar = 'standard'
    accuTRMMData.conventions = 'COARDS'
    # dimensions
    accuTRMMData.createDimension('time', None)
    accuTRMMData.createDimension('lat', nygrdTRMM)
    accuTRMMData.createDimension('lon', nxgrdTRMM)

    # variables
    TRMMprecip = ('time', 'lat', 'lon',)
    times = accuTRMMData.createVariable('time', 'f8', ('time',))
    times.units = 'hours since ' + \
        str(thisNode['cloudElementTime']).replace(" ", "_")[:-6]
    latitude = accuTRMMData.createVariable('latitude', 'f8', ('lat',))
    longitude = accuTRMMData.createVariable('longitude', 'f8', ('lon',))
    rainFallacc = accuTRMMData.createVariable(
        'precipitation_Accumulation', 'f8', TRMMprecip)
    rainFallacc.units = 'mm'

    longitude[:] = LONTRMM[0, :]
    longitude.units = "degrees_east"
    longitude.long_name = "Longitude"

    latitude[:] = LATTRMM[:, 0]
    latitude.units = "degrees_north"
    latitude.long_name = "Latitude"

    rainFallacc[:] = accuPrecipRate[:]

    accuTRMMData.close()

    # do plot
    plotTitle = 'TRMM Accumulated [mm]'
    createPrecipPlot(np.squeeze(accuPrecipRate, axis=0),
                     LATTRMM[:, 0], LONTRMM[0, :], plotTitle, imgFilename)

    return
#******************************************************************


def plotAccuInTimeRange(starttime, endtime):
    '''
    Purpose::
            Create accumulated precip plot within a time range given using all CEs

    Input::
            starttime: a string representing the time to start the accumulations format yyyy-mm-dd_hh:mm:ss
            endtime: a string representing the time to end the accumulations format yyyy-mm-dd_hh:mm:ss

    Output::
            a netcdf file containing the accumulated precip for specified times
            a 2D matlablibplot

    '''

    os.chdir((MAINDIRECTORY + '/TRMMnetcdfCEs/'))

    imgFilename = ''
    firstPartName = ''
    firstTime = True

    fileList = []
    sTime = datetime.strptime(starttime.replace("_", " "), '%Y-%m-%d %H:%M:%S')
    eTime = datetime.strptime(endtime.replace("_", " "), '%Y-%m-%d %H:%M:%S')
    thisTime = sTime

    while thisTime <= eTime:
        fileList = list(
            filter(
                os.path.isfile,
                glob.glob(
                    ('TRMM' +
                     str(thisTime).replace(
                         " ",
                         "_") +
                        '*' +
                        '.nc'))))
        for fname in fileList:
            TRMMCEData = Dataset(fname, 'r', format='NETCDF4')
            precipRate = TRMMCEData.variables['precipitation_Accumulation'][:]
            lats = TRMMCEData.variables['latitude'][:]
            lons = TRMMCEData.variables['longitude'][:]
            LONTRMM, LATTRMM = np.meshgrid(lons, lats)
            nygrdTRMM = len(LATTRMM[:, 0])
            nxgrdTRMM = len(LONTRMM[0, :])
            precipRate = ma.masked_array(precipRate, mask=(precipRate < 0.0))
            TRMMCEData.close()

            if firstTime:
                accuPrecipRate = ma.zeros((precipRate.shape))
                firstTime = False

            accuPrecipRate += precipRate

        # increment the time
        thisTime += timedelta(hours=TRES)

    # create new netCDF file
    accuTRMMFile = MAINDIRECTORY + '/TRMMnetcdfCEs/accu' + \
        starttime + '-' + endtime + '.nc'
    print("accuTRMMFile %s" % accuTRMMFile)
    # write the file
    accuTRMMData = Dataset(accuTRMMFile, 'w', format='NETCDF4')
    accuTRMMData.description = 'Accumulated precipitation data'
    accuTRMMData.calendar = 'standard'
    accuTRMMData.conventions = 'COARDS'
    # dimensions
    accuTRMMData.createDimension('time', None)
    accuTRMMData.createDimension('lat', nygrdTRMM)
    accuTRMMData.createDimension('lon', nxgrdTRMM)

    # variables
    TRMMprecip = ('time', 'lat', 'lon',)
    times = accuTRMMData.createVariable('time', 'f8', ('time',))
    times.units = 'hours since ' + starttime[:-6]
    latitude = accuTRMMData.createVariable('latitude', 'f8', ('lat',))
    longitude = accuTRMMData.createVariable('longitude', 'f8', ('lon',))
    rainFallacc = accuTRMMData.createVariable(
        'precipitation_Accumulation', 'f8', TRMMprecip)
    rainFallacc.units = 'mm'

    longitude[:] = LONTRMM[0, :]
    longitude.units = "degrees_east"
    longitude.long_name = "Longitude"

    latitude[:] = LATTRMM[:, 0]
    latitude.units = "degrees_north"
    latitude.long_name = "Latitude"
    rainFallacc[:] = accuPrecipRate[:]

    accuTRMMData.close()

    # plot the stuff
    imgFilename = MAINDIRECTORY + '/images/accu' + starttime + '-' + endtime + '.gif'
    plotTitle = "TRMM Accumulated Precipitation [mm] " + \
        starttime + '-' + endtime
    createPrecipPlot(np.squeeze(accuPrecipRate, axis=0),
                     LATTRMM[:, 0], LONTRMM[0, :], plotTitle, imgFilename)

    return
#******************************************************************


def createPrecipPlot(dataset, lats, lons, plotTitle, imgFilename):
    '''
    Purpose::
            To create the actual plots for precip data only

    Input::
            dataset: a 2d numpy (lon,lat) dataset of the precip data to be plotted
            domainDict: a dictionary with the domain (lons and lats) details required
            plotTitle: a string representing the title for the plot
            imgFilename: a string representing the string (including path) of where the plot is to be showed

    Output::
            A 2D plot with precipitation using the NWS precipitation colormap (from matlplotlib Basemap)

    '''

    fig, ax = plt.subplots(
        1, facecolor='white', figsize=(
            8.5, 11.))  # , dpi=300)
    latmin = np.min(lats)
    latmax = np.max(lats)
    lonmin = np.min(lons)
    lonmax = np.max(lons)

    m = Basemap(
        projection='merc',
        llcrnrlon=lonmin,
        urcrnrlon=lonmax,
        llcrnrlat=latmin,
        urcrnrlat=latmax,
        resolution='l',
        ax=ax)
    m.drawcoastlines(linewidth=1)
    m.drawcountries(linewidth=0.75)
    # draw meridians
    meridians = np.arange(180., 360., 5.)
    m.drawmeridians(
        meridians,
        labels=[
            0,
            0,
            0,
            1],
        linewidth=0.75,
        fontsize=10)
    # draw parallels
    parallels = np.arange(0., 90, 5.)
    m.drawparallels(
        parallels,
        labels=[
            1,
            0,
            0,
            0],
        linewidth=0.75,
        fontsize=10)

    # projecting on to the correct map grid
    longitudes, latitudes = m.makegrid(lons.shape[0], lats.shape[0])
    x, y = m(longitudes, latitudes)

    # draw filled contours.
    clevs = [0, 1, 2.5, 5, 7.5, 10, 15, 20, 30, 40, 50,
             70, 100, 150, 200, 250, 300, 400, 500, 600, 750]

    # actually print the map
    cs = m.contourf(x, y, dataset, clevs, cmap=cmbm.s3pcpn)

    # add colorbar
    cbar = m.colorbar(cs, location='bottom', pad="10%")
    cbar.set_label('mm')

    plt.title(plotTitle)

    plt.savefig(imgFilename, facecolor=fig.get_facecolor(), transparent=True)
    return
#******************************************************************


def createTextFile(finalMCCList, identifier):
    '''
    Purpose::
            Create a text file with information about the MCS
            This function is expected to be especially of use regarding long term record checks

    Input::
            finalMCCList: a list of dictionaries representing a list of nodes representing a MCC
            identifier: an integer representing the type of list that has been entered...this is for creating file purposes
                    1 - MCCList; 2- MCSList

    Output::
            a user readable text file with all information about each MCS
            a user readable text file with the summary of the MCS

    Assumptions::
    '''

    durations = 0.0
    startTimes = []
    endTimes = []
    averagePropagationSpeed = 0.0
    speedCounter = 0
    maxArea = 0.0
    amax = 0.0
    avgMaxArea = []
    maxAreaCounter = 0.0
    maxAreaTime = ''
    eccentricity = 0.0
    firstTime = True
    matureFlag = True
    timeMCSMatures = ''
    maxCEprecipRate = 0.0
    minCEprecipRate = 0.0
    averageArea = 0.0
    averageAreaCounter = 0
    durationOfMatureMCC = 0
    avgMaxPrecipRate = 0.0
    avgMaxPrecipRateCounter = 0
    avgMinPrecipRate = 0.0
    avgMinPrecipRateCounter = 0
    CEspeed = 0.0
    MCSspeed = 0.0
    MCSspeedCounter = 0
    MCSPrecipTotal = 0.0
    avgMCSPrecipTotalCounter = 0
    bigPtotal = 0.0
    bigPtotalCounter = 0
    allPropagationSpeeds = []
    averageAreas = []
    areaAvg = 0.0
    avgPrecipTotal = 0.0
    avgPrecipTotalCounter = 0
    avgMaxMCSPrecipRate = 0.0
    avgMaxMCSPrecipRateCounter = 0
    avgMinMCSPrecipRate = 0.0
    avgMinMCSPrecipRateCounter = 0
    minMax = []
    avgPrecipArea = []
    location = []
    avgPrecipAreaPercent = 0.0
    precipArea = 0.0
    precipAreaPercent = 0.0
    precipPercent = []
    precipCounter = 0
    precipAreaAvg = 0.0
    minSpeed = 0.0
    maxSpeed = 0.0

    if identifier == 1:
        MCSUserFile = open(
            (MAINDIRECTORY + '/textFiles/MCCsUserFile.txt'), 'wb')
        MCSSummaryFile = open(
            (MAINDIRECTORY + '/textFiles/MCCSummary.txt'), 'wb')
        MCSPostFile = open(
            (MAINDIRECTORY + '/textFiles/MCCPostPrecessing.txt'), 'wb')

    if identifier == 2:
        MCSUserFile = open(
            (MAINDIRECTORY + '/textFiles/MCSsUserFile.txt'), 'wb')
        MCSSummaryFile = open(
            (MAINDIRECTORY + '/textFiles/MCSSummary.txt'), 'wb')
        MCSPostFile = open(
            (MAINDIRECTORY + '/textFiles/MCSPostPrecessing.txt'), 'wb')

    for eachPath in finalMCCList:
        eachPath.sort(key=lambda nodeID: (
            len(nodeID.split('C')[0]), nodeID.split('C')[0], nodeID.split('CE')[1]))
        MCSPostFile.write("\n %s" % eachPath)

        startTime = thisDict(eachPath[0])['cloudElementTime']
        endTime = thisDict(eachPath[-1])['cloudElementTime']
        duration = (endTime - startTime) + timedelta(hours=TRES)

        # convert datatime duration to seconds and add to the total for the
        # average duration of all MCS in finalMCCList
        durations += (duration.total_seconds())

        #durations += duration
        startTimes.append(startTime)
        endTimes.append(endTime)

        # get the precip info

        for eachNode in eachPath:

            thisNode = thisDict(eachNode)

            # set first time min "fake" values
            if firstTime:
                minCEprecipRate = thisNode['CETRMMmin']
                avgMinMCSPrecipRate += thisNode['CETRMMmin']
                firstTime = False

            # calculate the speed
            if thisNode['cloudElementArea'] >= OUTER_CLOUD_SHIELD_AREA:
                averagePropagationSpeed += findCESpeed(eachNode, eachPath)
                speedCounter += 1

            # Amax: find max area
            if thisNode['cloudElementArea'] > maxArea:
                maxArea = thisNode['cloudElementArea']
                maxAreaTime = str(thisNode['cloudElementTime'])
                eccentricity = thisNode['cloudElementEccentricity']
                location = thisNode['cloudElementCenter']

                # determine the time the feature matures
                if matureFlag:
                    timeMCSMatures = str(thisNode['cloudElementTime'])
                    matureFlag = False

            # find min and max precip rate
            if thisNode['CETRMMmin'] < minCEprecipRate:
                minCEprecipRate = thisNode['CETRMMmin']

            if thisNode['CETRMMmax'] > maxCEprecipRate:
                maxCEprecipRate = thisNode['CETRMMmax']

            # calculations for only the mature stage
            if thisNode['nodeMCSIdentifier'] == 'M':
                # calculate average area of the maturity feature only
                averageArea += thisNode['cloudElementArea']
                averageAreaCounter += 1
                durationOfMatureMCC += 1
                avgMaxPrecipRate += thisNode['CETRMMmax']
                avgMaxPrecipRateCounter += 1
                avgMinPrecipRate += thisNode['CETRMMmin']
                avgMinPrecipRateCounter += 1
                avgMaxMCSPrecipRate += thisNode['CETRMMmax']
                avgMaxMCSPrecipRateCounter += 1
                avgMinMCSPrecipRate += thisNode['CETRMMmin']
                avgMinMCSPrecipRateCounter += 1

                # the precip percentage (TRMM area/CE area)
                if thisNode['cloudElementArea'] >= 0.0 and thisNode['TRMMArea'] >= 0.0:
                    precipArea += thisNode['TRMMArea']
                    avgPrecipArea.append(thisNode['TRMMArea'])
                    avgPrecipAreaPercent += (
                        thisNode['TRMMArea'] / thisNode['cloudElementArea'])
                    precipPercent.append(
                        (thisNode['TRMMArea'] / thisNode['cloudElementArea']))
                    precipCounter += 1

                # system speed for only mature stage
                CEspeed = findCESpeed(eachNode, eachPath)
                if CEspeed > 0.0:
                    MCSspeed += CEspeed
                    MCSspeedCounter += 1

            # find accumulated precip
            if thisNode['cloudElementPrecipTotal'] > 0.0:
                MCSPrecipTotal += thisNode['cloudElementPrecipTotal']
                avgMCSPrecipTotalCounter += 1

        # A: calculate the average Area of the (mature) MCS
        if averageAreaCounter > 0:  # and averageAreaCounter > 0:
            averageArea /= averageAreaCounter
            averageAreas.append(averageArea)

        # v: MCS speed
        if MCSspeedCounter > 0:  # and MCSspeed > 0.0:
            MCSspeed /= MCSspeedCounter

        # smallP_max: calculate the average max precip rate (mm/h)
        if avgMaxMCSPrecipRateCounter > 0:  # and avgMaxPrecipRate > 0.0:
            avgMaxMCSPrecipRate /= avgMaxMCSPrecipRateCounter

        # smallP_min: calculate the average min precip rate (mm/h)
        if avgMinMCSPrecipRateCounter > 0:  # and avgMinPrecipRate > 0.0:
            avgMinMCSPrecipRate /= avgMinMCSPrecipRateCounter

        # smallP_avg: calculate the average precipitation (mm hr-1)
        if MCSPrecipTotal > 0.0:  # and avgMCSPrecipTotalCounter> 0:
            avgMCSPrecipTotal = MCSPrecipTotal / avgMCSPrecipTotalCounter
            avgPrecipTotal += avgMCSPrecipTotal
            avgPrecipTotalCounter += 1

        #smallP_total = MCSPrecipTotal
        # precip over the MCS lifetime prep for bigP_total
        if MCSPrecipTotal > 0.0:
            bigPtotal += MCSPrecipTotal
            bigPtotalCounter += 1

        if maxArea > 0.0:
            avgMaxArea.append(maxArea)
            maxAreaCounter += 1

        # verage precipate area precentage (TRMM/CE area)
        if precipCounter > 0:
            avgPrecipAreaPercent /= precipCounter
            precipArea /= precipCounter

        # write stuff to file
        MCSUserFile.write("\n\n\nStarttime is: %s " % (str(startTime)))
        MCSUserFile.write("\nEndtime is: %s " % (str(endTime)))
        MCSUserFile.write("\nLife duration is %s hrs" % (str(duration)))
        MCSUserFile.write("\nTime of maturity is %s " % (timeMCSMatures))
        MCSUserFile.write(
            "\nDuration mature stage is: %s " %
            durationOfMatureMCC * TRES)
        MCSUserFile.write("\nAverage area is: %.4f km^2 " % (averageArea))
        MCSUserFile.write("\nMax area is: %.4f km^2 " % (maxArea))
        MCSUserFile.write("\nMax area time is: %s " % (maxAreaTime))
        MCSUserFile.write(
            "\nEccentricity at max area is: %.4f " %
            (eccentricity))
        MCSUserFile.write(
            "\nCenter (lat,lon) at max area is: %.2f\t%.2f" %
            (location[0], location[1]))
        MCSUserFile.write("\nPropagation speed is %.4f " % (MCSspeed))
        MCSUserFile.write(
            "\nMCS minimum preicip rate is %.4f mmh^-1" %
            (minCEprecipRate))
        MCSUserFile.write(
            "\nMCS maximum preicip rate is %.4f mmh^-1" %
            (maxCEprecipRate))
        MCSUserFile.write(
            "\nTotal precipitation during MCS is %.4f mm/lifetime" %
            (MCSPrecipTotal))
        MCSUserFile.write(
            "\nAverage MCS precipitation is %.4f mm" %
            (avgMCSPrecipTotal))
        MCSUserFile.write(
            "\nAverage MCS maximum precipitation is %.4f mmh^-1" %
            (avgMaxMCSPrecipRate))
        MCSUserFile.write(
            "\nAverage MCS minimum precipitation is %.4f mmh^-1" %
            (avgMinMCSPrecipRate))
        MCSUserFile.write(
            "\nAverage precipitation area is %.4f km^2 " %
            (precipArea))
        MCSUserFile.write(
            "\nPrecipitation area percentage of mature system %.4f percent " %
            (avgPrecipAreaPercent * 100))

        # append stuff to lists for the summary file
        if MCSspeed > 0.0:
            allPropagationSpeeds.append(MCSspeed)
            averagePropagationSpeed += MCSspeed
            speedCounter += 1

        # reset vars for next MCS in list
        aaveragePropagationSpeed = 0.0
        averageArea = 0.0
        averageAreaCounter = 0
        durationOfMatureMCC = 0
        MCSspeed = 0.0
        MCSspeedCounter = 0
        MCSPrecipTotal = 0.0
        avgMaxMCSPrecipRate = 0.0
        avgMaxMCSPrecipRateCounter = 0
        avgMinMCSPrecipRate = 0.0
        avgMinMCSPrecipRateCounter = 0
        firstTime = True
        matureFlag = True
        avgMCSPrecipTotalCounter = 0
        avgPrecipAreaPercent = 0.0
        precipArea = 0.0
        precipCounter = 0
        maxArea = 0.0
        maxAreaTime = ''
        eccentricity = 0.0
        timeMCSMatures = ''
        maxCEprecipRate = 0.0
        minCEprecipRate = 0.0
        location = []

    # LD: average duration
    if len(finalMCCList) > 1:
        durations /= len(finalMCCList)
        durations /= 3600.0  # convert to hours

        # A: average area
        areaAvg = sum(averageAreas) / len(finalMCCList)
    # create histogram plot here
    if len(averageAreas) > 1:
        imgFilename = MAINDIRECTORY + '/images/averageAreas.gif'
        plotter.draw_histogram(
            averageAreas, [
                "Average Area [km^2]", "Area [km^2]"], imgFilename, 10)

    # Amax: average maximum area
    if maxAreaCounter > 0.0:  # and avgMaxArea > 0.0 :
        amax = sum(avgMaxArea) / maxAreaCounter
        # create histogram plot here
        if len(avgMaxArea) > 1:
            imgFilename = MAINDIRECTORY + '/images/avgMaxArea.gif'
            plotter.draw_histogram(
                avgMaxArea, [
                    "Maximum Area [km^2]", "Area [km^2]"], imgFilename, 10)

    # v_avg: calculate the average propagation speed
    if speedCounter > 0:  # and averagePropagationSpeed > 0.0
        averagePropagationSpeed /= speedCounter

    # bigP_min: calculate the min rate in mature system
    if avgMinPrecipRate > 0.0:  # and avgMinPrecipRateCounter > 0.0:
        avgMinPrecipRate /= avgMinPrecipRateCounter

    # bigP_max: calculate the max rate in mature system
    if avgMinPrecipRateCounter > 0.0:  # and avgMaxPrecipRate >  0.0:
        avgMaxPrecipRate /= avgMaxPrecipRateCounter

    # bigP_avg: average total preicip rate mm/hr
    if avgPrecipTotalCounter > 0.0:  # and avgPrecipTotal > 0.0:
        avgPrecipTotal /= avgPrecipTotalCounter

    # bigP_total: total precip rate mm/LD
    if bigPtotalCounter > 0.0:  # and bigPtotal > 0.0:
        bigPtotal /= bigPtotalCounter

    # precipitation area percentage
    if len(precipPercent) > 0:
        precipAreaPercent = (sum(precipPercent) / len(precipPercent)) * 100.0

    # average precipitation area
    if len(avgPrecipArea) > 0:
        precipAreaAvg = sum(avgPrecipArea) / len(avgPrecipArea)
        if len(avgPrecipArea) > 1:
            imgFilename = MAINDIRECTORY + '/images/avgPrecipArea.gif'
            plotter.draw_histogram(
                avgPrecipArea, [
                    "Average Rainfall Area [km^2]", "Area [km^2]"], imgFilename, 10)

    sTime = str(averageTime(startTimes))
    eTime = str(averageTime(endTimes))
    if len(allPropagationSpeeds) > 1:
        maxSpeed = max(allPropagationSpeeds)
        minSpeed = min(allPropagationSpeeds)

    # write stuff to the summary file
    MCSSummaryFile.write("\nNumber of features is %d " % (len(finalMCCList)))
    MCSSummaryFile.write("\nAverage duration is %.4f hrs " % (durations))
    MCSSummaryFile.write("\nAverage startTime is %s " % (sTime[-8:]))
    MCSSummaryFile.write("\nAverage endTime is %s " % (eTime[-8:]))
    MCSSummaryFile.write("\nAverage size is %.4f km^2 " % (areaAvg))
    MCSSummaryFile.write(
        "\nAverage precipitation area is %.4f km^2 " %
        (precipAreaAvg))
    MCSSummaryFile.write("\nAverage maximum size is %.4f km^2 " % (amax))
    MCSSummaryFile.write(
        "\nAverage propagation speed is %.4f ms^-1" %
        (averagePropagationSpeed))
    MCSSummaryFile.write(
        "\nMaximum propagation speed is %.4f ms^-1 " %
        (maxSpeed))
    MCSSummaryFile.write(
        "\nMinimum propagation speed is %.4f ms^-1 " %
        (minSpeed))
    MCSSummaryFile.write(
        "\nAverage minimum precipitation rate is %.4f mmh^-1" %
        (avgMinPrecipRate))
    MCSSummaryFile.write(
        "\nAverage maximum precipitation rate is %.4f mm h^-1" %
        (avgMaxPrecipRate))
    MCSSummaryFile.write(
        "\nAverage precipitation is %.4f mm h^-1 " %
        (avgPrecipTotal))
    MCSSummaryFile.write(
        "\nAverage total precipitation during MCSs is %.4f mm/LD " %
        (bigPtotal))
    MCSSummaryFile.write(
        "\nAverage precipitation area percentage is %.4f percent " %
        (precipAreaPercent))

    MCSUserFile.close
    MCSSummaryFile.close
    MCSPostFile.close
    return
#******************************************************************
#			PLOTTING UTIL SCRIPTS
#******************************************************************


def to_percent(y, position):
    '''
    Purpose::
            Utility script for generating the y-axis for plots
    '''
    return (str(100 * y) + '%')
#******************************************************************


def colorbar_index(ncolors, nlabels, cmap):
    '''
    Purpose::
            Utility script for crating a colorbar
            Taken from http://stackoverflow.com/questions/18704353/correcting-matplotlib-colorbar-ticks

    '''
    cmap = cmap_discretize(cmap, ncolors)
    mappable = cm.ScalarMappable(cmap=cmap)
    mappable.set_array([])
    mappable.set_clim(-0.5, ncolors + 0.5)
    colorbar = plt.colorbar(mappable)  # , orientation='horizontal')
    colorbar.set_ticks(np.linspace(0, ncolors, ncolors))
    colorbar.set_ticklabels(nlabels)
    return
#******************************************************************


def cmap_discretize(cmap, N):
    '''
    Taken from: http://stackoverflow.com/questions/18704353/correcting-matplotlib-colorbar-ticks
    http://wiki.scipy.org/Cookbook/Matplotlib/ColormapTransformations
    Return a discrete colormap from the continuous colormap cmap.

            cmap: colormap instance, eg. cm.jet.
            N: number of colors.

    Example
            x = resize(arange(100), (5,100))
            djet = cmap_discretize(cm.jet, 5)
            imshow(x, cmap=djet)
    '''

    if isinstance(cmap, str):
        cmap = plt.get_cmap(cmap)
    colors_i = np.concatenate((np.linspace(0, 1., N), (0., 0., 0., 0.)))
    colors_rgba = cmap(colors_i)
    indices = np.linspace(0, 1., N + 1)
    cdict = {}
    for ki, key in enumerate(('red', 'green', 'blue')):
        cdict[key] = [(indices[i], colors_rgba[i - 1, ki], colors_rgba[i, ki])
                      for i in range(N + 1)]
    # Return colormap object.
    return mcolors.LinearSegmentedColormap(cmap.name + "_%d" % N, cdict, 1024)
#******************************************************************
# def preprocessingMERG(MERGdirname):
# 	'''
# 	Purpose::
# 		Utility script for unzipping and converting the merg*.Z files from Mirador to
# 		NETCDF format. The files end up in a folder called mergNETCDF in the directory
# 		where the raw MERG data is
# 		NOTE: VERY RAW AND DIRTY
# 	Input::
# 		Directory to the location of the raw MERG files, preferably zipped
# 	Output::
# 	   none
# 	Assumptions::
# 	   1 GrADS (http://www.iges.org/grads/gadoc/) and lats4D (http://opengrads.org/doc/scripts/lats4d/)
# 		 have been installed on the system and the user can access
# 	   2 User can write files in location where script is being called
# 	   3 the files havent been unzipped
# 	'''

# 	os.chdir((MERGdirname+'/'))
# 	imgFilename = ''

# 	#Just incase the X11 server is giving problems
# 	subprocess.call('export DISPLAY=:0.0', shell=True)

# 	for files in glob.glob("*-pixel"):
# 	#for files in glob.glob("*.Z"):
# 		fname = os.path.splitext(files)[0]

# 		#unzip it
# 		bash_cmd = 'gunzip ' + files
# 		subprocess.call(bash_cmd, shell=True)

# 		#determine the time from the filename
# 		ftime = re.search('\_(.*)\_',fname).group(1)

# 		yy = ftime[0:4]
# 		mm = ftime[4:6]
# 		day = ftime[6:8]
# 		hr = ftime [8:10]

# 		#TODO: must be something more efficient!

# 		if mm=='01':
# 			mth = 'Jan'
# 		if mm == '02':
# 			mth = 'Feb'
# 		if mm == '03':
# 			mth = 'Mar'
# 		if mm == '04':
# 			mth = 'Apr'
# 		if mm == '05':
# 			mth = 'May'
# 		if mm == '06':
# 			mth = 'Jun'
# 		if mm == '07':
# 			mth = 'Jul'
# 		if mm == '08':
# 			mth = 'Aug'
# 		if mm == '09':
# 			mth = 'Sep'
# 		if mm == '10':
# 			mth = 'Oct'
# 		if mm == '11':
# 			mth = 'Nov'
# 		if mm == '12':
# 			mth = 'Dec'


# 		subprocess.call('rm merg.ctl', shell=True)
# 		subprocess.call('touch merg.ctl', shell=True)
# 		replaceExpDset = 'echo DSET ' + fname +' >> merg.ctl'
# 		replaceExpTdef = 'echo TDEF 99999 LINEAR '+hr+'z'+day+mth+yy +' 30mn' +' >> merg.ctl'
# 		subprocess.call(replaceExpDset, shell=True)
# 		subprocess.call('echo "OPTIONS yrev little_endian template" >> merg.ctl', shell=True)
# 		subprocess.call('echo "UNDEF  330" >> merg.ctl', shell=True)
# 		subprocess.call('echo "TITLE  globally merged IR data" >> merg.ctl', shell=True)
# 		subprocess.call('echo "XDEF 9896 LINEAR   0.0182 0.036378335" >> merg.ctl', shell=True)
# 		subprocess.call('echo "YDEF 3298 LINEAR   -59.982 0.036383683" >> merg.ctl', shell=True)
# 		subprocess.call('echo "ZDEF   01 LEVELS 1" >> merg.ctl', shell=True)
# 		subprocess.call(replaceExpTdef, shell=True)
# 		subprocess.call('echo "VARS 1" >> merg.ctl', shell=True)
# 		subprocess.call('echo "ch4  1  -1,40,1,-1 IR BT  (add  "75" to this value)" >> merg.ctl', shell=True)
# 		subprocess.call('echo "ENDVARS" >> merg.ctl', shell=True)

# 		#generate the lats4D command for GrADS
# 		lats4D = 'lats4d -v -q -lat '+LATMIN + ' ' +LATMAX +' -lon ' +LONMIN +' ' +LONMAX +' -time '+hr+'Z'+day+mth+yy + ' -func @+75 ' + '-i merg.ctl' + ' -o ' + fname

# 		#lats4D = 'lats4d -v -q -lat -40 -15 -lon 10 40 -time '+hr+'Z'+day+mth+yy + ' -func @+75 ' + '-i merg.ctl' + ' -o ' + fname
# 		#lats4D = 'lats4d -v -q -lat -5 40 -lon -90 60 -func @+75 ' + '-i merg.ctl' + ' -o ' + fname

# 		gradscmd = 'grads -blc ' + '\'' +lats4D + '\''
# 		#run grads and lats4d command
# 		subprocess.call(gradscmd, shell=True)
# 		imgFilename = hr+'Z'+day+mth+yy+'.gif'
# 		tempMaskedImages(imgFilename)

# 	#when all the files have benn converted, mv the netcdf files
# 	subprocess.call('mkdir mergNETCDF', shell=True)
# 	subprocess.call('mv *.nc mergNETCDF', shell=True)
# 	#mv all images
# 	subprocess.call('mkdir mergImgs', shell=True)
# 	subprocess.call('mv *.gif mergImgs', shell=True)
# 	return