def calc_area_weighted_spatial_average(dataset, area_weight=False):
    '''Calculate area weighted average of the values in OCW dataset

    :param dataset: Dataset object 
    :type dataset: :class:`dataset.Dataset`

    :returns: time series for the dataset of shape (nT)
    '''

    if dataset.lats.ndim ==1:
        lons, lats = np.meshgrid(dataset.lons, dataset.lats)
    else:
        lons = dataset.lons
        lats = dataset.lats
    weights = np.cos(lats*np.pi/180.) 

    nt, ny, nx = dataset.values.shape
    spatial_average = ma.zeros(nt)
    for it in np.arange(nt):
        if area_weight:
            spatial_average[it] = ma.average(dataset.values[it,:], weights = weights)
        else:
            spatial_average[it] = ma.average(dataset.values[it,:])

    return spatial_average

def mean_wind(prof, pbot=850, ptop=250, dp=-1, stu=0, stv=0):
    '''
    Calculates a pressure-weighted mean wind through a layer. The default
    layer is 850 to 200 hPa.

    Parameters
    ----------
    prof: profile object
        Profile object
    pbot : number (optional; default 850 hPa)
        Pressure of the bottom level (hPa)
    ptop : number (optional; default 250 hPa)
        Pressure of the top level (hPa)
    dp : negative integer (optional; default -1)
        The pressure increment for the interpolated sounding
    stu : number (optional; default 0)
        U-component of storm-motion vector
    stv : number (optional; default 0)
        V-component of storm-motion vector

    Returns
    -------
    mnu : number
        U-component
    mnv : number
        V-component

    '''
    if dp > 0: dp = -dp
    ps = np.arange(pbot, ptop+dp, dp)
    u, v = interp.components(prof, ps)
    # u -= stu; v -= stv
    return ma.average(u, weights=ps)-stu, ma.average(v, weights=ps)-stv

    def testAttrs(self):
        # Same export as testAnaNetwork, but check that the
        # attributes are synchronised
        query = dballe.Record()
        query["var"] = "B10004"
        query["datetime"] = datetime.datetime(2007, 1, 1, 0, 0, 0)
        vars = read(self.db.query_data(query), (AnaIndex(), NetworkIndex()), attributes=True)
        self.assertEqual(len(vars), 1)
        self.assertCountEqual(vars.keys(), ["B10004"])
        data = vars["B10004"]
        self.assertEqual(len(data.attrs), 2)
        self.assertCountEqual(sorted(data.attrs.keys()), ['B33007', 'B33040'])

        for net, a in ('synop', 'B33007'), ('temp', 'B33040'):
            self.assertEqual(data.dims, data.attrs[a].dims)
            self.assertEqual(data.vals.size, data.attrs[a].vals.size)
            self.assertEqual(data.vals.shape, data.attrs[a].vals.shape)

            # Find what is the network dimension where we have the attributes
            netidx = -1
            for idx, n in enumerate(data.dims[1]):
                if n == net:
                    netidx = idx
                    break
            self.assertNotEqual(netidx, -1)

            # No attrs in the other network
            self.assertEqual([x for x in data.attrs[a].vals.mask[:,1-netidx].flat], [True]*len(data.attrs[a].vals.mask[:,1-netidx].flat))
            # Same attrs as values in this network
            self.assertEqual([x for x in data.vals.mask[:,netidx].flat], [x for x in data.attrs[a].vals.mask[:,netidx].flat])
        self.assertEqual(round(ma.average(data.attrs['B33007'].vals)), 32)
        self.assertEqual(round(ma.average(data.attrs['B33040'].vals)), 54)

    def avgprices(self, stockweighted=False):
        """Return a masked array of the average price by element"""
        p = ma.array(self.prices, mask=self.prices <= 0)

        if stockweighted:
            s = ma.array(self.stock, mask=self.stock <= 0)
            avgprices = ma.average(p, weights=s, axis=1)
        else:
            #avgprices = p.sum(axis=1)/(p > 0).sum(axis=1) #denominator sums the non-zero values
            avgprices = ma.average(p, axis=1)
        return avgprices

