def plot_field(self, projection='lambert', contour=True, geopolygons=None, showfig=True, vmin=None, vmax=None):
     """Plot data with contour
     """
     m=self._get_basemap(projection=projection, geopolygons=geopolygons)
     x, y=m(self.lonArr, self.latArr)
     try:
         evx, evy=m(self.evlo, self.evla)
         m.plot(evx, evy, 'yo', markersize=10)
     except:
         pass
     
     try:
         stx, sty=m(self.lonArrIn, self.latArrIn)
         m.plot(stx, sty, 'y^', markersize=10)
     except:
         pass
     im=m.pcolormesh(x, y, self.Zarr, cmap='gist_ncar_r', shading='gouraud', vmin=vmin, vmax=vmax)
     cb = m.colorbar(im, "bottom", size="3%", pad='2%')
     cb.ax.tick_params(labelsize=10)
     if self.fieldtype=='Tph' or self.fieldtype=='Tgr':
         cb.set_label('sec', fontsize=12, rotation=0)
     if self.fieldtype=='Amp':
         cb.set_label('nm', fontsize=12, rotation=0)
     if contour:
         # levels=np.linspace(ma.getdata(self.Zarr).min(), ma.getdata(self.Zarr).max(), 20)
         levels=np.linspace(ma.getdata(self.Zarr).min(), ma.getdata(self.Zarr).max(), 60)
         m.contour(x, y, self.Zarr, colors='k', levels=levels, linewidths=0.5)
     if showfig:
         plt.show()
     return

def clipdata(y, xm, x0, x1, minFlag=False):
    mx = ma.getdata(mask(xm, xm, x0, x1))
    my = ma.getdata(mask(y, xm, x0, x1))
    if minFlag:  # allow clipping from first minimum after the start time
        u = ma.argmin(mask(y, xm, x0, x1))
        mx = mx[u:-1]
        my = my[u:-1]
    return (mx, my)

def create_ac9_hs6_plot(netcdf_file_path):
    """create a png from a ac9/hs6 netcdf file"""
    plot_output_filepath = os.path.splitext(netcdf_file_path)[0] + '.png'
    dataset              = Dataset(netcdf_file_path)
    n_wavelength         = dataset.variables['wavelength'].shape[0]

    # look for main variable
    for varname in dataset.variables.keys():
        dim = dataset.variables[varname].dimensions
        if 'obs' in dim and 'wavelength' in dim:
            main_data = dataset.variables[varname]

    fig = figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')
    labels = []
    # only one profile per file in AC9 HS6 files
    for i_wl in range(n_wavelength):
        main_var_val   = main_data[:, i_wl]  # for i_prof
        depth_val      = dataset.variables['DEPTH'][:]
        wavelength_val = dataset.variables['wavelength'][i_wl]
        plot(main_var_val, depth_val)#, c=np.random.rand(3, 1))
        labels.append('%s nm' % wavelength_val)

    station_name = ''.join(ma.getdata(dataset.variables['station_name'][dataset.variables['station_index'][0] - 1]))
    title('%s\nCruise: %s\n Station %s' % (dataset.source, dataset.cruise_id, station_name))
    gca().invert_yaxis()
    xlabel('%s: %s in %s' % (main_data.name, main_data.long_name, main_data.units))
    ylabel('%s in %s; positive %s' % (dataset.variables['DEPTH'].long_name,
                                      dataset.variables['DEPTH'].units,
                                      dataset.variables['DEPTH'].positive))
    legend(labels, loc='upper right', prop=fontP, title='Wavelength')
    savefig(plot_output_filepath)
    plt.close(fig)

    def make_gene_map_1(self):
        count = 0
        self.iterator = itertools.product(range(self.time_len), range(self.lat_len), range(self.lon_len))
        for x_valid in self.iterator:
            binary_string = bin(count)[2:]
            while len(binary_string) < self.string_length:
                binary_string = "0" + binary_string
                # self.gene_map[binary_string] = {}
                # self.gene_map[binary_string]['coordinate'] = tuple(x_valid)
            if ma.getdata(self.array[x_valid]) < 100:  # chooses unmasked points?
                self.gene_map[binary_string] = {}
                self.gene_map[binary_string]["coordinate"] = tuple(x_valid)
                self.gene_map[binary_string]["value"] = self.array[x_valid]
                # count += 1
                print count, binary_string, self.gene_map[binary_string]["value"]  # debug
            else:
                # self.gene_map[binary_string]['value'] = 1E06
                pass
            count += 1
            # print count, binary_string, self.gene_map[binary_string]['value'] # debug
        self.last_valid_binary_string = binary_string
        not_valid_first = eval("0b" + binary_string) + 1
        not_valid_last = eval("0b" + "1" * self.string_length)
        print not_valid_last
        print x_valid
        print ma.isMA(self.array)

        for x_not_valid in range(not_valid_first, not_valid_last + 1):
            binary_string = bin(x_not_valid)[2:]
            self.gene_map[binary_string] = {}
            self.gene_map[binary_string]["coordinate"] = (999, 999, 999)
            self.gene_map[binary_string]["value"] = 1e06
            print x_not_valid, binary_string, self.gene_map[binary_string]["value"]

