def noUnit(field):
  """ remove the units from a grid """
  import visad
  from visad import CommonUnit
  newunit = CommonUnit.promiscuous
  rt = GridUtil.getParamType(field).getRealComponents()[0]
  newType = Util.makeRealType(rt.getName(), visad.CommonUnit.promiscuous)
  return GridUtil.setParamType(field, newType,0)

def mycombineRGB(red, green, blue):
    """Three Color (RGB) Image (Auto-scale) formula."""
    global uniqueID
    uniqueID += 1
    red = GridUtil.setParamType(red, makeRealType("redimage%d" % uniqueID), 0)
    green = GridUtil.setParamType(green, makeRealType("greenimage%d" % uniqueID), 0)
    blue = GridUtil.setParamType(blue, makeRealType("blueimage%d" % uniqueID), 0)
    return DerivedGridFactory.combineGrids([red, green, blue], 1)

def holeFilter(sdataset, user_brkpoint1=0, user_brkpoint2=1, user_stretchval='Contrast', user_britlo=0, user_brithi=255):
    """ hole filter from McIDAS-X - searches for missing data and fills the holes
          using the surrounding element values
        brkpoint1 - low end breakpoint value (default = minimum sdataset value)
        brkpoint2 - high end breakpoint value (default = maximum sdataset value)
        user_britlo    - minimum brightness value for the calibration
        user_brithi    - maximum brightness value for the calibration
    """
    
    data = sdataset.clone()
    brkpoint1 = int(user_brkpoint1)
    brkpoint2 = int(user_brkpoint2)
    stretch = user_stretchval
    britlo = int(user_britlo)
    brithi = int(user_brithi)
    
    for t in xrange(data.getDomainSet().getLength()):
        rangeObject = data.getSample(t)
        vals = rangeObject.getFloats(0)
        in_low = min(vals[0])
        in_hi = max(vals[0])
        minVal = min([brkpoint1, brkpoint2])
        maxVal = max([brkpoint1, brkpoint2])
        domain = GridUtil.getSpatialDomain(rangeObject)
        [element_size, line_size] = domain.getLengths()
        
        for i in range(line_size):
            for j in range(element_size)[1:-1]:
                curVal = vals[0][i * element_size + j]
                """ search line for bad values """
                if curVal >= minVal and curVal <= maxVal:
                    """ look for the next good value """
                    doFill = 0
                    for k in range(element_size)[j:]:
                        nextVal = vals[0][i * element_size + k]
                        if nextVal < minVal or nextVal > maxVal:
                            doFill = 1
                            break
                            
                    if doFill == 1:
                        for fill in range(element_size)[j:k]:
                            vals[0][i * element_size + fill] = (vals[0][i * element_size + j - 1] + vals[0][i * element_size + k]) / 2
                            
        for i in range(line_size):
            for j in range(element_size):
                vals[0][i * element_size + j] = scaleOutsideVal(vals[0][i * element_size + j], britlo, brithi)
                   
        filt_low = int(min(vals[0]))
        filt_hi = int(max(vals[0]))
        
        if stretch == 'Contrast':
            lookup = contrast(filt_low, filt_hi, britlo, brithi, filt_low, filt_hi)
        elif stretch == 'Histogram':
            h = makeHistogram(vals, element_size, line_size, filt_low, brithi - britlo)
            lookup = histoStretch(filt_low, filt_hi, in_low, in_hi, h)
            
        vals = modify(vals, element_size, line_size, filt_low, lookup)
        rangeObject.setSamples(vals)
        
    return data

 def wrapper(*args, **kwargs):
     wrappedArgs = []
     returnFlatField = False
     for i, arg in enumerate(args):
         # print '%s: %s' % (i, arg)
         if isinstance(arg, FlatField):
             returnFlatField = True
             arg = makeFlatFieldSequence([arg])
         wrappedArgs.append(arg)
     # print [type(a) for a in wrappedArgs]
     wrappedKwargs = {}
     for keyword in kwargs:
         keywordValue = kwargs[keyword]
         if isinstance(keywordValue, FlatField):
             returnFlatField = True # ??? not sure about kwarg case
             keywordValue = makeFlatFieldSequence([keywordValue])
         wrappedKwargs[keyword] = keywordValue
     # print [type(wrappedKwargs[a]) for a in wrappedKwargs]
     result = func(*wrappedArgs, **wrappedKwargs)
     # print 'result type=%s' % (type(result))
     if GridUtil.isTimeSequence(result) and len(result) == 1 and returnFlatField:
         # print 'attempting conversion...'
         result = result.getSample(0)
     # else:
         # print 'not a time sequence! isseq=%s seqtype=%s' % (GridUtil.isSequence(result), GridUtil.getSequenceType(result))
     # print 'returning type=%s' % (type(result))
     return result