def servTime (acqSeqL, relSeqL, dim=0):
    servTimeList = map (partial(avgWaitTime, dim=dim), acqSeqL, relSeqL)
    servTimes, servTimesSq, counts = zip (*servTimeList)
    servTimesMtx = listOfArrToMtx (servTimes)
    servTimesSqMtx = listOfArrToMtx (servTimesSq)
    countMtx = listOfArrToMtx (counts)
    # norm of columns
    norms = normalizeRowWise(countMtx.T).T
    ma_servTimesMtx = ma.array(servTimesMtx, mask = servTimesMtx == 0)
    ma_servTimesSqMtx = ma.array(servTimesSqMtx, mask = servTimesSqMtx == 0)
    return ma.average (ma_servTimesMtx, axis=0, weights=norms), ma.average(ma_servTimesSqMtx, axis=0, weights=norms)

    def testAnaTrangeNetwork(self):
        # 3 dimensions: ana, timerange, network
        # 2 variables
        query = dballe.Record(datetime=datetime.datetime(2007, 1, 1, 0, 0, 0))
        vars = read(self.db.query_data(query), (AnaIndex(), TimeRangeIndex(shared=False), NetworkIndex()))
        self.assertEqual(len(vars), 2)
        self.assertEqual(sorted(vars.keys()), ["B10004", "B13011"])

        data = vars["B10004"]
        self.assertEqual(data.name, "B10004")
        self.assertEqual(len(data.attrs), 0)
        self.assertEqual(len(data.dims), 3)
        self.assertEqual(len(data.dims[0]), 6)
        self.assertEqual(len(data.dims[1]), 1)
        self.assertEqual(len(data.dims[2]), 2)
        self.assertEqual(data.vals.size, 12)
        self.assertEqual(data.vals.shape, (6, 1, 2))
        self.assertEqual(sum(data.vals.mask.flat), 1)
        self.assertEqual(round(ma.average(data.vals)), 83185)
        self.assertEqual(data.dims[0][0], (1, 10., 15., None))
        self.assertEqual(data.dims[0][1], (2, 10., 25., None))
        self.assertEqual(data.dims[0][2], (3, 20., 15., None))
        self.assertEqual(data.dims[0][3], (4, 20., 25., None))
        self.assertEqual(data.dims[0][4], (5, 30., 15., None))
        self.assertEqual(data.dims[0][5], (6, 30., 25., None))
        self.assertEqual(data.dims[1][0], (0, None, None))
        self.assertEqual(set(data.dims[2]), set(("temp", "synop")))

        data = vars["B13011"]
        self.assertEqual(data.name, "B13011")
        self.assertEqual(len(data.attrs), 0)
        self.assertEqual(len(data.dims), 3)
        self.assertEqual(len(data.dims[0]), 6)
        self.assertEqual(len(data.dims[1]), 2)
        self.assertEqual(len(data.dims[2]), 2)
        self.assertEqual(data.vals.size, 24)
        self.assertEqual(data.vals.shape, (6, 2, 2))
        self.assertEqual(sum(data.vals.mask.flat), 0)
        self.assertAlmostEqual(ma.average(data.vals), 5.325, 6)
        self.assertEqual(data.dims[0][0], (1, 10., 15., None))
        self.assertEqual(data.dims[0][1], (2, 10., 25., None))
        self.assertEqual(data.dims[0][2], (3, 20., 15., None))
        self.assertEqual(data.dims[0][3], (4, 20., 25., None))
        self.assertEqual(data.dims[0][4], (5, 30., 15., None))
        self.assertEqual(data.dims[0][5], (6, 30., 25., None))
        self.assertEqual(data.dims[1][0], (4, -43200, 0))
        self.assertEqual(data.dims[1][1], (4, -21600, 0))
        self.assertEqual(set(data.dims[2]), set(("temp", "synop")))

        self.assertEqual(vars["B10004"].dims[0], vars["B13011"].dims[0])
        self.assertNotEqual(vars["B10004"].dims[1], vars["B13011"].dims[1])
        self.assertEqual(vars["B10004"].dims[2], vars["B13011"].dims[2])

    def average_combine(self):
        """Average combine together a set of arrays.   A CCDData object is
           returned with the data property set to the average of the arrays.
           If the data was masked or any data have been rejected, those pixels
           will not be included in the median.   A mask will be returned, and
           if a pixel has been rejected in all images, it will be masked.   The
           uncertainty of the combined image is set by the standard deviation
           of the input images.

           Returns
           -------
           combined_image: CCDData object
               CCDData object based on the combined input of CCDData objects.

        """
        #set up the data
        data, wei = ma.average(self.data_arr, axis=0, weights=self.weights,
                               returned=True)

        #set up the mask
        mask = self.data_arr.mask.sum(axis=0)
        mask = (mask == len(self.data_arr))

        #set up the variance
        uncertainty = ma.std(self.data_arr, axis=0)

        #create the combined image
        combined_image = CCDData(data.data, mask=mask, unit=self.unit,
                                 uncertainty=StdDevUncertainty(uncertainty))

        #update the meta data
        combined_image.meta['NCOMBINE'] = len(self.data_arr)

        #return the combined image
        return combined_image

