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