def highPass2DFilter(sdataset, user_stretchval='Contrast', user_britlo=0, user_brithi=255):
    """Two dimensional high pass filter from McIDAS-X.
    
    Equation for each sdataset element is 
    (sdataset - (sample average) + (sample midpoint))
    
    Args:
        user_britlo: Minimum brightness value for the calibration.
        user_brithi: Maximum brightness value for the calibration.
    """
    newData = sdataset.clone()
    stretch = user_stretchval
    britlo = int(user_britlo)
    brithi = int(user_brithi)
    
    for t in xrange(newData.getDomainSet().getLength()):
        rangeObject = newData.getSample(t)
        vals = rangeObject.getFloats(0)
        in_hi = max(vals[0])
        in_low = min(vals[0])
        domain = GridUtil.getSpatialDomain(rangeObject)
        [element_size, line_size] = domain.getLengths()
        
        # first and last 2 lines of the image do not change
        firstLine = vals[0][0:element_size]
        last2Lines = vals[0][(line_size - 2) * element_size:line_size * element_size]
        
        # do the filter using 3 lines at a time
        for i in range(line_size)[:-3]:
            for j in range(element_size)[1:-1]:
                midValue = vals[0][(i + 1) * element_size + j]
                
                val = (vals[0][(i + 2) * element_size + j] + vals[0][i * element_size + j] + \
                      vals[0][(i + 1) * element_size + j + 1] + vals[0][(i + 1) * element_size + j - 1]) - \
                      4 * midValue
                
                if (val < midValue):
                    vals[0][i * element_size + j] = scaleOutsideVal(midValue - val, britlo, brithi)
                else:
                    vals[0][i * element_size + j] = 0
                    
        vals[0][0:line_size * element_size] = firstLine + vals[0][0:(line_size - 3) * element_size] + last2Lines
        
        for i in range(line_size):
            for j in range(element_size):
                vals[0][i * element_size + j] = scaleOutsideVal(vals[0][i * element_size + j], britlo, brithi)
                
        filt_low = int(min(vals[0]))
        filt_hi = int(max(vals[0]))
        
        if stretch == 'Contrast':
            lookup = contrast(filt_low, filt_hi, britlo, brithi, filt_low, filt_hi)
        elif stretch == 'Histogram':
            h = makeHistogram(vals, element_size, line_size, filt_low, brithi - britlo)
            lookup = histoStretch(filt_low, filt_hi, in_low, in_hi, h)
            
        vals = modify(vals, element_size, line_size, filt_low, lookup)
        rangeObject.setSamples(vals)
        
    return newData

def getSliceAtLevel(fieldimpl, level) :
  """Extract a 2D horizontal slice from a 3D grid at "Level."
     level is a real number; must be appropriate for the grid.
     param fieldimpl is a grid which may have one or more time steps.  """
  level = float(level)
  ff = GridUtil.sliceAtLevel(fieldimpl, level)
  return ff