def binner(x, y, w_sta, nbins, rang = None, ebar = False, per = None) :
	from numpy import array, digitize, lexsort, linspace
	from numpy.ma import average, median

	ind    = lexsort((y, x))
	xs, ys = x[ind], y[ind]

	if rang is None : mn, mx = min(xs), max(xs)
	else            : mn, mx = rang
	
	bins  = linspace(mn, mx, nbins + 1)
	x_cen = (bins[: - 1] + bins[1:])*0.5
	bins  = linspace(mn, mx, nbins)
	ibins = digitize(xs, bins)

	if w_sta   == "median" : y_sta = array([median(ys[ibins == i]) for i in range(1, bins.size + 1)])
	elif w_sta == "mean"   : y_sta = array([average(ys[ibins == i]) for i in range(1, bins.size + 1)])
	elif w_sta == "mode"   : y_sta = array([mode(ys[ibins == i])[0] for i in range(1, bins.size + 1)])

	if ebar   == False                : return x_cen, y_sta
	elif ebar == True and per == None :
		myer = abs(array([scoreatpercentile(ys[ibins == i], 15.8) for i in range(1, bins.size + 1)]) - y_sta)
		pyer = abs(array([scoreatpercentile(ys[ibins == i], 84.0) for i in range(1, bins.size + 1)]) - y_sta)
		yer  = array([myer, pyer])
		return x_cen, y_sta, yer

	elif ebar == True and per != None :
		myer = abs(array([scoreatpercentile(ys[ibins == i], per[0]) for i in range(1, bins.size + 1)]) - y_sta)
		pyer = abs(array([scoreatpercentile(ys[ibins == i], per[1]) for i in range(1, bins.size + 1)]) - y_sta)
		yer = array([myer, pyer])
		return x_cen, y_sta, yer