def apply_on_fields(series, func, *args, **kwargs):
    """
    Apply the function ``func`` to each field of ``series``.

    Parameters
    ----------
    series : array-like
        A (subclass) of ndarray.
    func : function
        Function to apply.
    args : var
        Additional arguments to ``func``.
    kwargs : var
        Additional optional parameters to ``func``
    """
    names = series.dtype.names
    if names is None:
        return func(series, *args, **kwargs)
    results = [func(series[f], *args, **kwargs) for f in names]
    rdtype = [(f[0], r.dtype) for (f, r) in zip(series.dtype.descr, results)]
    rdata = [ma.getdata(r) for r in results]
    rmask = [ma.getmaskarray(r) for r in results]
    isarray = isinstance(results[0], np.ndarray)
    if isarray:
        output = ma.array(zip(*rdata), mask=zip(*rmask),
                          dtype=rdtype).view(type=type(series))
        output._update_from(series)
    else:
        output = ma.array(tuple(rdata), mask=tuple(rmask), dtype=rdtype)
    return output

def putdata_raw(imgname, data, clone=None):
    """Store (overwrite) data in an existing CASA image.
       See getdata_raw(imgname) for the reverse operation.

       Parameters
       ----------
       imagename : str
           The (absolute) CASA image filename.  It should exist
           already, unless **clone** was given.

       data : 2D numpy array or a list of 2D numpy arrays
           The data...

       clone : str, optional
           An optional filename from which to clone the image
           for output. It needs to be an absolute filename.
  
    """
    if clone != None:
        taskinit.ia.fromimage(infile=clone,outfile=imgname,overwrite=True) 
        taskinit.ia.close()
    # @todo this seems circumvent to have to borrow the odd dimensions (nx,ny,1,1,1) shape was seen
    if type(data) == type([]):
        # @todo since this needs to extend the axes, the single plane clone and replace data doesn't work here
        raise Exception,"Not Implemented Yet"
        bigim = ia.imageconcat(outfile=imgname, infiles=infiles, axis=2, relax=T, tempclose=F, overwrite=T)
        bigim.close()
    else:
        taskinit.tb.open(imgname,nomodify=False)
        d = taskinit.tb.getcol('map')
        pdata = ma.getdata(data).reshape(d.shape)
        taskinit.tb.putcol('map',pdata)
        taskinit.tb.flush()
        taskinit.tb.close()
    return

 def ma2np(self):
     """Convert all the maksed data array to numpy array
     """
     self.Zarr=ma.getdata(self.Zarr)
     try:
         self.diffaArr=ma.getdata(self.diffaArr)
     except:
         pass
     try:
         self.appV=ma.getdata(self.appV)
     except:
         pass
     try:
         self.lplc=ma.getdata(self.lplc)
     except:
         pass
     return

    def _calculateStats(self, di, dof=0.):
        """Calculate the core QA statistics on a difference image"""
        mask = di.getMask()
        maskArr = di.getMask().getArray()

        # Create a mask using BAD, SAT, NO_DATA, EDGE bits.  Keep detections
        maskArr &= mask.getPlaneBitMask(["BAD", "SAT", "NO_DATA", "EDGE"])

        # Mask out values based on maskArr
        diArr = ma.array(di.getImage().getArray(), mask=maskArr)
        varArr = ma.array(di.getVariance().getArray(), mask=maskArr)

        # Normalize by sqrt variance, units are in sigma
        diArr /= np.sqrt(varArr)
        mean = diArr.mean()

        # This is the maximum-likelihood extimate of the variance stdev**2
        stdev = diArr.std()
        median = ma.extras.median(diArr)

        # Compute IQR of just un-masked data
        data = ma.getdata(diArr[~diArr.mask])
        iqr = np.percentile(data, 75.) - np.percentile(data, 25.)

        #Calculte chisquare of the residual
        chisq=np.sum(np.power(data, 2.))

        # Mean squared error: variance + bias**2
        # Bias = |data - model| = mean of diffim
        # Variance = |(data - model)**2| = mean of diffim**2
        bias = mean
        variance = np.power(data, 2.).mean()
        mseResids = bias**2 + variance

        # If scipy is not set up, return zero for the stats
        try:
            #In try block because of risk of divide by zero
            rchisq=chisq/(len(data)-1-dof)
            # K-S test on the diffim to a Normal distribution
            import scipy.stats
            D, prob = scipy.stats.kstest(data, 'norm')

            A2, crit, sig = scipy.stats.anderson(data, 'norm')
            # Anderson Darling statistic cand be inf for really non-Gaussian distributions.
            if np.isinf(A2) or np.isnan(A2):
                A2 = 9999.
        except ZeroDivisionError:
            D = 0.
            prob = 0.
            A2 = 0.
            crit = np.zeros(5)
            sig = np.zeros(5)
            rchisq = 0

        return {"mean": mean, "stdev": stdev, "median": median, "iqr": iqr,
                "D": D, "prob": prob, "A2": A2, "crit": crit, "sig": sig,
                "rchisq": rchisq, "mseResids": mseResids}

 def plot(self, ds=1000, unit='km', cmap='seismic_r', vmin=None, vmax=None, zsize=10):
     """Plot velocity model
     =============================================================================
     Input Parameters:
     ds              - grid spacing
     unit            - unit
     vmin, vmax      - vmin,vmax for colorbar
     =============================================================================
     """
     
     # XLength=self.xmax-self.xmin
     # ZLength=self.zmax-self.zmin
     # xsize=zsize*(XLength/ZLength)
     # fig = plt.figure(figsize=(xsize, zsize))
     if cmap=='ses3d':
         cmap = colormaps.make_colormap({0.0:[0.1,0.0,0.0], 0.2:[0.8,0.0,0.0], 0.3:[1.0,0.7,0.0],0.48:[0.92,0.92,0.92],
             0.5:[0.92,0.92,0.92], 0.52:[0.92,0.92,0.92], 0.7:[0.0,0.6,0.7], 0.8:[0.0,0.0,0.8], 1.0:[0.0,0.0,0.1]})
     if self.plotflag==False:
         raise ValueError('No plot array!')
     # plt.figure(figsize=(16,13))
     if self.regular==True:
         im=plt.pcolormesh(self.XArrPlot/ds, self.ZArrPlot/ds, self.VsArrPlot/ds, cmap=cmap, vmin=vmin, vmax=vmax)
     else:
         xi = np.linspace(self.xmin, self.xmax, self.Nx*10)
         zi = np.linspace(self.zmin, self.zmax, self.Nz*10)
         self.xi, self.zi = np.meshgrid(xi, zi)
         #-- Interpolating at the points in xi, yi
         self.vi = griddata(self.XArr, self.ZArr, self.VsArr, self.xi, self.zi, 'linear')
         im=plt.pcolormesh(self.xi/ds, self.zi/ds, ma.getdata(self.vi)/ds, cmap=cmap, vmin=vmin, vmax=vmax, shading='gouraud')
     ##########################################
     plt.plot(500., 1000 , 'r*', markersize=30, lw=3)
     # plt.plot( [0., 4000.], [1000, 1000] , 'b--', lw=3)
     # plt.plot( [500., 500.], [700., 1300.] , 'g-', lw=3)
     # plt.plot( [500., 3500.], [700., 700.] , 'g-', lw=3)
     # plt.plot( [500., 3500.], [1300., 1300.] , 'g-', lw=3)
     # plt.plot( [3500., 3500.], [700., 1300.] , 'g-', lw=3)
     # 
     # plt.plot( [0., 0.], [0., 2000.] , 'k-', lw=3)
     # plt.plot( [0., 4000.], [0., 0.] , 'k-', lw=3)
     # plt.plot( [4000., 4000.], [0., 2000.] , 'k-', lw=3)
     # plt.plot( [0., 4000.], [2000., 2000.] , 'k-', lw=3)
     ##########################################
     plt.xlabel('x('+unit+')', fontsize=35)
     plt.ylabel('z('+unit+')', fontsize=35)
     # plt.axis([self.xmin/ds, self.xmax/ds, self.zmin/ds, self.zmax/ds])
     plt.axis('scaled')
     cb=plt.colorbar(shrink=0.8)#, size="3%", pad='2%')
     cb.set_label('Vs (km/s)', fontsize=20, rotation=90)
     plt.yticks(fontsize=30)
     plt.xticks(fontsize=30)
     ########################
     # plt.ylim([-100, 2100])
     # plt.xlim([-100, 4100])
     ########################
     plt.show()
     return

    def test_set_fields(self):
        # Tests setting fields.
        base = self.base.copy()
        mbase = base.view(mrecarray)
        mbase = mbase.copy()
        mbase.fill_value = (999999, 1e20, 'N/A')
        # Change the data, the mask should be conserved
        mbase.a._data[:] = 5
        assert_equal(mbase['a']._data, [5, 5, 5, 5, 5])
        assert_equal(mbase['a']._mask, [0, 1, 0, 0, 1])
        # Change the elements, and the mask will follow
        mbase.a = 1
        assert_equal(mbase['a']._data, [1]*5)
        assert_equal(ma.getmaskarray(mbase['a']), [0]*5)
        # Use to be _mask, now it's recordmask
        assert_equal(mbase.recordmask, [False]*5)
        assert_equal(mbase._mask.tolist(),
                     np.array([(0, 0, 0),
                               (0, 1, 1),
                               (0, 0, 0),
                               (0, 0, 0),
                               (0, 1, 1)],
                              dtype=bool))
        # Set a field to mask ........................
        mbase.c = masked
        # Use to be mask, and now it's still mask !
        assert_equal(mbase.c.mask, [1]*5)
        assert_equal(mbase.c.recordmask, [1]*5)
        assert_equal(ma.getmaskarray(mbase['c']), [1]*5)
        assert_equal(ma.getdata(mbase['c']), [asbytes('N/A')]*5)
        assert_equal(mbase._mask.tolist(),
                     np.array([(0, 0, 1),
                               (0, 1, 1),
                               (0, 0, 1),
                               (0, 0, 1),
                               (0, 1, 1)],
                              dtype=bool))
        # Set fields by slices .......................
        mbase = base.view(mrecarray).copy()
        mbase.a[3:] = 5
        assert_equal(mbase.a, [1, 2, 3, 5, 5])
        assert_equal(mbase.a._mask, [0, 1, 0, 0, 0])
        mbase.b[3:] = masked
        assert_equal(mbase.b, base['b'])
        assert_equal(mbase.b._mask, [0, 1, 0, 1, 1])
        # Set fields globally..........................
        ndtype = [('alpha', '|S1'), ('num', int)]
        data = ma.array([('a', 1), ('b', 2), ('c', 3)], dtype=ndtype)
        rdata = data.view(MaskedRecords)
        val = ma.array([10, 20, 30], mask=[1, 0, 0])

        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            rdata['num'] = val
            assert_equal(rdata.num, val)
            assert_equal(rdata.num.mask, [1, 0, 0])