def shotFilter(sdataset, user_bline='Default', user_eline='Default', user_pdiff=15, user_stretchval='Contrast', user_britlo=0, user_brithi=255):
    """ shot noise filter from McIDAS-X
        bline - beginning line in the source image to clean (default=first line)
        eline - ending line in the source image to clean (default = last line)
        pdiff - maximum percentage of the product range to allow before a new value for the pixel is derived using the
                average of two adjacent pixels
        user_britlo    - minimum brightness value for the calibration
        user_brithi    - maximum brightness value for the calibration
    """
    newData = sdataset.clone()
    bline = user_bline
    eline = user_eline
    britlo = int(user_britlo)
    brithi = int(user_brithi)
       
    if bline != 'Default':
        bline = int(bline)
    else:
        bline = 0
        
    if eline != 'Default':
        eline = int(eline)
      
    filter_diff = int(user_pdiff)
    stretch = user_stretchval
    
    for t in range(newData.getDomainSet().getLength()):
        rangeObject = newData.getSample(t)
        vals = rangeObject.getFloats(0)
        in_hi = int(max(vals[0]))
        in_low = int(min(vals[0]))
        """ the next four lines are to make sure the point_diff value is floored """
        p = (in_hi - in_low + 1) * (filter_diff / 100.0)
        a = field((p,))
        b = a.floor().getValues(0)
        point_diff = b[0]
        
        domain = GridUtil.getSpatialDomain(rangeObject)
        [element_size, line_size] = domain.getLengths()
        if (eline == 'Default'):
            eline = line_size
            
        vals = shotMain(vals, bline, eline, element_size, line_size, point_diff)
        
        for i in range(line_size):
            for j in range(element_size):
                vals[0][i * element_size + j] = scaleOutsideVal(vals[0][i * element_size + j], britlo, brithi)
                
        filt_low = int(min(vals[0]))
        filt_hi = int(max(vals[0]))
        if stretch == 'Contrast':
            lookup = contrast(filt_low, filt_hi, britlo, brithi, filt_low, filt_hi)
        elif stretch == 'Histogram':
            h = makeHistogram(vals, element_size, line_size, filt_low, brithi - britlo)
            lookup = histoStretch(filt_low, filt_hi, britlo, brithi, h)
            
        vals = modify(vals, element_size, line_size, filt_low, lookup)
        rangeObject.setSamples(vals)
        
    return newData

def getAltitude(z):
   """ Change units from geopotential meters to meters """
   import ucar.visad.quantities.GeopotentialAltitude as ga
   import ucar.visad.quantities.Gravity as gr
   zUnit = GridUtil.getParamType(z).getRealComponents()[0].getDefaultUnit()
   if zUnit.equals(ga.getGeopotentialMeter()):
      z = z.divide(gr.newReal())
   return z

def printVals(sdataset):
    data = sdataset.clone()
    
    for t in xrange(data.getDomainSet().getLength()):
        rangeObj = data.getSample(t)
        vals = rangeObj.getFloats(0)
        domain = GridUtil.getSpatialDomain(rangeObj)
        [element_size, line_size] = domain.getLengths()
        
        for i in range(line_size):
            for j in range(element_size):
                print i, j, vals[0][i * element_size + j]

def applyToRange(function,data):
    """ Apply the function name to each timestep of the data """
    newData = data.clone()
    f = function +'(rangeValue)'
    if (GridUtil.isTimeSequence(newData)):
        for t in range(newData.getDomainSet().getLength()):
            rangeValue = newData.getSample(t)
            result = eval(f)
            newData.setSample(t,result,0)
    else:
        rangeValue = newData
        newData = eval(f)
    return newData

def oldaverageOverTime(field,makeTimes = 0):
    """@deprecated Average the values in each time step
    If makeTimes is true (1) then we return a field mapping all of the times
    to the average. Else we just return the average """
    if (GridUtil.isTimeSequence(field)==0):
        return field;
    cnt = 0;
    domainSet = field.getDomainSet()
    current = None;
    for t in range(domainSet.getLength()):
        cnt=cnt+1
        rangeValue = field.getSample(t)
        if(current is None):
            current = rangeValue.clone();
        else:
            current = current+rangeValue;
    if(cnt == 0):
        return None;
    current = current/cnt;
    if(makeTimes):
        return Util.makeTimeField(current, GridUtil.getDateTimeList(field))
    return current

def getSliceAtAltitude(fieldimpl, alt, unit="m") :
  """ Extract a 2D horizontal slice from a 3D grid at the given altitude;
      level is a real number; if unit is supplied, it must
      be compatible with meters (ft, fathoms, etc)
      param fieldimpl is a grid which may have
      one or more time steps.  """
  #import methods from
  from visad import RealType
  from visad import Real
  alt = float(alt)
  unit = Util.parseUnit(unit)
  altitude = Real(RealType.Altitude, alt, unit)
  ff = GridUtil.sliceAtLevel(fieldimpl, altitude)
  return ff

def printValueDiff(sdataset1, sdataset2):
    data1 = sdataset1.clone()
    data2 = sdataset2.clone()
    
    for t in xrange(data1.getDomainSet().getLength()):
        rangeObj1 = data1.getSample(t)
        rangeObj2 = data2.getSample(t)
        vals1 = rangeObj1.getFloats(0)
        """ vals2 = rangeObj1.getFloats(0)  TODO: bug? """
        vals2 = rangeObj2.getFloats(0)
        domain = GridUtil.getSpatialDomain(rangeObj1)
        [element_size, line_size] = domain.getLengths()
        
        for i in xrange(line_size):
            for j in xrange(element_size):
                print i, j, vals1[0][i * element_size + j] - vals2[0][i * element_size + j]