def temporal_rebin_with_time_index(target_dataset, nt_average):
    """ Rebin a Dataset to a new temporal resolution

    :param target_dataset: Dataset object that needs temporal rebinned
    :type target_dataset: :class:`dataset.Dataset`

    :param nt_average: Time resolution for the output datasets. 
     It is the same as the number of time indicies to be averaged. (length of time dimension in the rebinned dataset) = (original time dimension length/nt_average)
    :type temporal_resolution: integer

    :returns: A new temporally rebinned Dataset
    :rtype: :class:`dataset.Dataset`
    """
    nt = target_dataset.times.size
    if nt % nt_average !=0:
        print 'Warning: length of time dimension must be a multiple of nt_average'
    # nt2 is the length of time dimension in the rebinned dataset
    nt2 = nt/nt_average
    binned_dates = target_dataset.times[np.arange(nt2)*nt_average]
    binned_values = ma.zeros(np.insert(target_dataset.values.shape[1:],0,nt2))
    for it in np.arange(nt2):
        binned_values[it,:] = ma.average(target_dataset.values[nt_average*it:nt_average*it+nt_average,:], axis=0)
    new_dataset = ds.Dataset(target_dataset.lats,
                             target_dataset.lons,
                             binned_dates,
                             binned_values,
                             variable=target_dataset.variable,
                             units=target_dataset.units,
                             name=target_dataset.name,
                             origin=target_dataset.origin)
    return new_dataset

    def testSomeAttrs(self):
        # Same export as testAnaNetwork, but check that the
        # attributes are synchronised
        query = dballe.Record()
        query["var"] = "B10004"
        query["datetime"] = datetime.datetime(2007, 1, 1, 0, 0, 0)
        vars = read(self.db.query_data(query), (AnaIndex(), NetworkIndex()), attributes=('B33040',))
        self.assertEqual(len(vars), 1)
        self.assertCountEqual(vars.keys(), ["B10004"])
        data = vars["B10004"]
        self.assertEqual(len(data.attrs), 1)
        self.assertCountEqual(data.attrs.keys(), ['B33040'])

        a = data.attrs['B33040']
        self.assertEqual(data.dims, a.dims)
        self.assertEqual(data.vals.size, a.vals.size)
        self.assertEqual(data.vals.shape, a.vals.shape)

        # Find the temp index
        netidx = -1
        for idx, n in enumerate(data.dims[1]):
            if n == "temp":
                netidx = idx
                break
        self.assertNotEqual(netidx, -1)

        # Only compare the values on the temp index
        self.assertEqual([x for x in a.vals.mask[:,1-netidx].flat], [True]*len(a.vals.mask[:,1-netidx].flat))
        self.assertEqual([x for x in data.vals.mask[:,netidx].flat], [x for x in a.vals.mask[:,netidx].flat])
        self.assertEqual(round(ma.average(a.vals)), 54)

 def testAnaNetwork(self):
     # Ana in one dimension, network in the other
     query = dballe.Record()
     query["var"] = "B10004"
     query["datetime"] = datetime.datetime(2007, 1, 1, 0, 0, 0)
     vars = read(self.db.query_data(query), (AnaIndex(), NetworkIndex()))
     self.assertEqual(len(vars), 1)
     self.assertCountEqual(vars.keys(), ["B10004"])
     data = vars["B10004"]
     self.assertEqual(data.name, "B10004")
     self.assertEqual(len(data.attrs), 0)
     self.assertEqual(len(data.dims), 2)
     self.assertEqual(len(data.dims[0]), 6)
     self.assertEqual(len(data.dims[1]), 2)
     self.assertEqual(data.vals.size, 12)
     self.assertEqual(data.vals.shape, (6, 2))
     self.assertEqual(sum(data.vals.mask.flat), 1)
     self.assertEqual(round(ma.average(data.vals)), 83185)
     self.assertEqual(data.dims[0][0], (1, 10., 15., None))
     self.assertEqual(data.dims[0][1], (2, 10., 25., None))
     self.assertEqual(data.dims[0][2], (3, 20., 15., None))
     self.assertEqual(data.dims[0][3], (4, 20., 25., None))
     self.assertEqual(data.dims[0][4], (5, 30., 15., None))
     self.assertEqual(data.dims[0][5], (6, 30., 25., None))
     self.assertEqual(set(data.dims[1]), set(("temp", "synop")))

    def executeOperations(self, task, inputs):
        available_inputIds = [ inputId.split('-')[0] for inputId in inputs ]
        data_inputIds = [ inputId.split('-')[0] for inputId in task.inputs ]
        wids = [ inputId + "_WEIGHTS_" for inputId in data_inputIds ]
        weight_inputIds = [ ( wid if (wid in available_inputIds) else None ) for wid in wids ]
        inputs_with_weights = zip( data_inputIds, weight_inputIds )
        self.logger.info("@@@@ data_inputIds = " + str(data_inputIds) + ", weight_inputIds = " + str(weight_inputIds) + ", inputs = " + str(inputs) )
        results = []
        for input_pair in inputs_with_weights:
            input = inputs.get( input_pair[0] )  # npArray
            if( input == None ): raise Exception( "Can't find input " + input_pair[0] + " in numpyModule.WeightedAverageKernel")
            else :
                weights = inputs.get( input_pair[1] ).array if( input_pair[1] != None ) else None
                axes = self.getOrderedAxes(task,input)
                self.logger.info("\n Executing average, input: " + str( input_pair[0] ) + ", shape = " + str(input.array.shape)+ ", task metadata = " + str(task.metadata) + " Input metadata: " + str( input.metadata ) )
                t0 = time.time()
                result = input.array
                for axis in axes:
                    current_shape = list( result.shape )
                    self.logger.info(" --> Exec: axis: " + str(axis) + ", shape: " + str(current_shape) )