def forward_fill(marr, maxgap=None):
    """
    Forward fills masked values in a 1-d array when there are less ``maxgap``
    consecutive masked values.

    Parameters
    ----------
    marr : MaskedArray
        Series to fill
    maxgap : {int}, optional
        Maximum gap between consecutive masked values.
        If ``maxgap`` is not specified, all masked values are forward-filled.


    Examples
    --------
    >>> x = ma.arange(20)
    >>> x[(x%5)!=0] = ma.masked
    >>> print x
    [0 -- -- -- -- 5 -- -- -- -- 10 -- -- -- -- 15 -- -- -- --]
    >>> print forward_fill(x)
    [0 0 0 0 0 5 5 5 5 5 10 10 10 10 10 15 15 15 15 15]

    """
    # !!!: We should probably port that to C.
    # Initialization ..................
    if np.ndim(marr) > 1:
        raise ValueError,"The input array should be 1D only!"
    a = ma.array(marr, copy=True)
    amask = getmask(a)
    if amask is nomask or a.size == 0:
        return a
    #
    adata = getdata(a)
    # Get the indices of the masked values (except a[0])
    idxtofill = amask[1:].nonzero()[0] + 1
    currGap = 0
    if maxgap is not None:
        previdx = -1
        for i in idxtofill:
            if i != previdx + 1:
                currGap = 0
            currGap += 1
            if currGap <= maxgap and not amask[i-1]:
                adata[i] = adata[i-1]
                amask[i] = False
                previdx = i
            else:
                amask[i-maxgap:i] = True
    else:
        for i in idxtofill:
            if not amask[i-1]:
                adata[i] = adata[i-1]
                amask[i] = False
    return a