def printValueDiff(sdataset1, sdataset2):
    """Debug function that prints value differences between two datasets."""
    data1 = sdataset1.clone()
    data2 = sdataset2.clone()
    
    for t in xrange(data1.getDomainSet().getLength()):
        rangeObj1 = data1.getSample(t)
        rangeObj2 = data2.getSample(t)
        vals1 = rangeObj1.getFloats(0)
        # vals2 = rangeObj1.getFloats(0)  TODO: bug?
        vals2 = rangeObj2.getFloats(0)
        domain = GridUtil.getSpatialDomain(rangeObj1)
        [element_size, line_size] = domain.getLengths()
        
        for i in xrange(line_size):
            for j in xrange(element_size):
                print i, j, vals1[0][i * element_size + j] - vals2[0][i * element_size + j]

 def wrapper(*args):
     wrappedArgs = []
     for i, arg in enumerate(args):
         # print '%s: %s' % (i, arg)
         if isinstance(arg, SingleBandedImage):
             arg = ImageSequenceImpl([arg])
         wrappedArgs.append(arg)
     # print [type(a) for a in wrappedArgs]
     result = func(*wrappedArgs)
     # print 'result type=%s' % (type(result))
     if GridUtil.isTimeSequence(result) and len(result) == 1:
         # print 'attempting conversion...'
         result = result.getImage(0)
     # else:
         # print 'not a time sequence! isseq=%s seqtype=%s' % (GridUtil.isSequence(result), GridUtil.getSequenceType(result))
     # print 'returning type=%s' % (type(result))
     return result

def applyToRangeValues(function,data):
    """ Apply the function name to each value in each timestep of the data """
    newData = data.clone()
    f = function +'(values,step=step,rangeObject=rangeObject,field=field)'
    step=0
    if (GridUtil.isTimeSequence(newData)):
        for t in range(newData.getDomainSet().getLength()):
            rangeObject = newData.getSample(t)
            values = rangeObject.getFloats(0)
            values = eval(f)
            rangeObject.setSamples(values,1)
            step= step+1
    else:
        rangeObject = newData
        values = rangeObject.getFloats(0)
        values = eval(f)
        rangeObject.setSamples(values,1)
    return newData

def gradientFilter(sdataset, user_stretchval='Contrast', user_britlo=0, user_brithi=255):
    """Gradient filter from McIDAS-X.
    
    Args:
        user_britlo: Minimum brightness value for the calibration.
        user_brithi: Maximum brightness value for the calibration.
    """
    newData = sdataset.clone()
    stretch = user_stretchval
    britlo = int(user_britlo)
    brithi = int(user_brithi)
    
    for t in range(newData.getDomainSet().getLength()):
        rangeObject = newData.getSample(t)
        vals = rangeObject.getFloats(0)
        in_hi = max(vals[0])
        in_low = min(vals[0])
        domain = GridUtil.getSpatialDomain(rangeObject)
        
        [element_size, line_size] = domain.getLengths()
        
        for i in range(line_size):
            for j in range(element_size)[:-1]:
                vals[0][i * element_size + j] = int(abs(vals[0][i * element_size + j] - vals[0][i * element_size + j + 1]))
                
            # set last value to zero
            vals[0][i * element_size + j + 1] = 0
            
        for i in range(line_size):
            for j in range(element_size):
                vals[0][i * element_size + j] = scaleOutsideVal(vals[0][i * element_size + j], britlo, brithi)
                
        filt_low = int(min(vals[0]))
        filt_hi = int(max(vals[0]))
        if stretch == 'Contrast':
            lookup = contrast(filt_low, filt_hi, britlo, brithi, filt_low, filt_hi)
        elif stretch == 'Histogram':
            h = makeHistogram(vals, element_size, line_size, filt_low, brithi - britlo)
            lookup = histoStretch(filt_low, filt_hi, in_low, in_hi, h)
            
        vals = modify(vals, element_size, line_size, filt_low, lookup)
        rangeObject.setSamples(vals)
        
    return newData

