def make_same_size(day_list,tseries_list):
    """
    Make all time-series the same size as longest time-series and
    combine them into one [ncases,ntime,...] array for plotting.
    """
    # find beg and end of period that encompasses all time-series
    day_beg = min([d.min() for d in day_list])
    day_end = max([d.max() for d in day_list])

    # find indices in 1000 year time array
    thousand_years = N.resize(dpm,12*1000).cumsum()
    ibeg = thousand_years.searchsorted(day_beg)
    iend = thousand_years.searchsorted(day_end)

    # set new array dimensions and create array
    day    = thousand_years[ibeg:iend+1]
    ncases = len(tseries_list)
    ntime  = iend -  ibeg + 1
    if tseries_list[0].ndim == 2: # regional time-series [ntime,nreg]
        data = MA.masked_all((ncases,ntime,tseries_list[0].shape[-1]),float)
    else:
        data = MA.masked_all((ncases,ntime),float)

    # fill new array according to time index
    for n in range(ncases):
        i1 = day.searchsorted(day_list[n].min())
        i2 = day.searchsorted(day_list[n].max()) + 1
        data[n,i1:i2,...] = tseries_list[n]

    return day, data

def ba_ratio_histograms(ba_files, ind_files, indices_names,minmax) :
    """computes histogram of ratio of MODIS BA to 0.5x0.5 deg occurrence
    
    Considering each day an independent measurement of the entire study area, 
    the ratio of total MODIS BA counts to total 0.5x0.5 deg cells is computed 
    for each parameter bin. Each parameter bin gets at most one observation per 
    day, and this observation embodies all 0.5x0.5deg cells in that bin for that
    day.
    """ 
    num_years = len(ind_files)
    max_days = 365
    histo_shape = zip(*minmax)[2]
    ratio_shape = histo_shape + (max_days,num_years)
    ratios = ma.masked_all(ratio_shape)
    halfdeg_counts = ma.masked_all(ratio_shape)
    
    
    ca = gca.GeoCompressedAxes(ind_files[0], 'land') 
    ca.set_clip_box(42.5, 66.5, 22, 130)
    
    for i_year in range(len(ind_files)) : 
        indfile = ind_files[i_year]
        bafile  = ba_files[i_year]
        count   = bafile.variables['count']
        timelim = len(indfile.dimensions['days'])-1
        filevars = [ indfile.variables[iname] for iname in indices_names ] 
        for i_day in range(10,timelim) : 
            print i_day
            day_data = [ f[i_day,:] for f in filevars ]
            i_conditions = zip(*day_data)
            ba_day = count[...,i_day]
            ba_total = np.sum(ba_day, axis=2)
            ba_total_cmp = ca.compress(ba_total)
            
            
            # per bin ba totals (units of modis pixels)
            burned_total = ah.AccumulatingHistogramdd(minmax=minmax)
            for i_tot,ba_tot in enumerate(ba_total_cmp) : 
                if ba_tot is ma.masked  :
                    continue
                if ba_tot > 0 : 
                    burned_total.put_record(i_conditions[i_tot], weight=ba_tot)
            
            # per bin occurrence totals (units of 0.5 deg cells)
            occurrence = ah.AccumulatingHistogramdd(minmax=minmax)
            for i_window,mask in enumerate(ca.get_vec_mask()) : 
                if not mask : 
                    occurrence.put_record(i_conditions[i_window])
            
            # calculate ratio
            i_occurrence = np.where(occurrence.H > 0)
            num_occurrence = len(i_occurrence[0])
            i_occ_oneday = i_occurrence + ( np.array([i_day]*num_occurrence), np.array([i_year]*num_occurrence))
            ratios[i_occ_oneday] = burned_total.H[i_occurrence]/occurrence.H[i_occurrence]
            halfdeg_counts[...,i_day,i_year] = occurrence.H

    ratio_histogram = compute_ratio_histo(ratios, minmax)
            
    return (ratios, halfdeg_counts, ratio_histogram)