def fromarrays(
    arraylist,
    dtype=None,
    shape=None,
    formats=None,
    names=None,
    titles=None,
    aligned=False,
    byteorder=None,
    fill_value=None,
):
    """Creates a mrecarray from a (flat) list of masked arrays.

    Parameters
    ----------
    arraylist : sequence
        A list of (masked) arrays. Each element of the sequence is first converted
        to a masked array if needed. If a 2D array is passed as argument, it is
        processed line by line
    dtype : {None, dtype}, optional
        Data type descriptor.
    shape : {None, integer}, optional
        Number of records. If None, shape is defined from the shape of the
        first array in the list.
    formats : {None, sequence}, optional
        Sequence of formats for each individual field. If None, the formats will
        be autodetected by inspecting the fields and selecting the highest dtype
        possible.
    names : {None, sequence}, optional
        Sequence of the names of each field.
    fill_value : {None, sequence}, optional
        Sequence of data to be used as filling values.

    Notes
    -----
    Lists of tuples should be preferred over lists of lists for faster processing.
    """
    datalist = [getdata(x) for x in arraylist]
    masklist = [np.atleast_1d(getmaskarray(x)) for x in arraylist]
    _array = recfromarrays(
        datalist,
        dtype=dtype,
        shape=shape,
        formats=formats,
        names=names,
        titles=titles,
        aligned=aligned,
        byteorder=byteorder,
    ).view(mrecarray)
    _array._mask.flat = zip(*masklist)
    if fill_value is not None:
        _array.fill_value = fill_value
    return _array

 def natgridInterp(self, interp='nn', copy=True, bounds_error=False, fill_value=np.nan):
     if self.Xarr.size == self.XarrIn.size:
         if (np.abs(self.Xarr.reshape(self.Nx*self.Ny)-self.XarrIn)).sum() < 0.01\
             and (np.abs(self.Yarr.reshape(self.Nx*self.Ny)-self.YarrIn)).sum() < 0.01:
             print 'No need to interpolate!'
             self.Zarr=self.ZarrIn
             return
     # self.Zarr = ma.getdata(griddata(self.XarrIn, self.YarrIn, self.ZarrIn, self.Xarr, self.Yarr, interp=interp))
     self.Zarr = griddata(self.XarrIn, self.YarrIn, self.ZarrIn, self.Xarr, self.Yarr, interp=interp)
     self.Zarr = ma.getdata(self.Zarr)
     # self.Zarr=self.Zarr.reshape(self.Nx*self.Ny)
     return

