import java.lang.Float as Float import jarray, inspect def searchNeighbors(gridValues, index, currentPixel, radius, assigned, badPixels, numAreaLines, numAreaElements, domainSet,group): """ searchNeigbors: accept pixel, find neighbors within a certain radius, find out those neighbors are in the list of badPixelss group the pixels. gridValues: Should I make this a global? This is a rank 2 array which returns line/elems values based on indices. index: The index within the badPixels and assigned array of the currentPixel currentPixel: The actual indice of the bad pixel within the area data. (badPixels(index)) radius: Numer of pixels to search around current pixel for clustering assigned: The array which contains the grouping numbers. 0: unassigned, -1: lat/lon nan, -2: isolated 1-n: cluster groups badPixels: the list of pixels which fell within a certain criteria. numAreaLines, numAreaElements: from the area metadata. The number of lines/elems in this area domainSet: 2D domain set retrieved from the area. group: current group number Not sure if I am doing this right, it can't be run without retrieving many of the inputs in the code previous. """ useLine=[] useElem=[] # assign the current grid point to the first value in the array to check for matches useLine.append(gridValues[1][index]) useElem.append(gridValues[0][index]) # get the corresponding grid values and check neighbors within radius for incrementLine in range(-radius,radius+1): for incrementElem in range(-radius, radius+1): # no need to check current pixel (where incrementLine + incrementElem = 0) if (incrementElem == 0 and incrementLine == 0): pass else: checkLine = useLine[0] + incrementLine checkElem = useElem[0] + incrementElem #make sure list is not empty and not in first or last row/col if (checkLine > -1 and int(checkLine) < numAreaLines and checkElem > -1 and int(checkElem) < numAreaElements): useLine.append(checkLine) useElem.append(checkElem) matches = domainSet.valueToIndex(domainSet.gridToValue([useElem, useLine])) # Convert the matching grid values back to indices, so that the surrounding pixels can be checked for existence in the indices found. # Do this by taking the line/elem pairs and converting to an array of R^2 values (using gridToValue) to 1-D indices (using valueToIndex) if (len(useLine) > 1): matches = domainSet.valueToIndex(domainSet.gridToValue([useElem[1:], useLine[1:]])) # loop through matches array to group data. usedMatch=[] for nMatch in matches: try: matchedPixel = badPixels.index(nMatch) if ((assigned[matchedPixel] == 0) or (assigned[matchedPixel] == -2)): #assigned[index] = assigned[matchedPixel] # otherwise, assign a new group to both the currentPixel and the match #else: #if (assigned[index] < 0): #incrementGroup = True assigned[index] = group assigned[matchedPixel] = group # if this match was used, will check neighbors of this pixel next usedMatch.append(nMatch) #print 'Build', usedMatch, nMatch except (ValueError, IndexError): pass return assigned, usedMatch #************************************************************************************** dataSet = 'GOES-12' imageType = 'Validation' desc = getDescriptor(dataSet, imageType) # test with this data addeParms = dict( debug=True, server='localhost:8112', coordinateSystem=CoordinateSystems.LATLON, place=Places.CENTER, location=(-16.6, -62.3), size=(20,16), mag=(1,1), ) metaData, band2Temp=getADDEImage(dataset=dataSet, descriptor=desc, band=2, unit='TEMP', **addeParms) #Find out how many data points are in the line (number of elements in the satellite image) numAreaLines = metaData['directory-block'][8] numAreaElements = metaData['directory-block'][9] print band2Temp # find the indices where the values match the criteria badPixels=findOutsideRange(band2Temp,160,265) numBadPixels = badPixels.__len__() # get domain of this data set band2Domain = band2Temp.getDomainSet() # use methods in Gridded1Dset/Subclass Linear1Dset to convert from indices to grid values band2IndexSet = band2Domain.indexToValue(badPixels) # gridValues is an array of grid values which correspond to the indices array gridValues = band2Domain.valueToGrid(band2IndexSet) # Set up an array of zeros the same size as the number of indices found (zero will be not grouped) assigned = jarray.zeros(numBadPixels,'i') # get the lat/lon values for the entire domain of the data. fieldLatLonValues = GridUtil.getLatLon(band2Domain) # create the Area Field from the flat field (to be used later when calculating area affected by bad pixels band2AreaField = createAreaField(band2Temp) #loop through the indices, check for valid lat/lon values, group neighboring pixels group = 1 # set the radius around a pixel for grouping to 1 pixel radius = 1 for index, currentPixel in enumerate(badPixels): if (assigned[index] == 0): # check for valid lat/lon values if (Float.isNaN(fieldLatLonValues[0][currentPixel]) or Float.isNaN(fieldLatLonValues[1][currentPixel])): assigned[index]=-1 # assign a group. start with isolated pixel assignment and change if neighboring pixels are found if (assigned[index] == 0): assigned[index] = -2 # find the pixels around this current pixel assigned, matchedList = searchNeighbors(gridValues, index, currentPixel, radius, assigned, badPixels, numAreaLines, numAreaElements, band2Domain,group) if (len(matchedList) > 0): incrementGroup = True # loop through list of matches, checking neighbors and grouping while (len(matchedList) > 0): checkThis = matchedList.pop(0) # force checkList to an integer. When pop is used, value is returned as a list] if (type(checkThis).__name__ == 'list'): checkThis = checkThis[0] checkIndex = badPixels.index(checkThis) #print 'Checking', index, currentPixel, checkThis, checkIndex secondaryMatches=[] assigned, secondaryMatches = searchNeighbors(gridValues, checkIndex, checkThis, radius, assigned, badPixels, numAreaLines, numAreaElements, band2Domain,group) if (len(secondaryMatches) > 0): for extraMatches in secondaryMatches: matchedList.append(extraMatches) #print 'In matching loop', index, currentPixel, matchedList # if this group number has been used for this set, increment the group number if (incrementGroup == True): group = group+1 incrementGroup = False print assigned print badPixels # Now loop through array and assigned group lists to a dictionary. Also, replace original data with group numbers when exists. grouping = {} for index, group in enumerate(assigned): if grouping.has_key(group): grouping[group].append(badPixels[index]) else: grouping.setifabsent(group, [badPixels[index]]) # Now compute the area coverage of bad pixels... areaDict = {} areaDict['total'] = computeSum(band2AreaField, badPixels) for key in grouping: a = computeSum(band2AreaField, grouping[key]) areaDict.setifabsent(key, a) # now write to a file