def substitute(data, low, high, newValue):
    """change values in data  between low/high to newvalue """
    newData = data.clone();
    if (GridUtil.isTimeSequence(newData)):
        for t in range(newData.getDomainSet().getLength()):
            rangeObject = newData.getSample(t)
            values = rangeObject.getFloats(0);
            for i in range(len(values)):
               for j in range(len(values[0])):
		           if values[i][j]>=low:
		             if values[i][j]<=high:  values[i][j] = newValue;
            rangeObject.setSamples(values,1);
    else:
        rangeObject = newData;
        values = rangeObject.getFloats(0);
        for i in range(len(values)):
          for j in range(len(values[0])):
             if values[i][j]>=low:
                if values[i][j]<=high: values[0][i] = newValue;
        rangeObject.setSamples(values,1);
    return newData;

def coreFilter(sdataset1, sdataset2, user_brkpoint1='Default', user_brkpoint2='Default', user_replace1='Default', user_replace2='Default', user_stretchval='Contrast', user_britlo=0, user_brithi=255):
    """Core filter from McIDAS-X.
    
    Requires 2 source datasets; resulting image has only 2 values.
    
    Args:
        user_brkpoint1: Sdataset1 breakpoint value.
                        (default=minimum value in either source dataset)
        user_brkpoint2: Sdataset2 breakpoint value.
                        (default=maximum value in either source dataset)
        user_replace1: Success condition replacement value.
                       (default=maximum value in either source dataset)
        user_replace2: Failure condition replacement value.
                       (default=minimum vlaue in either source dataset)
        user_britlo: Minimum brightness value for the calibration.
        user_brithi: Maximum brightness value for the calibration.
    """
    data1 = sdataset1.clone()
    data2 = sdataset2.clone()
    brkpoint1 = user_brkpoint1
    brkpoint2 = user_brkpoint2
    replace1 = user_replace1
    replace2 = user_replace2
    britlo = int(user_britlo)
    brithi = int(user_brithi)
    
    if brkpoint1 != 'Default':
        brkpoint1=int(brkpoint1)
        
    if brkpoint2 != 'Default':
        brkpoint2=int(brkpoint2)
        
    if replace1 != 'Default':
        replace1=int(replace1)
        
    if replace2 != 'Default':
        replace2 = int(replace2)
        
    stretch = user_stretchval
    
    for t in range(data1.getDomainSet().getLength()):
        range1 = data1.getSample(t)
        range2 = data2.getSample(t)
        vals1 = range1.getFloats(0)
        max1 = max(vals1[0])
        min1 = min(vals1[0])
        vals2 = range2.getFloats(0)
        max2 = max(vals2[0])
        min2 = min(vals2[0])
        in_low = min([min1, min2])
        in_hi = max([max1, max2])
        
        if brkpoint1 == 'Default':
            brkpoint1 = in_low
            
        if brkpoint2 == 'Default':
            brkpoint2 = in_hi
            
        if replace1 == 'Default':
            replace1 = brkpoint2
            
        if replace2 == 'Default':
            replace2 = brkpoint1
            
        domain = GridUtil.getSpatialDomain(range1)
        [element_size, line_size] = domain.getLengths()
        for i in range(line_size):
            for j in range(element_size):
                if (vals1[0][i * element_size + j] > brkpoint1 and vals2[0][i * element_size + j] > brkpoint2):
                    vals1[0][i * element_size + j] = replace1
                else:
                    vals1[0][i * element_size + j]=replace2
                    
        for i in range(line_size):
            for j in range(element_size):
                vals1[0][i * element_size + j] = scaleOutsideVal(vals1[0][i * element_size + j], britlo, brithi)
                
        filt_low = int(min([min(vals1[0]), min2]))
        filt_hi = int(max([max(vals1[0]), max2]))
        if stretch == 'Contrast':
            lookup = contrast(filt_low, filt_hi, in_low, in_hi, filt_low, filt_hi)
        elif stretch == 'Histogram':
            # make a histogram from both datasets
            v = []
            v.append(vals1[0])
            v.append(vals2[0])
            h = makeHistogram(v, element_size, line_size, filt_low, brithi - britlo)
            lookup = histoStretch(filt_low, filt_hi, in_low, in_hi, h)
            
        vals1 = modify(vals1, element_size, line_size, filt_low, lookup)
        range1.setSamples(vals1)
        
    return data1

def exportGridToNetcdf(grid,filename='grid.nc'):
    """Write out the grid data to CF compliant netCDF file"""
    GridUtil.exportGridToNetcdf(grid,filename)
    return grid