def create_pigment_tss_plot(netcdf_file_path):
    """create a png from a pigment netcdf file"""
    plot_output_filepath = os.path.splitext(netcdf_file_path)[0] + '.png'
    dataset              = Dataset(netcdf_file_path)

    profiles          = dataset.variables['profile']
    n_obs_per_profile = dataset.variables['row_size']
    n_profiles        = len(profiles)

    if 'CPHL_a' in dataset.variables.keys():
        main_data = dataset.variables['CPHL_a']
    elif 'SPM' in dataset.variables.keys():
        main_data = dataset.variables['SPM']
    else:
        find_main_var_name = (set(dataset.variables.keys()) - set([u'TIME', u'LATITUDE', u'LONGITUDE', u'station_name', u'station_index', u'profile', u'row_size', u'DEPTH'])).pop()
        main_data = dataset.variables[find_main_var_name]

    fig = figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')
    labels = []
    ylim([-1, max(dataset.variables['DEPTH'][:]) + 1])
    for i_prof in range(n_profiles):
        # we look for the observations indexes related to the choosen profile
        idx_obs_start = sum(n_obs_per_profile[0:i_prof])
        idx_obs_end   = idx_obs_start + n_obs_per_profile[i_prof] - 1
        idx_obs       = range(idx_obs_start, idx_obs_end + 1)

        main_var_val = main_data[idx_obs]  # for i_prof
        depth_val    = dataset.variables['DEPTH'][idx_obs]

        if len(main_var_val) == 1:
            scatter(main_var_val, depth_val)
        elif len(main_var_val) > 1:
            plot(main_var_val, depth_val, '--')#, c=np.random.rand(3, 1))
        else:
            scatter(main_var_val, depth_val, c=np.random.rand(3, 1))

        station_name = ''.join(ma.getdata(dataset.variables['station_name'][dataset.variables['station_index'][i_prof] - 1]))
        labels.append(station_name)

    gca().invert_yaxis()
    title('%s\nCruise: %s' % (dataset.source, dataset.cruise_id))
    xlabel('%s: %s in %s' % (main_data.name, main_data.long_name, main_data.units))
    ylabel('%s in %s; positive %s' % (dataset.variables['DEPTH'].long_name,
                                      dataset.variables['DEPTH'].units,
                                      dataset.variables['DEPTH'].positive))
    try:
        legend(labels, loc='upper left', prop=fontP, title='Station')
    except:
        pass

    savefig(plot_output_filepath)
    plt.close(fig)

 def fftDiff2(self, m, n):
     Nx=1<<(self.Nx-1).bit_length()
     Ny=1<<(self.Ny-1).bit_length()
     # h = np.fft.fft2(self.Zarr, s=[Nx, Ny] )
     h = np.fft.fft2(ma.getdata(self.Zarr), s=[Nx, Ny] )
     hshift = np.fft.fftshift(h)
     u = np.arange(Nx) - Nx/2.
     v = np.arange(Ny) - Ny/2.
     U,V = np.meshgrid(u,v)
     hdiff = ( (1j*2*np.pi*U/Nx)**m )*( (1j*2*np.pi*V/Ny)**n ) * hshift
     out_diff = np.real( np.fft.ifft2( np.fft.ifftshift(hdiff) ) )/(self.dx**m)/(self.dy**n)
     out_diff = out_diff[:self.Ny, :self.Nx]
     return out_diff