#                    ( result, weights ) = ma.average( result, axis, weights, True )
                    ( result, weights ) = ma.average( result, axis, np.broadcast_to( weights, current_shape ), True )
                    current_shape[axis] = 1
                    result = result.reshape( current_shape )
                    weights = weights.reshape( current_shape )

                results.append( npArray.createResult( task, input, result.filled( input.array.fill_value ) ) )
                t1 = time.time()
                self.logger.info( " ------------------------------- AVEW KERNEL: Operating on input '{0}', shape = {1}, origin = {2}, time = {3}".format( input.name, input.shape, input.origin, t1-t0 ))
        return results

def v_average_test():
    import numpy.ma as ma
    M=[[1,1,0],[0,0,1],[0,1,0]]
    Coins=[100000,200000,300000]
    Mat=numpy.matrix(M)
    Mean = ma.average(Mat, axis=0, weights=numpy.hstack(Coins))
    print(Mean)
    print(v_average(M, ReWeight(Coins)))

def mean(phi, lower=0., upper=120.):
    """wrapping the numpy.ma.average function for weighed average of masked arrays
    here weight = the image intensity,
    and the coordinates  X, Y = np.meshgrid(range(921), range(881))
        are to be averaged out.
    """
    phi1 = phi.view(ma.MaskedArray).copy()    # to be safe (slow?)
    try:
        phi1.mask += (phi1<lower) + (phi1>upper)  # masking the out-of-range regions
    except:
        phi1.mask  = (phi1<lower) + (phi1>upper)
    height, width = phi1.shape
    X, Y = np.meshgrid(range(width), range(height))
    I, J = Y, X             # always work with I, J internally
    Ibar = ma.average(I, weights=phi1)
    Jbar = ma.average(J, weights=phi1)
    return {'i':Ibar, 'j':Jbar}

    def applyOperation( self, input_variable, operation ):
        result = None
        try:
            self.setTimeBounds( input_variable )
            operator = None
