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