def interp_data(var, xpts, ypts, xpts2m, ypts2m, int_meth='cubic'):
#interpoalte data onto a regular 2d grid. Used by interp_uv
    data = ma.getdata(var)
    mask = ma.getmask(var)
    index = np.where(mask==False)
    data_masked = data[index]
    xpoints = xpts[index]
    ypoints = ypts[index]
    points = (xpoints, ypoints)
    var_int = griddata(points,data_masked,(xpts2m,ypts2m),method=int_meth)
    var_int = ma.masked_array(var_int,np.isnan(var_int))

    return var_int

def Fill2ThetaAzimuthMap(masks,TA,tam,image):
    'Needs a doc string'
    Zlim = masks['Thresholds'][1]
    rings = masks['Rings']
    arcs = masks['Arcs']
    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
    for tth,thick in rings:
        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
    for tth,azm,thick in arcs:
        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
        tam = ma.mask_or(tam.flatten(),tamt*tama)
    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
    tabs = np.ones_like(taz)
    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
    return tax,tay,taz,tad,tabs

 def generate_corrected_map(self, dataid, glbdir, outdir, pers=np.array([]), glbpfx='smpkolya_phv_R_', outpfx='smpkolya_phv_R_'):
     """
     Generate corrected global phave velocity map using a regional phase velocity map.
     =================================================================================================================
     Input Parameters:
     dataid              - dataid for regional phase velocity map
     glbdir              - location of global reference phase velocity map files
     outdir              - output directory
     pers                - period array for correction (default is 4)
     glbpfx              - prefix for global reference phase velocity map files
     outpfx              - prefix for output reference phase velocity map files
     -----------------------------------------------------------------------------------------------------------------
     Output format:
     outdir/outpfx+str(int(per))
     =================================================================================================================
     """
     if not os.path.isdir(outdir):
         os.makedirs(outdir)
     if pers.size==0:
         pers=np.append( np.arange(7.)*10.+40., np.arange(2.)*25.+125.)
     for per in pers:
         inglobalfname=glbdir+'/'+glbpfx+str(int(per))
         try:
             self.get_data4plot(dataid=dataid, period=per)
         except:
             print 'No regional data for period =',per,'sec'
             continue
         if not os.path.isfile(inglobalfname):
             print 'No global data for period =',per,'sec'
             continue
         outfname=outdir+'/'+outpfx+'%g' %(per)
         InglbArr=np.loadtxt(inglobalfname)
         outArr=InglbArr.copy()
         lonArr=self.lonArr.reshape(self.lonArr.size)
         latArr=self.latArr.reshape(self.latArr.size)
         vel_iso=ma.getdata(self.vel_iso)
         vel_iso=vel_iso.reshape(vel_iso.size)
         for i in xrange(InglbArr[:,0].size):
             lonG=InglbArr[i,0]
             latG=InglbArr[i,1]
             phVG=InglbArr[i,2]
             for j in xrange(lonArr.size):
                 lonR=lonArr[j]
                 latR=latArr[j]
                 phVR=vel_iso[j]
                 if abs(lonR-lonG)<0.05 and abs(latR-latG)<0.05 and phVR!=0:
                     outArr[i,2]=phVR
         np.savetxt(outfname, outArr, fmt='%g %g %.4f')
     return

 def plot_propagation(self, projection='lambert', inbasemap=None, factor=3, showfig=False):
     """Plot propagation direction
     """
     if inbasemap==None:
         m=self._get_basemap(projection=projection)
     else:
         m=inbasemap
     if self.lonArr.shape[0]-2==self.grad[0].shape[0] and self.lonArr.shape[1]-2==self.grad[0].shape[1]:
         self.cut_edge(1,1)
     elif self.lonArr.shape[0]!=self.grad[0].shape[0] or self.lonArr.shape[1]!=self.grad[0].shape[1]:
         raise ValueError('Incompatible shape for gradient and lon/lat array!')
     normArr = np.sqrt ( ma.getdata(self.grad[0] )** 2 + ma.getdata(self.grad[1]) ** 2)
     x, y=m(self.lonArr, self.latArr)
     U=self.grad[1]/normArr
     V=self.grad[0]/normArr
     if factor!=None:
         x=x[0:self.Nlat:factor, 0:self.Nlon:factor]
         y=y[0:self.Nlat:factor, 0:self.Nlon:factor]
         U=U[0:self.Nlat:factor, 0:self.Nlon:factor]
         V=V[0:self.Nlat:factor, 0:self.Nlon:factor]
     Q = m.quiver(x, y, U, V, scale=50, width=0.001)
     if showfig:
         plt.show()
     return

    def test_set_fields(self):
        "Tests setting fields."
        base = self.base.copy()
        mbase = base.view(mrecarray)
        mbase = mbase.copy()
        mbase.fill_value = (999999, 1e20, "N/A")
        # Change the data, the mask should be conserved
        mbase.a._data[:] = 5
        assert_equal(mbase["a"]._data, [5, 5, 5, 5, 5])
        assert_equal(mbase["a"]._mask, [0, 1, 0, 0, 1])
        # Change the elements, and the mask will follow
        mbase.a = 1
        assert_equal(mbase["a"]._data, [1] * 5)
        assert_equal(ma.getmaskarray(mbase["a"]), [0] * 5)
        assert_equal(mbase._mask, [False] * 5)
        assert_equal(
            mbase._fieldmask.tolist(), np.array([(0, 0, 0), (0, 1, 1), (0, 0, 0), (0, 0, 0), (0, 1, 1)], dtype=bool)
        )
        # Set a field to mask ........................
        mbase.c = masked
        assert_equal(mbase.c.mask, [1] * 5)
        assert_equal(ma.getmaskarray(mbase["c"]), [1] * 5)
        assert_equal(ma.getdata(mbase["c"]), ["N/A"] * 5)
        assert_equal(
            mbase._fieldmask.tolist(), np.array([(0, 0, 1), (0, 1, 1), (0, 0, 1), (0, 0, 1), (0, 1, 1)], dtype=bool)
        )
        # Set fields by slices .......................
        mbase = base.view(mrecarray).copy()
        mbase.a[3:] = 5
        assert_equal(mbase.a, [1, 2, 3, 5, 5])
        assert_equal(mbase.a._mask, [0, 1, 0, 0, 0])
        mbase.b[3:] = masked
        assert_equal(mbase.b, base["b"])
        assert_equal(mbase.b._mask, [0, 1, 0, 1, 1])
        # Set fields globally..........................
        ndtype = [("alpha", "|S1"), ("num", int)]
        data = ma.array([("a", 1), ("b", 2), ("c", 3)], dtype=ndtype)
        rdata = data.view(MaskedRecords)
        val = ma.array([10, 20, 30], mask=[1, 0, 0])
        #
        import warnings

        warnings.simplefilter("ignore")
        rdata["num"] = val
        assert_equal(rdata.num, val)
        assert_equal(rdata.num.mask, [1, 0, 0])