#            pydevd.settrace('localhost', port=8030, stdoutToServer=False, stderrToServer=True)
            wpsLog.debug( " $$$ ApplyOperation: %s " % str( operation ) )
            if operation is not None:
                type = operation.get('type','').lower()
                bounds = operation.get('bounds','').lower()
                op_start_time = time.clock() # time.time()
                if not bounds:
                    if type == 'departures':
                        ave = cdutil.averager( input_variable, axis='t', weights='equal' )
                        result = input_variable - ave
                    elif type == 'climatology':
                        result = cdutil.averager( input_variable, axis='t', weights='equal' )
                    else:
                        result = input_variable
                    time_axis = input_variable.getTime()
                elif bounds == 'np':
                    if   type == 'departures':
                        result = ma.anomalies( input_variable ).squeeze()
                    elif type == 'climatology':
                        result = ma.average( input_variable ).squeeze()
                    else:
                        result = input_variable
                    time_axis = input_variable.getTime()
                else:
                    if bounds == 'djf': operator = cdutil.DJF
                    elif bounds == 'mam': operator = cdutil.MAM
                    elif bounds == 'jja': operator = cdutil.JJA
                    elif bounds == 'son': operator = cdutil.SON
                    elif bounds == 'year':          operator = cdutil.YEAR
                    elif bounds == 'annualcycle':   operator = cdutil.ANNUALCYCLE
                    elif bounds == 'seasonalcycle': operator = cdutil.SEASONALCYCLE
                    if operator <> None:
                        if   type == 'departures':    result = operator.departures( input_variable ).squeeze()
                        elif type == 'climatology':   result = operator.climatology( input_variable ).squeeze()
                        else:                         result = operator( input_variable ).squeeze()
                    time_axis = result.getTime()
                op_end_time = time.clock() # time.time()
                wpsLog.debug( " ---> Base Operation Time: %.5f" % (op_end_time-op_start_time) )
            else:
                result = input_variable
                time_axis = input_variable.getTime()

            if isinstance( result, float ):
                result_data = [ result ]
            elif result is not None:
                if result.__class__.__name__ == 'TransientVariable':
                    result = ma.masked_equal( result.squeeze().getValue(), input_variable.getMissing() )
                result_data = result.tolist( numpy.nan )
            else: result_data = None
        except Exception, err:
            wpsLog.debug( "Exception applying Operation '%s':\n %s" % ( str(operation), traceback.format_exc() ) )
            return ( None, None )

 def seasonal_cycle( self, input_variable ):
     t0 = time.time()
     time_vals = input_variable.getTime().asComponentTime()
     season_index_array = np.array( [  self.season_def_array[tv.month] for tv in time_vals ] )
     squeezed_input = input_variable.squeeze()
     acycle = [ ma.average( get_subset( squeezed_input, season_index, season_index_array ) ) for season_index in range(0,4) ]
     t1 = time.time()
     wpsLog.debug( "Computed seasonal cycle, time = %.4f, result:\n %s" % ( (t1-t0), str(acycle) ) )
     return ma.array(acycle)

 def annual_cycle( self, input_variable ):
     t0 = time.time()
     time_vals = input_variable.getTime().asComponentTime()
     month_index_array = np.array( [  tv.month for tv in time_vals ] )
     squeezed_input = input_variable.squeeze()
     acycle = [ ma.average( get_subset( squeezed_input, month_index, month_index_array ) ) for month_index in range(1,13) ]
     t1 = time.time()
     wpsLog.debug( "Computed annual cycle, time = %.4f, result:\n %s" % ( (t1-t0), str(acycle) ) )
     return ma.array(acycle)

    def weighted_mean_(self):
        """
        Calculates the weighted mean of the image given the probabilistic
        segmentation. If binary, mean and weighted mean will give the same
        result

        :return:
        """
        masked_seg = np.tile(self.masked_seg, [self.img_channels, 1]).T
        return ma.average(self.masked_img, axis=0, weights=masked_seg).flatten()

    def weighted_average(self,axis=0,expaxis=None):
        """ Calculate weighted average of data along axis
            after optionally inserting a new dimension into the
            shape array at position expaxis
        """

        if expaxis is not None:
            vals = ma.expand_dims(self.vals,expaxis)
            dmin = ma.expand_dims(self.dmin,expaxis)
            dmax = ma.expand_dims(self.dmax,expaxis)
            wt = ma.expand_dims(self.wt,expaxis)
        else:
            vals = self.vals
            wt = self.wt
            dmin = self.dmin
            dmax = self.dmax
        
        # Get average value
        avg,norm = ma.average(vals,axis=axis,weights=wt,returned=True)
        avg_ex = ma.expand_dims(avg,0)

        # Calculate weighted uncertainty
        wtmax = ma.max(wt,axis=axis)
        neff = norm/wtmax       # Effective number of samples based on uncertainties

        # Seeking max deviation from the average; if above avg use max, if below use min
        term = np.empty_like(vals)
        
        indices = np.where(vals > avg_ex)
        i0 = indices[0]
        irest = indices[1:]
        ii = tuple(x for x in itertools.chain([i0],irest))
        jj = tuple(x for x in itertools.chain([np.zeros_like(i0)],irest))
        term[ii] = (dmax[ii] - avg_ex[jj])**2
        
        indices = np.where(vals <= avg_ex)
        i0 = indices[0]
        irest = indices[1:]
        ii = tuple(x for x in itertools.chain([i0],irest))
        jj = tuple(x for x in itertools.chain([np.zeros_like(i0)],irest))
        term[ii] = (avg_ex[jj] - dmin[ii])**2
        
        dsum = ma.sum(term*wt,axis=0)     # Sum for weighted average of deviations

        dev = 0.5*np.sqrt(dsum/(norm*neff))
        
        if isinstance(avg,(float,np.float)):
            avg = avg_ex

        tmp_min = avg - dev
        ii = np.where(tmp_min < 0)
        tmp_min[ii] = TOL*avg[ii]
        
        return UncertContainer(avg,tmp_min,avg+dev)

def calculate_mean(data, lats):
    """
    data - a 2d lat-lon array with latitude axis first
    lats - a 1d array containing the corresponding latitude values
    
    returns - a latitude-weighted mean of the entire data array
    """
    
    # Create a 2d-array containing the weights of each cell - i.e. the cosine of the latitude of that cell
    lat_weights = np.repeat(np.cos([lats * np.pi / 180.0]).T, np.shape(data)[1], axis=1)
    return ma.average(data, weights=lat_weights)