def woa_profile_from_dap(var, d, lat, lon, depth, cfg):
    """
    Monthly Climatologic Mean and Standard Deviation from WOA,
    used either for temperature or salinity.

    INPUTS
        time: [day of the year]
        lat: [-90<lat<90]
        lon: [-180<lon<180]
        depth: [meters]

    Reads the WOA Monthly Climatology NetCDF file and
    returns the corresponding WOA values of salinity or temperature mean and
    standard deviation for the given time, lat, lon, depth.
    """
    if lon < 0:
        lon = lon+360

    url = cfg['url']

    doy = int(d.strftime('%j'))
    dataset = open_url(url)

    dn = (np.abs(doy-dataset['time'][:])).argmin()
    xn = (np.abs(lon-dataset['lon'][:])).argmin()
    yn = (np.abs(lat-dataset['lat'][:])).argmin()

    if re.match("temperature\d?$", var):
        mn = ma.masked_values(dataset.t_mn.t_mn[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.t_mn.attributes['_FillValue'])
        sd = ma.masked_values(dataset.t_sd.t_sd[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.t_sd.attributes['_FillValue'])
        # se = ma.masked_values(dataset.t_se.t_se[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.t_se.attributes['_FillValue'])
        # Use this in the future. A minimum # of samples
        # dd = ma.masked_values(dataset.t_dd.t_dd[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.t_dd.attributes['_FillValue'])
    elif re.match("salinity\d?$", var):
        mn = ma.masked_values(dataset.s_mn.s_mn[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.s_mn.attributes['_FillValue'])
        sd = ma.masked_values(dataset.s_sd.s_sd[dn, :, yn, xn].reshape(
            dataset['depth'].shape[0]), dataset.s_sd.attributes['_FillValue'])
        # dd = ma.masked_values(dataset.s_dd.s_dd[dn, :, yn, xn].reshape(
        #    dataset['depth'].shape[0]), dataset.s_dd.attributes['_FillValue'])
    zwoa = ma.array(dataset.depth[:])

    ind = (depth <= zwoa.max()) & (depth >= zwoa.min())
    # Mean value profile
    f = interp1d(zwoa[~ma.getmaskarray(mn)].compressed(), mn.compressed())
    mn_interp = ma.masked_all(depth.shape)
    mn_interp[ind] = f(depth[ind])
    # The stdev profile
    f = interp1d(zwoa[~ma.getmaskarray(sd)].compressed(), sd.compressed())
    sd_interp = ma.masked_all(depth.shape)
    sd_interp[ind] = f(depth[ind])

    output = {'woa_an': mn_interp, 'woa_sd': sd_interp}

    return output

    def __init__(self, inputdir, inputpattern=".*\.cnv",
            cfg=None, saveauxiliary=False, timeout=60):
        """
        """
        self.name = "ProfilesQCCollection"

        self.inputfiles = make_file_list(inputdir, inputpattern)

        self.profiles = process_profiles(self.inputfiles, cfg, saveauxiliary,
                timeout=timeout)
        # self.profiles = process_profiles_serial(self.inputfiles, cfg,
        #        saveauxiliary)

        self.data = {'id': [], 'profileid': [], 'profilename': []}
        self.flags = {}
        if saveauxiliary is True:
            self.auxiliary = {}

        offset = 0
        for p in self.profiles:
            N = p['timeS'].size

            # Be sure that all have the same lenght.
            for v in p.keys():
                assert p[v].size == N
            ids = offset + np.arange(N)
            self.data['id'] = np.append(self.data['id'],
                    ids).astype('i')
            profileid = [p.attributes['md5']] * N
            self.data['profileid'] = np.append(self.data['profileid'],
                    profileid)
            profilename = [p.attributes['filename']] * N
            self.data['profilename'] = np.append(self.data['profilename'],
                    profilename)
            for v in p.keys():
                if v not in self.data:
                    self.data[v] = ma.masked_all(offset)
                self.data[v] = ma.append(self.data[v], p[v])

            # ---- Dealing with the flags --------------------------------
            for v in p.flags.keys():
                if v not in self.flags:
                    self.flags[v] = {'id': [], 'profileid': []}
                self.flags[v]['id'] = np.append(self.flags[v]['id'],
                        ids).astype('i')
                self.flags[v]['profileid'] = np.append(
                        self.flags[v]['profileid'], profileid)
                for t in p.flags[v]:
                    if t not in self.flags[v]:
                        self.flags[v][t] = ma.masked_all(offset)
                    self.flags[v][t] = ma.append(self.flags[v][t],
                            p.flags[v][t])
            offset += N

        return

    def draw(self, size=1, spaces=None):
        """
        Draw samples from the transdimensional distribution.
        """
        if spaces is not None:
            if len(spaces) != size:
                raise ValueError('Sample size inconsistent with number of spaces saved')
            space_inds = np.empty(size)
            for space_id, space in enumerate(self.spaces):
                subspace = np.all(spaces == space, axis=1)
                space_inds[subspace] = space_id

        else:
            # Draws spaces randomly with the assigned weights
            cumulative_weights = np.cumsum(np.exp(self._logweights))
            space_inds = np.searchsorted(cumulative_weights, np.random.rand(size))

        draws = ma.masked_all((size, self._max_ndim))
        for space_id in range(len(self.spaces)):
            sel = space_inds == space_id
            n_fixedd = np.count_nonzero(sel)
            if n_fixedd > 0:
                # Populate only the valid entries for this parameter space
                draws[np.ix_(sel, self._spaces[space_id])] = self.kdes[space_id].draw(n_fixedd)

        return draws

def load_data():
    # Read the log file
    module_dir = os.path.dirname(__file__)
    filepath = os.path.join(module_dir, 'data', 'apm.log')
    lines = open(filepath).read().splitlines()
    
    # Parse the data
    data = dict(MAG=[], IMU=[], ATT=[])
    for line in lines:
        msgid, *fields = re.split(',\s*', line)
        if msgid in data:
            data[msgid].append([float(f) for f in fields])
    data = {key: np.asarray(val) for key, val in data.items()}
    imu = data['IMU']
    mag = data['MAG']
    
    # Build the output array
    t = np.sort(np.hstack((imu[:, 0], mag[:, 0])))
    imu_inds = np.array([tk in imu[:, 0] for tk in t])
    mag_inds = np.array([tk in mag[:, 0] for tk in t])
    y = ma.masked_all((t.size, GeneratedDTModel.ny))
    y[imu_inds, :6] = imu[:, [4, 5, 6, 1, 2, 3]]
    y[mag_inds, 6:] = mag[:, [1, 2, 3]]
    t *= 1e-3
    
    # Select the experiment interval
    range_ = np.s_[905:1800]#np.s_[900:1800]
    t = t[range_]
    y = y[range_]
    assert np.unique(t).size == t.size
    return t, y, data

    def __init__(self, nwalkers, ndim, lnpostfn, transd=False,
                 processes=None, pool=None):
        self.nwalkers = nwalkers
        self.dim = ndim
        self._kde = None
        self._kde_size = self.nwalkers
        self.updates = np.array([])

        self._get_lnpost = lnpostfn

        self.iterations = 0
        self.stored_iterations = 0

        self.pool = pool
        self.processes = processes
        if self.processes != 1 and self.pool is None:
            self.pool = Pool(self.processes)

        self._transd = transd
        if self._transd:
            self._chain = ma.masked_all((0, self.nwalkers, self.dim))
        else:
            self._chain = np.zeros((0, self.nwalkers, self.dim))

        self._lnpost = np.empty((0, self.nwalkers))
        self._lnprop = np.empty((0, self.nwalkers))
        self._acceptance = np.zeros((0, self.nwalkers))
        self._blobs = []

        self._last_run_mcmc_result = None
        self._failed_p = None

 def set_window(self, data):
     """returns an array the size of the dataset, with data correctly 
     located in the window, and all other pixels masked.
     Data must be a 2D array exactly the size of the window."""
     ds = ma.masked_all(self._dataset_shape, dtype=data.dtype)
     ds[self._window] = data
     return ds

def spike(x):
    y = ma.masked_all(x.shape, dtype=x.dtype)
    y[1:-1] = np.abs(x[1:-1] - (x[:-2] + x[2:])/2.0) - \
                np.abs((x[2:] - x[:-2])/2.0)
    # ATENTION, temporary solution
    #y[0]=0; y[-1]=0
    return y

 def test_convert_to_annual(self):
     "Test convert_to_annual"
     base = dict(D=1, H=24, T=24 * 60, S=24 * 3600)
     #for fq in ('D', 'H', 'T', 'S'):
     # Don't test for minuTe and Second frequency, too time consuming.
     for fq in ('D', 'H'):
         dates = date_array(start_date=Date(fq, '2001-01-01 00:00:00'),
                            end_date=Date(fq, '2004-12-31 23:59:59'))
         bq = base[fq]
         series = time_series(range(365 * bq) * 3 + range(366 * bq),
                              dates=dates)
         control = ma.masked_all((4, 366 * bq), dtype=series.dtype)
         control[0, :58 * bq] = range(58 * bq)
         control[0, 59 * bq:] = range(58 * bq, 365 * bq)
         control[[1, 2]] = control[0]
         control[3] = range(366 * bq)
         test = convert_to_annual(series)
         assert_equal(test, control)
     #
     series = time_series(range(59, 365) + range(366) + range(365),
                          start_date=Date('D', '2003-03-01'))
     test = convert_to_annual(series)
     assert_equal(test[:, 59:62],
                  ma.masked_values([[-1, 59, 60], [59, 60, 61], [-1, 59, 60]],
                                   - 1))

def observations(nobs, ndim):
    a = ma.masked_all((n,n), dtype=float)
    for i in range(nobs):
        (i,j) = np.random.randint(0, n, 2)
        a[(i,j)] = np.random.random_sample()
        a[(j,i)] = a[(i,j)]
    return a

def calc_uv(xvel, yvel, lons):
# Script to convert vectors from xy to uv
    u_z = ma.masked_all((xvel.shape[0],xvel.shape[1]))
    v_m = ma.masked_all((xvel.shape[0],xvel.shape[1]))
    mask = ma.getmask(xvel)
    #index = np.where(mask==False)

    for i in xrange(xvel.shape[0]):
        for j in xrange(xvel.shape[1]):
            #TO TRANSPOSE OR NOT?..
            alpha = (lons[i, j])*pi/180.
            if (mask[i, j]==False):
                u_z[i, j] = yvel[i, j]*sin(alpha) + xvel[i, j]*cos(alpha)
                v_m[i, j] = yvel[ i, j]*cos(alpha) - xvel[ i, j]*sin(alpha) 

    return u_z, v_m 

def wmean_bandpass_1D_serial(data, lshorterpass, llongerpass, t=None,
        method='hann', axis=0):
    """ Equivalent to wmean_1D_serial, but it is a bandpass

        Input:
            - data: np.array or ma.maked_array, nD
            - lshorterpass: The size of the highpass filter, i.e. shorter
                wavelenghts are preserved. It is in the same unit of t.
            - llongerpass: The size of the lowpass filter, i.e.longer
                wavelenghts are preserved. It is in the same unit of t.
	    - t: is the scale of the choosed axis, 1D. If not
                defined, it will be considered a sequence.
            - method: ['hann', 'hamming', 'blackman']
                Defines the weight function type
            - axis: Dimension which the filter will be applied
    """
    assert False, "There is a BUG here"

    assert axis <= data.ndim, "Invalid axis!"

    # If necessary, move the axis to be filtered for the first axis
    if axis != 0:
        data_smooth = wmean_bandpass_1D_serial(data.swapaxes(0, axis),
                lshorterpass = lshorterpass,
                llongerpass = llongerpass,
                t = t,
                method = method,
                axis = 0)

        return data_smooth.swapaxes(0, axis)
    # Below here, the filter will be always applied on axis=0

    # If t is not given, creates a regularly spaced t
    if t is None:
        print "The scale along the choosed axis weren't defined. I'll consider a constant sequence."
	t = np.arange(data.shape[axis])

    assert t.shape == (data.shape[axis],), "Invalid size of t."

    # ----
    winfunc = window_func(method)

    data_smooth = ma.masked_all(data.shape)

    if data.ndim==1:
        (I,) = np.nonzero(~ma.getmaskarray(data))
        for i in I:
            # First remove the high frequency
            tmp = _convolve_1D(t[i], t, llongerpass, winfunc, data)
            # Then remove the low frequency
            data_smooth[i] = tmp - \
                    _convolve_1D(t[i], t, lshorterpass, winfunc, tmp)

    else:
        I = data.shape[1]
        for i in range(I):
            data_smooth[:,i] = wmean_bandpass_1D_serial(data[:,i],
                    lshorterpass, llongerpass, t, method, axis)

    return data_smooth

def get_halfpower_period(data, filtered, dt):
    """ Returns the gain per frequency
    """
    nt, ni, nj = data.shape
    gain = ma.masked_all((nt, ni, nj))
    for i in range(ni):
        for j in range(nj):
	    if ~filtered[:,i,j].mask.all():
	        gain[:,i,j] = np.absolute(np.fft.fft(filtered[:,i,j]-filtered[:,i,j].mean())) / np.absolute(np.fft.fft(data[:,i,j]-data[:,i,j].mean()))
    gain_median = ma.masked_all(nt)
    gain_25 = ma.masked_all(nt)
    gain_75 = ma.masked_all(nt)
    # Run for each frequency, which are in the same number of timesteps
    from scipy.stats import scoreatpercentile
    for t in range(nt):
        #gain_median[t] = numpy.median(gain[t,:,:].compressed()[numpy.isfinite(gain[t,:,:].compressed())])
        #tmp = gain[t,:,:].compressed()[numpy.isfinite(gain[t,:,:].compressed())]
        #gain_median[t] = scoreatpercentile(tmp,50)
        gain_median[t] = ma.median(gain[t])
        #gain_25[t] = scoreatpercentile(tmp,25)
        #gain_75[t] = scoreatpercentile(tmp,75)

    freq = np.fft.fftfreq(nt)/dt.days

    halfpower_period = 1./freq[np.absolute(gain_median-0.5).argmin()]


    #from scipy.interpolate import UnivariateSpline
    #s = UnivariateSpline(gain_median[numpy.ceil(nt/2.):], -freq[numpy.ceil(nt/2.):], s=1)
    #xs = -freq[numpy.ceil(nt/2.):]
    #ys = s(xs)

    #import rpy2.robjects as robjects
    #smooth = robjects.r['smooth.spline'](robjects.FloatVector(gain_median[numpy.ceil(nt/2.):]),robjects.FloatVector(-freq[numpy.ceil(nt/2.):]),spar=.4)
    ##smooth = robjects.r['smooth.spline'](robjects.FloatVector(-freq[numpy.ceil(nt/2.):]),robjects.FloatVector(gain_median[numpy.ceil(nt/2.):]),spar=.4)
    #s_interp = robjects.r['predict'](smooth,x=0.5)
    ##halfpower_period = 1./s_interp.rx2['y'][0]
    #halfpower_period = 1./s_interp.rx2(2)[0]

    #smooth = robjects.r['smooth.spline'](robjects.FloatVector(-freq[numpy.ceil(nt/2.):]),robjects.FloatVector(gain_median[numpy.ceil(nt/2.):]),spar=.4)
    #s_interp = robjects.r['predict'](smooth, x = robjects.FloatVector(-freq[numpy.ceil(nt/2.):]))

    #print "Filter half window size: %s" % l
    #print "Half Power Period: %s" % halfpower_period
    #self.halfpower_period = halfpower_period

    return halfpower_period

def wmean_1D(data, l, t=None, method='hann', axis=0, interp = False):
    """ A moving window mean filter, not necessarily a regular grid.

        It is equivalent to wmean_1D_serial but run in parallel with
          multiprocesing for higher efficiency.

        Check wmean_1D_serial documentation for the inputs and other
          details.
    """
    #assert type(data)) in [np.ndarray, ma.MaskedArray]
    assert axis <= data.ndim, "Invalid axis!"

    # If necessary, move the axis to be filtered for the first axis
    if axis != 0:
        data_smooth = wmean_1D(data.swapaxes(0, axis),
            l = l,
            t = t,
            method = method,
            axis = 0,
            interp = interp)

        return data_smooth.swapaxes(0, axis)
    # Below here, the filter will be always applied on axis=0

    # If t is not given, creates a regularly spaced t
    if t is None:
        print "The scale along the choosed axis weren't defined. I'll consider a constant sequence."
	t = np.arange(data.shape[axis])

    assert t.shape == (data.shape[axis],), "Invalid size of t."

    # ----
    # Only one dimensions usually means overhead to run in parallel.
    if data.ndim==1:
        data_smooth = wmean_1D_serial(data, l, t=t, method=method, axis=axis,
                interp=interp)
        return data_smooth
    # ----

    npes = 2 * mp.cpu_count()
    pool = mp.Pool(npes)
    results = []
    I = data.shape[1]
    for i in range(I):
        results.append(pool.apply_async(wmean_1D_serial, \
                (data[:,i], l, t, method, 0, interp)))
    pool.close()

    # Collecting the results.
    if type(data) is np.ndarray:
        data_smooth = np.empty(data.shape)
    else:
        data_smooth = ma.masked_all(data.shape)

    for i, r in enumerate(results):
        data_smooth[:,i] = r.get()
    pool.terminate()

    return data_smooth

def window_mean(y,x=None,x_out=None,method="rectangular",boxsize=None):
    """Windowed means along 1-D array

    Input:
        - x [0,1,2,3,...] =>
        - x_out [x] =>
        - method [rectangular]:
            + rectangular => All data in window have same weight
        - boxsize [mean space] =>
    Output:

    Apply windowed means on a 1-D data array. Selecting adequate x_out
    and boxsize could define boxmeans or smooth filters. Method defines
    the weight method.

    An important point of this function is the ability to work with 
    unhomogenious spaced samples. Data ([1,2,3]) colected at [1,2,4] 
    times would be different if was collected at [1,2,3].
    """
    if(x==None):
        x=N.arange(N.size(y))

    if(x_out==None):
        x_out=x

    #y_out = N.zeros(N.size(x_out),N.float)
    y_out = ma.masked_all(x_out.shape)

    if(boxsize==None):
        # !!! Improve it! A better way than *1. ?!
        boxsize =(max(x)-min(x))/(N.size(x_out)*1.)

    half_boxsize = boxsize/2.

    #for x_i in x_out:
    for i in range(N.size(x_out)):
        x_i = x_out[i]


        # Higher window limit
        hi_limit = x_i+half_boxsize
        # Lower window limit
        lo_limit = x_i-half_boxsize
        # index of values inside window
        index = N.less_equal(x,hi_limit)*N.greater_equal(x,lo_limit)

        # !!! INSERT some type of check for minimum number of samples to be considered

        # x values on window around x_i
        x_tmp = N.compress(index,x)-x_i
        # y values on window
        y_tmp = N.compress(index,y)

        # weights in window according to x position
        weight = window_weight(x_tmp,boxsize,method)

        y_out[i] = N.sum(y_tmp*weight)

    return y_out

    def interpolate(self, lat, lon, var):
        """ Interpolate each var on the coordinates requested

        """

        subset, dims = self.subset(lat, lon, var)

        if np.all([y in dims['lat'] for y in lat]) & \
                np.all([x in dims['lon'] for x in lon]):
                    yn = np.nonzero([y in lat for y in dims['lat']])[0]
                    xn = np.nonzero([x in lon for x in dims['lon']])[0]
                    output = {}
                    for v in subset:
                        # output[v] = subset[v][dn, zn, yn, xn]
                        # Seriously that this is the way to do it?!!??
                        output[v] = subset[v][:, xn][yn]
                    return output

        # The output coordinates shall be created only once.
        points_out = []
        for latn in lat:
            for lonn in lon:
                points_out.append([latn, lonn])
        points_out = np.array(points_out)

        output = {}
        for v in var:
            output[v] = ma.masked_all(
                    (lat.size, lon.size),
                    dtype=subset[v].dtype)

            # The valid data
            idx = np.nonzero(~ma.getmaskarray(subset[v]))

            if idx[0].size > 0:
                points = np.array([
                    dims['lat'][idx[0]], dims['lon'][idx[1]]]).T
                values = subset[v][idx]

                # Interpolate along the dimensions that have more than one
                #   position, otherwise it means that the output is exactly
                #   on that coordinate.
                ind = np.array(
                        [np.unique(points[:, i]).size > 1 for i in
                            range(points.shape[1])])
                assert ind.any()

                values_out = griddata(
                        np.atleast_1d(np.squeeze(points[:, ind])),
                        values,
                        np.atleast_1d(np.squeeze(points_out[:, ind]))
                        )

                # Remap the interpolated value back into a 4D array
                idx = np.isfinite(values_out)
                for [y, x], out in zip(points_out[idx], values_out[idx]):
                    output[v][y==lat, x==lon] = out

        return output

 def extend_interp(datafield):
   # add masked values at southernmost end
   southernlimitmask = ma.masked_all(len(self.olon))
   olat_ext          = np.append(-82.1,self.olat)
   dfield_ext = ma.concatenate([ma.column_stack(southernlimitmask), datafield], 0)
   # f = interp2d(self.olon, olat_ext, dfield_ext)
   # return f(self.pismlon, self.pismlat)
   return interp(dfield_ext, self.olon, olat_ext, self.pismlon, self.pismlat)

def window_1Dmean_grid(data, l, method='hann', axis=0, parallel=False):
    """ A moving window mean filter applied to a regular grid.

        1D means that the filter is applied to along only one
          of the dimensions, but in the whole array. For example
          in a 3D array, each latXlon point is filtered along the
          time.

        The other types of filter consider the scales of each
          dimension. On this case it's considered a regular
          grid, so the filter can be based on the number of
          elements, and so be much optimized.

        l is in number of cells around the point being evaluated.
    """
    assert axis <= data.ndim, "Invalid axis!"

    if axis != 0:
        data_smooth = window_1Dmean_grid(data.swapaxes(0, axis), 
                l = l, 
                method = method, 
                axis = 0, 
                parallel = parallel)

        return data_smooth.swapaxes(0, axis)

    winfunc = window_func(method)
    r = np.arange(-np.floor(l/2),np.floor(l/2)+1)
    w = winfunc(r, l)

    data_smooth = ma.masked_all(data.shape)

    I = data.shape[0]
    norm = np.convolve(np.ones(I), w ,mode='same')
    if len(data.shape)==1:
	norm=numpy.convolve(numpy.ones(I),w,mode='same')
	data_smooth[:] = numpy.convolve(data[:],w,mode='same')/norm

    elif len(data.shape) == 2:
        I, J = data.shape
        for j in range(J):
            data_smooth[:,j] = np.convolve(data[:,j], w, mode='same')/norm

    elif len(data.shape) >2:
        I, J = data.shape[:2]
        for j in range(J):
                data_smooth[:,j] = window_1Dmean_grid(data[:,j], 
                        l = l, 
                        method = method, 
                        axis=0, 
                        parallel = parallel)

    try:
        data_smooth.mask = data.mask
    except:
        pass

    return data_smooth

def _weight_triangular(r,l):
    """
    """
    w = ma.masked_all(r.shape)
    ind = np.abs(r)<l/2.
    ind2 = np.abs(r)>l/2.
    w[ind] = 1-np.abs(2*r[ind]/l)
    w[ind2] = 0
    return w