astrobase/lcproc/catalogs.py
def filter_lclist(lc_catalog,
                  objectidcol='objectid',
                  racol='ra',
                  declcol='decl',
                  xmatchexternal=None,
                  xmatchdistarcsec=3.0,
                  externalcolnums=(0,1,2),
                  externalcolnames=['objectid','ra','decl'],
                  externalcoldtypes='U20,f8,f8',
                  externalcolsep=None,
                  externalcommentchar='#',
                  conesearch=None,
                  conesearchworkers=1,
                  columnfilters=None,
                  field_fitsfile=None,
                  field_wcsfrom=None,
                  field_scale=ZScaleInterval(),
                  field_stretch=LinearStretch(),
                  field_colormap=plt.cm.gray_r,
                  field_findersize=None,
                  field_pltopts={'marker':'o',
                                 'markersize':10.0,
                                 'markerfacecolor':'none',
                                 'markeredgewidth':2.0,
                                 'markeredgecolor':'red'},
                  field_grid=False,
                  field_gridcolor='k',
                  field_zoomcontain=True,
                  copylcsto=None):

    '''This is used to perform cone-search, cross-match, and column-filter
    operations on a light curve catalog generated by `make_lclist`.

    Uses the output of `make_lclist` above. This function returns a list of
    light curves matching various criteria specified by the `xmatchexternal`,
    `conesearch`, and `columnfilters kwargs`. Use this function to generate
    input lists for other lcproc functions,
    e.g. `lcproc.lcvfeatures.parallel_varfeatures`,
    `lcproc.periodfinding.parallel_pf`, and `lcproc.lcbin.parallel_timebin`,
    among others.

    The operations are applied in this order if more than one is specified:
    `xmatchexternal` -> `conesearch` -> `columnfilters`. All results from these
    operations are joined using a logical AND operation.

    Parameters
    ----------

    objectidcol : str
        This is the name of the object ID column in the light curve catalog.

    racol : str
        This is the name of the RA column in the light curve catalog.

    declcol : str
        This is the name of the Dec column in the light curve catalog.

    xmatchexternal : str or None
        If provided, this is the filename of a text file containing objectids,
        ras and decs to match the objects in the light curve catalog to by their
        positions.

    xmatchdistarcsec : float
        This is the distance in arcseconds to use when cross-matching to the
        external catalog in `xmatchexternal`.

    externalcolnums : sequence of int
        This a list of the zero-indexed column numbers of columns to extract
        from the external catalog file.

    externalcolnames : sequence of str
        This is a list of names of columns that will be extracted from the
        external catalog file. This is the same length as
        `externalcolnums`. These must contain the names provided as the
        `objectid`, `ra`, and `decl` column names so this function knows which
        column numbers correspond to those columns and can use them to set up
        the cross-match.

    externalcoldtypes : str
        This is a CSV string containing numpy dtype definitions for all columns
        listed to extract from the external catalog file. The number of dtype
        definitions should be equal to the number of columns to extract.

    externalcolsep : str or None
        The column separator to use when extracting columns from the external
        catalog file. If None, any whitespace between columns is used as the
        separator.

    externalcommentchar : str
        The character indicating that a line in the external catalog file is to
        be ignored.

    conesearch : list of float

        This is used to specify cone-search parameters. It should be a three
        element list:

        [center_ra_deg, center_decl_deg, search_radius_deg]

    conesearchworkers : int
        The number of parallel workers to launch for the cone-search operation.

    columnfilters : list of str

        This is a list of strings indicating any filters to apply on each column
        in the light curve catalog. All column filters are applied in the
        specified sequence and are combined with a logical AND operator. The
        format of each filter string should be:

        '<lc_catalog column>|<operator>|<operand>'

        where:

        - <lc_catalog column> is a column in the lc_catalog pickle file
        - <operator> is one of: 'lt', 'gt', 'le', 'ge', 'eq', 'ne', which
          correspond to the usual operators: <, >, <=, >=, ==, != respectively.
        - <operand> is a float, int, or string.

    field_fitsfile : str or None
        If this is not None, it should be the path to a FITS image containing
        the objects these light curves are for. If this is provided,
        `make_lclist` will use the WCS information in the FITS itself if
        `field_wcsfrom` is None (or from a WCS header file pointed to by
        `field_wcsfrom`) to obtain x and y pixel coordinates for all of the
        objects in the field. A finder chart will also be made using
        `astrobase.plotbase.fits_finder_chart` using the corresponding
        `field_scale`, `_stretch`, `_colormap`, `_findersize`, `_pltopts`,
        `_grid`, and `_gridcolors` kwargs for that function, reproduced here to
        enable customization of the finder chart plot.

    field_wcsfrom : str or None
        If `wcsfrom` is None, the WCS to transform the RA/Dec to pixel x/y will
        be taken from the FITS header of `fitsfile`. If this is not None, it
        must be a FITS or similar file that contains a WCS header in its first
        extension.

    field_scale : astropy.visualization.Interval object
        `scale` sets the normalization for the FITS pixel values. This is an
        astropy.visualization Interval object.
        See http://docs.astropy.org/en/stable/visualization/normalization.html
        for details on `scale` and `stretch` objects.

    field_stretch : astropy.visualization.Stretch object
        `stretch` sets the stretch function for mapping FITS pixel values to
        output pixel values. This is an astropy.visualization Stretch object.
        See http://docs.astropy.org/en/stable/visualization/normalization.html
        for details on `scale` and `stretch` objects.

    field_colormap : matplotlib Colormap object
        `colormap` is a matplotlib color map object to use for the output image.

    field_findersize : None or tuple of two ints
        If `findersize` is None, the output image size will be set by the NAXIS1
        and NAXIS2 keywords in the input `fitsfile` FITS header. Otherwise,
        `findersize` must be a tuple with the intended x and y size of the image
        in inches (all output images will use a DPI = 100).

    field_pltopts : dict
        `field_pltopts` controls how the overlay points will be plotted. This
        a dict with standard matplotlib marker, etc. kwargs as key-val pairs,
        e.g. 'markersize', 'markerfacecolor', etc. The default options make red
        outline circles at the location of each object in the overlay.

    field_grid : bool
        `grid` sets if a grid will be made on the output image.

    field_gridcolor : str
        `gridcolor` sets the color of the grid lines. This is a usual matplotib
        color spec string.

    field_zoomcontain : bool
        `field_zoomcontain` controls if the finder chart will be zoomed to
        just contain the overlayed points. Everything outside the footprint of
        these points will be discarded.

    copylcsto : str
        If this is provided, it is interpreted as a directory target to copy
        all the light curves that match the specified conditions.

    Returns
    -------

    tuple
        Returns a two elem tuple: (matching_object_lcfiles, matching_objectids)
        if conesearch and/or column filters are used. If `xmatchexternal` is
        also used, a three-elem tuple is returned: (matching_object_lcfiles,
        matching_objectids, extcat_matched_objectids).

    '''

    with open(lc_catalog,'rb') as infd:
        lclist = pickle.load(infd)

    # generate numpy arrays of the matching object indexes. we do it this way so
    # we can AND everything at the end, instead of having to look up the objects
    # at these indices and running the columnfilter on them
    xmatch_matching_index = np.full_like(lclist['objects'][objectidcol],
                                         False,
                                         dtype=np.bool)
    conesearch_matching_index = np.full_like(lclist['objects'][objectidcol],
                                             False,
                                             dtype=np.bool)

    # do the xmatch first
    ext_matches = []
    ext_matching_objects = []

    if (xmatchexternal and
        isinstance(xmatchexternal, str) and
        os.path.exists(xmatchexternal)):

        try:

            # read in the external file
            extcat = np.genfromtxt(xmatchexternal,
                                   usecols=externalcolnums,
                                   delimiter=externalcolsep,
                                   names=externalcolnames,
                                   dtype=externalcoldtypes,
                                   comments=externalcommentchar)

            ext_cosdecl = np.cos(np.radians(extcat['decl']))
            ext_sindecl = np.sin(np.radians(extcat['decl']))
            ext_cosra = np.cos(np.radians(extcat['ra']))
            ext_sinra = np.sin(np.radians(extcat['ra']))

            ext_xyz = np.column_stack((ext_cosra*ext_cosdecl,
                                       ext_sinra*ext_cosdecl,
                                       ext_sindecl))
            ext_xyzdist = 2.0 * np.sin(np.radians(xmatchdistarcsec/3600.0)/2.0)

            # get our kdtree
            our_kdt = lclist['kdtree']

            # get the external kdtree
            ext_kdt = sps.cKDTree(ext_xyz)

            # do a query_ball_tree
            extkd_matchinds = ext_kdt.query_ball_tree(our_kdt, ext_xyzdist)

            for extind, mind in enumerate(extkd_matchinds):
                if len(mind) > 0:
                    ext_matches.append(mind[0])

                    # get the whole matching row for the ext objects recarray
                    ext_matching_objects.append(extcat[extind])

            ext_matches = np.array(ext_matches)

            if ext_matches.size > 0:

                # update the xmatch_matching_index
                xmatch_matching_index[ext_matches] = True

                LOGINFO('xmatch: objects matched to %s within %.1f arcsec: %s' %
                        (xmatchexternal, xmatchdistarcsec, ext_matches.size))

            else:

                LOGERROR("xmatch: no objects were cross-matched to external "
                         "catalog spec: %s, can't continue" % xmatchexternal)
                return None, None, None


        except Exception as e:

            LOGEXCEPTION('could not match to external catalog spec: %s' %
                         repr(xmatchexternal))
            raise


    # do the cone search next
    if (conesearch and
        isinstance(conesearch, (list, tuple)) and
        len(conesearch) == 3):

        try:

            racenter, declcenter, searchradius = conesearch
            cosdecl = np.cos(np.radians(declcenter))
            sindecl = np.sin(np.radians(declcenter))
            cosra = np.cos(np.radians(racenter))
            sinra = np.sin(np.radians(racenter))

            # this is the search distance in xyz unit vectors
            xyzdist = 2.0 * np.sin(np.radians(searchradius)/2.0)

            # get the kdtree
            our_kdt = lclist['kdtree']

            # look up the coordinates
            kdtindices = our_kdt.query_ball_point([cosra*cosdecl,
                                                   sinra*cosdecl,
                                                   sindecl],
                                                  xyzdist,
                                                  n_jobs=conesearchworkers)

            if kdtindices and len(kdtindices) > 0:

                LOGINFO('cone search: objects within %.4f deg '
                        'of (%.3f, %.3f): %s' %
                        (searchradius, racenter, declcenter, len(kdtindices)))

                # update the conesearch_matching_index
                matchingind = kdtindices
                conesearch_matching_index[np.array(matchingind)] = True

            # we fail immediately if we found nothing. this assumes the user
            # cares more about the cone-search than the regular column filters
            else:

                LOGERROR("cone-search: no objects were found within "
                         "%.4f deg of (%.3f, %.3f): %s, can't continue" %
                         (searchradius, racenter, declcenter, len(kdtindices)))
                return None, None


        except Exception as e:

            LOGEXCEPTION('cone-search: could not run a cone-search, '
                         'is there a kdtree present in %s?' % lc_catalog)
            raise


    # now that we're done with cone-search, do the column filtering
    allfilterinds = []
    if columnfilters and isinstance(columnfilters, list):

        # go through each filter
        for cfilt in columnfilters:

            try:

                fcol, foperator, foperand = cfilt.split('|')
                foperator = FILTEROPS[foperator]

                # generate the evalstring
                filterstr = (
                    "np.isfinite(lclist['objects']['%s']) & "
                    "(lclist['objects']['%s'] %s %s)"
                ) % (fcol, fcol, foperator, foperand)
                filterind = eval(filterstr)

                ngood = lclist['objects'][objectidcol][filterind].size
                LOGINFO('filter: %s -> objects matching: %s ' % (cfilt, ngood))

                allfilterinds.append(filterind)

            except Exception as e:

                LOGEXCEPTION('filter: could not understand filter spec: %s'
                             % cfilt)
                LOGWARNING('filter: not applying this broken filter')


    # now that we have all the filter indices good to go
    # logical-AND all the things

    # make sure we only do filtering if we were told to do so
    if (xmatchexternal or conesearch or columnfilters):

        filterstack = []
        if xmatchexternal:
            filterstack.append(xmatch_matching_index)
        if conesearch:
            filterstack.append(conesearch_matching_index)
        if columnfilters:
            filterstack.extend(allfilterinds)

        finalfilterind = np.column_stack(filterstack)
        finalfilterind = np.all(finalfilterind, axis=1)

        # get the filtered object light curves and object names
        filteredobjectids = lclist['objects'][objectidcol][finalfilterind]
        filteredlcfnames = lclist['objects']['lcfname'][finalfilterind]

    else:

        filteredobjectids = lclist['objects'][objectidcol]
        filteredlcfnames = lclist['objects']['lcfname']


    # if we're told to make a finder chart with the selected objects
    if field_fitsfile is not None and os.path.exists(field_fitsfile):

        # get the RA and DEC of the matching objects
        matching_ra = lclist['objects'][racol][finalfilterind]
        matching_decl = lclist['objects'][declcol][finalfilterind]

        matching_postfix = []

        if xmatchexternal is not None:
            matching_postfix.append(
                'xmatch_%s' %
                os.path.splitext(os.path.basename(xmatchexternal))[0]
            )
        if conesearch is not None:
            matching_postfix.append('conesearch_RA%.3f_DEC%.3f_RAD%.5f' %
                                    tuple(conesearch))

        if columnfilters is not None:
            for cfi, cf in enumerate(columnfilters):
                if cfi == 0:
                    matching_postfix.append('filter_%s_%s_%s' %
                                            tuple(cf.split('|')))
                else:
                    matching_postfix.append('_and_%s_%s_%s' %
                                            tuple(cf.split('|')))

        if len(matching_postfix) > 0:
            matching_postfix = '-%s' % '_'.join(matching_postfix)
        else:
            matching_postfix = ''

        # next, we'll make a PNG plot for the finder
        finder_outfile = os.path.join(
            os.path.dirname(lc_catalog),
            '%s%s.png' %
            (os.path.splitext(os.path.basename(lc_catalog))[0],
             matching_postfix)
        )

        finder_png = fits_finder_chart(
            field_fitsfile,
            finder_outfile,
            wcsfrom=field_wcsfrom,
            scale=field_scale,
            stretch=field_stretch,
            colormap=field_colormap,
            findersize=field_findersize,
            overlay_ra=matching_ra,
            overlay_decl=matching_decl,
            overlay_pltopts=field_pltopts,
            field_zoomcontain=field_zoomcontain,
            grid=field_grid,
            gridcolor=field_gridcolor
        )

        if finder_png is not None:
            LOGINFO('generated a finder PNG '
                    'with an object position overlay '
                    'for this filtered LC list: %s' % finder_png)



    # if copylcsto is not None, copy LCs over to it
    if copylcsto is not None:

        if not os.path.exists(copylcsto):
            os.mkdir(copylcsto)

        if TQDM:
            lciter = tqdm(filteredlcfnames)
        else:
            lciter = filteredlcfnames

        LOGINFO('copying matching light curves to %s' % copylcsto)

        for lc in lciter:
            shutil.copy(lc, copylcsto)


    LOGINFO('done. objects matching all filters: %s' % filteredobjectids.size)

    if xmatchexternal and len(ext_matching_objects) > 0:
        return filteredlcfnames, filteredobjectids, ext_matching_objects
    else:
        return filteredlcfnames, filteredobjectids
def filter_lclist(lc_catalog,
                  objectidcol='objectid',
                  racol='ra',
                  declcol='decl',
                  xmatchexternal=None,
                  xmatchdistarcsec=3.0,
                  externalcolnums=(0,1,2),
                  externalcolnames=['objectid','ra','decl'],
                  externalcoldtypes='U20,f8,f8',
                  externalcolsep=None,
                  externalcommentchar='#',
                  conesearch=None,
                  conesearchworkers=1,
                  columnfilters=None,
                  field_fitsfile=None,
                  field_wcsfrom=None,
                  field_scale=ZScaleInterval(),
                  field_stretch=LinearStretch(),
                  field_colormap=plt.cm.gray_r,
                  field_findersize=None,
                  field_pltopts={'marker':'o',
                                 'markersize':10.0,
                                 'markerfacecolor':'none',
                                 'markeredgewidth':2.0,
                                 'markeredgecolor':'red'},
                  field_grid=False,
                  field_gridcolor='k',
                  field_zoomcontain=True,
                  copylcsto=None):

    '''This is used to perform cone-search, cross-match, and column-filter
    operations on a light curve catalog generated by `make_lclist`.

    Uses the output of `make_lclist` above. This function returns a list of
    light curves matching various criteria specified by the `xmatchexternal`,
    `conesearch`, and `columnfilters kwargs`. Use this function to generate
    input lists for other lcproc functions,
    e.g. `lcproc.lcvfeatures.parallel_varfeatures`,
    `lcproc.periodfinding.parallel_pf`, and `lcproc.lcbin.parallel_timebin`,
    among others.

    The operations are applied in this order if more than one is specified:
    `xmatchexternal` -> `conesearch` -> `columnfilters`. All results from these
    operations are joined using a logical AND operation.

    Parameters
    ----------

    objectidcol : str
        This is the name of the object ID column in the light curve catalog.

    racol : str
        This is the name of the RA column in the light curve catalog.

    declcol : str
        This is the name of the Dec column in the light curve catalog.

    xmatchexternal : str or None
        If provided, this is the filename of a text file containing objectids,
        ras and decs to match the objects in the light curve catalog to by their
        positions.

    xmatchdistarcsec : float
        This is the distance in arcseconds to use when cross-matching to the
        external catalog in `xmatchexternal`.

    externalcolnums : sequence of int
        This a list of the zero-indexed column numbers of columns to extract
        from the external catalog file.

    externalcolnames : sequence of str
        This is a list of names of columns that will be extracted from the
        external catalog file. This is the same length as
        `externalcolnums`. These must contain the names provided as the
        `objectid`, `ra`, and `decl` column names so this function knows which
        column numbers correspond to those columns and can use them to set up
        the cross-match.

    externalcoldtypes : str
        This is a CSV string containing numpy dtype definitions for all columns
        listed to extract from the external catalog file. The number of dtype
        definitions should be equal to the number of columns to extract.

    externalcolsep : str or None
        The column separator to use when extracting columns from the external
        catalog file. If None, any whitespace between columns is used as the
        separator.

    externalcommentchar : str
        The character indicating that a line in the external catalog file is to
        be ignored.

    conesearch : list of float

        This is used to specify cone-search parameters. It should be a three
        element list:

        [center_ra_deg, center_decl_deg, search_radius_deg]

    conesearchworkers : int
        The number of parallel workers to launch for the cone-search operation.

    columnfilters : list of str

        This is a list of strings indicating any filters to apply on each column
        in the light curve catalog. All column filters are applied in the
        specified sequence and are combined with a logical AND operator. The
        format of each filter string should be:

        '<lc_catalog column>|<operator>|<operand>'

        where:

        - <lc_catalog column> is a column in the lc_catalog pickle file
        - <operator> is one of: 'lt', 'gt', 'le', 'ge', 'eq', 'ne', which
          correspond to the usual operators: <, >, <=, >=, ==, != respectively.
        - <operand> is a float, int, or string.

    field_fitsfile : str or None
        If this is not None, it should be the path to a FITS image containing
        the objects these light curves are for. If this is provided,
        `make_lclist` will use the WCS information in the FITS itself if
        `field_wcsfrom` is None (or from a WCS header file pointed to by
        `field_wcsfrom`) to obtain x and y pixel coordinates for all of the
        objects in the field. A finder chart will also be made using
        `astrobase.plotbase.fits_finder_chart` using the corresponding
        `field_scale`, `_stretch`, `_colormap`, `_findersize`, `_pltopts`,
        `_grid`, and `_gridcolors` kwargs for that function, reproduced here to
        enable customization of the finder chart plot.

    field_wcsfrom : str or None
        If `wcsfrom` is None, the WCS to transform the RA/Dec to pixel x/y will
        be taken from the FITS header of `fitsfile`. If this is not None, it
        must be a FITS or similar file that contains a WCS header in its first
        extension.

    field_scale : astropy.visualization.Interval object
        `scale` sets the normalization for the FITS pixel values. This is an
        astropy.visualization Interval object.
        See http://docs.astropy.org/en/stable/visualization/normalization.html
        for details on `scale` and `stretch` objects.

    field_stretch : astropy.visualization.Stretch object
        `stretch` sets the stretch function for mapping FITS pixel values to
        output pixel values. This is an astropy.visualization Stretch object.
        See http://docs.astropy.org/en/stable/visualization/normalization.html
        for details on `scale` and `stretch` objects.

    field_colormap : matplotlib Colormap object
        `colormap` is a matplotlib color map object to use for the output image.

    field_findersize : None or tuple of two ints
        If `findersize` is None, the output image size will be set by the NAXIS1
        and NAXIS2 keywords in the input `fitsfile` FITS header. Otherwise,
        `findersize` must be a tuple with the intended x and y size of the image
        in inches (all output images will use a DPI = 100).

    field_pltopts : dict
        `field_pltopts` controls how the overlay points will be plotted. This
        a dict with standard matplotlib marker, etc. kwargs as key-val pairs,
        e.g. 'markersize', 'markerfacecolor', etc. The default options make red
        outline circles at the location of each object in the overlay.

    field_grid : bool
        `grid` sets if a grid will be made on the output image.

    field_gridcolor : str
        `gridcolor` sets the color of the grid lines. This is a usual matplotib
        color spec string.

    field_zoomcontain : bool
        `field_zoomcontain` controls if the finder chart will be zoomed to
        just contain the overlayed points. Everything outside the footprint of
        these points will be discarded.

    copylcsto : str
        If this is provided, it is interpreted as a directory target to copy
        all the light curves that match the specified conditions.

    Returns
    -------

    tuple
        Returns a two elem tuple: (matching_object_lcfiles, matching_objectids)
        if conesearch and/or column filters are used. If `xmatchexternal` is
        also used, a three-elem tuple is returned: (matching_object_lcfiles,
        matching_objectids, extcat_matched_objectids).

    '''

    with open(lc_catalog,'rb') as infd:
        lclist = pickle.load(infd)

    # generate numpy arrays of the matching object indexes. we do it this way so
    # we can AND everything at the end, instead of having to look up the objects
    # at these indices and running the columnfilter on them
    xmatch_matching_index = np.full_like(lclist['objects'][objectidcol],
                                         False,
                                         dtype=np.bool)
    conesearch_matching_index = np.full_like(lclist['objects'][objectidcol],
                                             False,
                                             dtype=np.bool)

    # do the xmatch first
    ext_matches = []
    ext_matching_objects = []

    if (xmatchexternal and
        isinstance(xmatchexternal, str) and
        os.path.exists(xmatchexternal)):

        try:

            # read in the external file
            extcat = np.genfromtxt(xmatchexternal,
                                   usecols=externalcolnums,
                                   delimiter=externalcolsep,
                                   names=externalcolnames,
                                   dtype=externalcoldtypes,
                                   comments=externalcommentchar)

            ext_cosdecl = np.cos(np.radians(extcat['decl']))
            ext_sindecl = np.sin(np.radians(extcat['decl']))
            ext_cosra = np.cos(np.radians(extcat['ra']))
            ext_sinra = np.sin(np.radians(extcat['ra']))

            ext_xyz = np.column_stack((ext_cosra*ext_cosdecl,
                                       ext_sinra*ext_cosdecl,
                                       ext_sindecl))
            ext_xyzdist = 2.0 * np.sin(np.radians(xmatchdistarcsec/3600.0)/2.0)

            # get our kdtree
            our_kdt = lclist['kdtree']

            # get the external kdtree
            ext_kdt = sps.cKDTree(ext_xyz)

            # do a query_ball_tree
            extkd_matchinds = ext_kdt.query_ball_tree(our_kdt, ext_xyzdist)

            for extind, mind in enumerate(extkd_matchinds):
                if len(mind) > 0:
                    ext_matches.append(mind[0])

                    # get the whole matching row for the ext objects recarray
                    ext_matching_objects.append(extcat[extind])

            ext_matches = np.array(ext_matches)

            if ext_matches.size > 0:

                # update the xmatch_matching_index
                xmatch_matching_index[ext_matches] = True

                LOGINFO('xmatch: objects matched to %s within %.1f arcsec: %s' %
                        (xmatchexternal, xmatchdistarcsec, ext_matches.size))

            else:

                LOGERROR("xmatch: no objects were cross-matched to external "
                         "catalog spec: %s, can't continue" % xmatchexternal)
                return None, None, None


        except Exception as e:

            LOGEXCEPTION('could not match to external catalog spec: %s' %
                         repr(xmatchexternal))
            raise


    # do the cone search next
    if (conesearch and
        isinstance(conesearch, (list, tuple)) and
        len(conesearch) == 3):

        try:

            racenter, declcenter, searchradius = conesearch
            cosdecl = np.cos(np.radians(declcenter))
            sindecl = np.sin(np.radians(declcenter))
            cosra = np.cos(np.radians(racenter))
            sinra = np.sin(np.radians(racenter))

            # this is the search distance in xyz unit vectors
            xyzdist = 2.0 * np.sin(np.radians(searchradius)/2.0)

            # get the kdtree
            our_kdt = lclist['kdtree']

            # look up the coordinates
            kdtindices = our_kdt.query_ball_point([cosra*cosdecl,
                                                   sinra*cosdecl,
                                                   sindecl],
                                                  xyzdist,
                                                  n_jobs=conesearchworkers)

            if kdtindices and len(kdtindices) > 0:

                LOGINFO('cone search: objects within %.4f deg '
                        'of (%.3f, %.3f): %s' %
                        (searchradius, racenter, declcenter, len(kdtindices)))

                # update the conesearch_matching_index
                matchingind = kdtindices
                conesearch_matching_index[np.array(matchingind)] = True

            # we fail immediately if we found nothing. this assumes the user
            # cares more about the cone-search than the regular column filters
            else:

                LOGERROR("cone-search: no objects were found within "
                         "%.4f deg of (%.3f, %.3f): %s, can't continue" %
                         (searchradius, racenter, declcenter, len(kdtindices)))
                return None, None


        except Exception as e:

            LOGEXCEPTION('cone-search: could not run a cone-search, '
                         'is there a kdtree present in %s?' % lc_catalog)
            raise


    # now that we're done with cone-search, do the column filtering
    allfilterinds = []
    if columnfilters and isinstance(columnfilters, list):

        # go through each filter
        for cfilt in columnfilters:

            try:

                fcol, foperator, foperand = cfilt.split('|')
                foperator = FILTEROPS[foperator]

                # generate the evalstring
                filterstr = (
                    "np.isfinite(lclist['objects']['%s']) & "
                    "(lclist['objects']['%s'] %s %s)"
                ) % (fcol, fcol, foperator, foperand)
                filterind = eval(filterstr)

                ngood = lclist['objects'][objectidcol][filterind].size
                LOGINFO('filter: %s -> objects matching: %s ' % (cfilt, ngood))

                allfilterinds.append(filterind)

            except Exception as e:

                LOGEXCEPTION('filter: could not understand filter spec: %s'
                             % cfilt)
                LOGWARNING('filter: not applying this broken filter')


    # now that we have all the filter indices good to go
    # logical-AND all the things

    # make sure we only do filtering if we were told to do so
    if (xmatchexternal or conesearch or columnfilters):

        filterstack = []
        if xmatchexternal:
            filterstack.append(xmatch_matching_index)
        if conesearch:
            filterstack.append(conesearch_matching_index)
        if columnfilters:
            filterstack.extend(allfilterinds)

        finalfilterind = np.column_stack(filterstack)
        finalfilterind = np.all(finalfilterind, axis=1)

        # get the filtered object light curves and object names
        filteredobjectids = lclist['objects'][objectidcol][finalfilterind]
        filteredlcfnames = lclist['objects']['lcfname'][finalfilterind]

    else:

        filteredobjectids = lclist['objects'][objectidcol]
        filteredlcfnames = lclist['objects']['lcfname']


    # if we're told to make a finder chart with the selected objects
    if field_fitsfile is not None and os.path.exists(field_fitsfile):

        # get the RA and DEC of the matching objects
        matching_ra = lclist['objects'][racol][finalfilterind]
        matching_decl = lclist['objects'][declcol][finalfilterind]

        matching_postfix = []

        if xmatchexternal is not None:
            matching_postfix.append(
                'xmatch_%s' %
                os.path.splitext(os.path.basename(xmatchexternal))[0]
            )
        if conesearch is not None:
            matching_postfix.append('conesearch_RA%.3f_DEC%.3f_RAD%.5f' %
                                    tuple(conesearch))

        if columnfilters is not None:
            for cfi, cf in enumerate(columnfilters):
                if cfi == 0:
                    matching_postfix.append('filter_%s_%s_%s' %
                                            tuple(cf.split('|')))
                else:
                    matching_postfix.append('_and_%s_%s_%s' %
                                            tuple(cf.split('|')))

        if len(matching_postfix) > 0:
            matching_postfix = '-%s' % '_'.join(matching_postfix)
        else:
            matching_postfix = ''

        # next, we'll make a PNG plot for the finder
        finder_outfile = os.path.join(
            os.path.dirname(lc_catalog),
            '%s%s.png' %
            (os.path.splitext(os.path.basename(lc_catalog))[0],
             matching_postfix)
        )

        finder_png = fits_finder_chart(
            field_fitsfile,
            finder_outfile,
            wcsfrom=field_wcsfrom,
            scale=field_scale,
            stretch=field_stretch,
            colormap=field_colormap,
            findersize=field_findersize,
            overlay_ra=matching_ra,
            overlay_decl=matching_decl,
            overlay_pltopts=field_pltopts,
            field_zoomcontain=field_zoomcontain,
            grid=field_grid,
            gridcolor=field_gridcolor
        )

        if finder_png is not None:
            LOGINFO('generated a finder PNG '
                    'with an object position overlay '
                    'for this filtered LC list: %s' % finder_png)



    # if copylcsto is not None, copy LCs over to it
    if copylcsto is not None:

        if not os.path.exists(copylcsto):
            os.mkdir(copylcsto)

        if TQDM:
            lciter = tqdm(filteredlcfnames)
        else:
            lciter = filteredlcfnames

        LOGINFO('copying matching light curves to %s' % copylcsto)

        for lc in lciter:
            shutil.copy(lc, copylcsto)


    LOGINFO('done. objects matching all filters: %s' % filteredobjectids.size)

    if xmatchexternal and len(ext_matching_objects) > 0:
        return filteredlcfnames, filteredobjectids, ext_matching_objects
    else:
        return filteredlcfnames, filteredobjectids
def filter_lclist ( lc_catalog , objectidcol = 'objectid' , racol = 'ra' , declcol = 'decl' , xmatchexternal = None , xmatchdistarcsec = 3.0 , externalcolnums = ( 0 , 1 , 2 ) , externalcolnames = [ 'objectid' , 'ra' , 'decl' ] , externalcoldtypes = 'U20,f8,f8' , externalcolsep = None , externalcommentchar = '#' , conesearch = None , conesearchworkers = 1 , columnfilters = None , field_fitsfile = None , field_wcsfrom = None , field_scale = ZScaleInterval ( ) , field_stretch = LinearStretch ( ) , field_colormap = plt . cm . gray_r , field_findersize = None , field_pltopts = { 'marker' : 'o' , 'markersize' : 10.0 , 'markerfacecolor' : 'none' , 'markeredgewidth' : 2.0 , 'markeredgecolor' : 'red' } , field_grid = False , field_gridcolor = 'k' , field_zoomcontain = True , copylcsto = None ) : with open ( lc_catalog , 'rb' ) as infd : lclist = pickle . load ( infd ) # generate numpy arrays of the matching object indexes. we do it this way so # we can AND everything at the end, instead of having to look up the objects # at these indices and running the columnfilter on them xmatch_matching_index = np . full_like ( lclist [ 'objects' ] [ objectidcol ] , False , dtype = np . bool ) conesearch_matching_index = np . full_like ( lclist [ 'objects' ] [ objectidcol ] , False , dtype = np . bool ) # do the xmatch first ext_matches = [ ] ext_matching_objects = [ ] if ( xmatchexternal and isinstance ( xmatchexternal , str ) and os . path . exists ( xmatchexternal ) ) : try : # read in the external file extcat = np . genfromtxt ( xmatchexternal , usecols = externalcolnums , delimiter = externalcolsep , names = externalcolnames , dtype = externalcoldtypes , comments = externalcommentchar ) ext_cosdecl = np . cos ( np . radians ( extcat [ 'decl' ] ) ) ext_sindecl = np . sin ( np . radians ( extcat [ 'decl' ] ) ) ext_cosra = np . cos ( np . radians ( extcat [ 'ra' ] ) ) ext_sinra = np . sin ( np . radians ( extcat [ 'ra' ] ) ) ext_xyz = np . column_stack ( ( ext_cosra * ext_cosdecl , ext_sinra * ext_cosdecl , ext_sindecl ) ) ext_xyzdist = 2.0 * np . sin ( np . radians ( xmatchdistarcsec / 3600.0 ) / 2.0 ) # get our kdtree our_kdt = lclist [ 'kdtree' ] # get the external kdtree ext_kdt = sps . cKDTree ( ext_xyz ) # do a query_ball_tree extkd_matchinds = ext_kdt . query_ball_tree ( our_kdt , ext_xyzdist ) for extind , mind in enumerate ( extkd_matchinds ) : if len ( mind ) > 0 : ext_matches . append ( mind [ 0 ] ) # get the whole matching row for the ext objects recarray ext_matching_objects . append ( extcat [ extind ] ) ext_matches = np . array ( ext_matches ) if ext_matches . size > 0 : # update the xmatch_matching_index xmatch_matching_index [ ext_matches ] = True LOGINFO ( 'xmatch: objects matched to %s within %.1f arcsec: %s' % ( xmatchexternal , xmatchdistarcsec , ext_matches . size ) ) else : LOGERROR ( "xmatch: no objects were cross-matched to external " "catalog spec: %s, can't continue" % xmatchexternal ) return None , None , None except Exception as e : LOGEXCEPTION ( 'could not match to external catalog spec: %s' % repr ( xmatchexternal ) ) raise # do the cone search next if ( conesearch and isinstance ( conesearch , ( list , tuple ) ) and len ( conesearch ) == 3 ) : try : racenter , declcenter , searchradius = conesearch cosdecl = np . cos ( np . radians ( declcenter ) ) sindecl = np . sin ( np . radians ( declcenter ) ) cosra = np . cos ( np . radians ( racenter ) ) sinra = np . sin ( np . radians ( racenter ) ) # this is the search distance in xyz unit vectors xyzdist = 2.0 * np . sin ( np . radians ( searchradius ) / 2.0 ) # get the kdtree our_kdt = lclist [ 'kdtree' ] # look up the coordinates kdtindices = our_kdt . query_ball_point ( [ cosra * cosdecl , sinra * cosdecl , sindecl ] , xyzdist , n_jobs = conesearchworkers ) if kdtindices and len ( kdtindices ) > 0 : LOGINFO ( 'cone search: objects within %.4f deg ' 'of (%.3f, %.3f): %s' % ( searchradius , racenter , declcenter , len ( kdtindices ) ) ) # update the conesearch_matching_index matchingind = kdtindices conesearch_matching_index [ np . array ( matchingind ) ] = True # we fail immediately if we found nothing. this assumes the user # cares more about the cone-search than the regular column filters else : LOGERROR ( "cone-search: no objects were found within " "%.4f deg of (%.3f, %.3f): %s, can't continue" % ( searchradius , racenter , declcenter , len ( kdtindices ) ) ) return None , None except Exception as e : LOGEXCEPTION ( 'cone-search: could not run a cone-search, ' 'is there a kdtree present in %s?' % lc_catalog ) raise # now that we're done with cone-search, do the column filtering allfilterinds = [ ] if columnfilters and isinstance ( columnfilters , list ) : # go through each filter for cfilt in columnfilters : try : fcol , foperator , foperand = cfilt . split ( '|' ) foperator = FILTEROPS [ foperator ] # generate the evalstring filterstr = ( "np.isfinite(lclist['objects']['%s']) & " "(lclist['objects']['%s'] %s %s)" ) % ( fcol , fcol , foperator , foperand ) filterind = eval ( filterstr ) ngood = lclist [ 'objects' ] [ objectidcol ] [ filterind ] . size LOGINFO ( 'filter: %s -> objects matching: %s ' % ( cfilt , ngood ) ) allfilterinds . append ( filterind ) except Exception as e : LOGEXCEPTION ( 'filter: could not understand filter spec: %s' % cfilt ) LOGWARNING ( 'filter: not applying this broken filter' ) # now that we have all the filter indices good to go # logical-AND all the things # make sure we only do filtering if we were told to do so if ( xmatchexternal or conesearch or columnfilters ) : filterstack = [ ] if xmatchexternal : filterstack . append ( xmatch_matching_index ) if conesearch : filterstack . append ( conesearch_matching_index ) if columnfilters : filterstack . extend ( allfilterinds ) finalfilterind = np . column_stack ( filterstack ) finalfilterind = np . all ( finalfilterind , axis = 1 ) # get the filtered object light curves and object names filteredobjectids = lclist [ 'objects' ] [ objectidcol ] [ finalfilterind ] filteredlcfnames = lclist [ 'objects' ] [ 'lcfname' ] [ finalfilterind ] else : filteredobjectids = lclist [ 'objects' ] [ objectidcol ] filteredlcfnames = lclist [ 'objects' ] [ 'lcfname' ] # if we're told to make a finder chart with the selected objects if field_fitsfile is not None and os . path . exists ( field_fitsfile ) : # get the RA and DEC of the matching objects matching_ra = lclist [ 'objects' ] [ racol ] [ finalfilterind ] matching_decl = lclist [ 'objects' ] [ declcol ] [ finalfilterind ] matching_postfix = [ ] if xmatchexternal is not None : matching_postfix . append ( 'xmatch_%s' % os . path . splitext ( os . path . basename ( xmatchexternal ) ) [ 0 ] ) if conesearch is not None : matching_postfix . append ( 'conesearch_RA%.3f_DEC%.3f_RAD%.5f' % tuple ( conesearch ) ) if columnfilters is not None : for cfi , cf in enumerate ( columnfilters ) : if cfi == 0 : matching_postfix . append ( 'filter_%s_%s_%s' % tuple ( cf . split ( '|' ) ) ) else : matching_postfix . append ( '_and_%s_%s_%s' % tuple ( cf . split ( '|' ) ) ) if len ( matching_postfix ) > 0 : matching_postfix = '-%s' % '_' . join ( matching_postfix ) else : matching_postfix = '' # next, we'll make a PNG plot for the finder finder_outfile = os . path . join ( os . path . dirname ( lc_catalog ) , '%s%s.png' % ( os . path . splitext ( os . path . basename ( lc_catalog ) ) [ 0 ] , matching_postfix ) ) finder_png = fits_finder_chart ( field_fitsfile , finder_outfile , wcsfrom = field_wcsfrom , scale = field_scale , stretch = field_stretch , colormap = field_colormap , findersize = field_findersize , overlay_ra = matching_ra , overlay_decl = matching_decl , overlay_pltopts = field_pltopts , field_zoomcontain = field_zoomcontain , grid = field_grid , gridcolor = field_gridcolor ) if finder_png is not None : LOGINFO ( 'generated a finder PNG ' 'with an object position overlay ' 'for this filtered LC list: %s' % finder_png ) # if copylcsto is not None, copy LCs over to it if copylcsto is not None : if not os . path . exists ( copylcsto ) : os . mkdir ( copylcsto ) if TQDM : lciter = tqdm ( filteredlcfnames ) else : lciter = filteredlcfnames LOGINFO ( 'copying matching light curves to %s' % copylcsto ) for lc in lciter : shutil . copy ( lc , copylcsto ) LOGINFO ( 'done. objects matching all filters: %s' % filteredobjectids . size ) if xmatchexternal and len ( ext_matching_objects ) > 0 : return filteredlcfnames , filteredobjectids , ext_matching_objects else : return filteredlcfnames , filteredobjectids
This is used to perform cone-search, cross-match, and column-filter
    operations on a light curve catalog generated by `make_lclist`.

    Uses the output of `make_lclist` above. This function returns a list of
    light curves matching various criteria specified by the `xmatchexternal`,
    `conesearch`, and `columnfilters kwargs`. Use this function to generate
    input lists for other lcproc functions,
    e.g. `lcproc.lcvfeatures.parallel_varfeatures`,
    `lcproc.periodfinding.parallel_pf`, and `lcproc.lcbin.parallel_timebin`,
    among others.

    The operations are applied in this order if more than one is specified:
    `xmatchexternal` -> `conesearch` -> `columnfilters`. All results from these
    operations are joined using a logical AND operation.

    Parameters
    ----------

    objectidcol : str
        This is the name of the object ID column in the light curve catalog.

    racol : str
        This is the name of the RA column in the light curve catalog.

    declcol : str
        This is the name of the Dec column in the light curve catalog.

    xmatchexternal : str or None
        If provided, this is the filename of a text file containing objectids,
        ras and decs to match the objects in the light curve catalog to by their
        positions.

    xmatchdistarcsec : float
        This is the distance in arcseconds to use when cross-matching to the
        external catalog in `xmatchexternal`.

    externalcolnums : sequence of int
        This a list of the zero-indexed column numbers of columns to extract
        from the external catalog file.

    externalcolnames : sequence of str
        This is a list of names of columns that will be extracted from the
        external catalog file. This is the same length as
        `externalcolnums`. These must contain the names provided as the
        `objectid`, `ra`, and `decl` column names so this function knows which
        column numbers correspond to those columns and can use them to set up
        the cross-match.

    externalcoldtypes : str
        This is a CSV string containing numpy dtype definitions for all columns
        listed to extract from the external catalog file. The number of dtype
        definitions should be equal to the number of columns to extract.

    externalcolsep : str or None
        The column separator to use when extracting columns from the external
        catalog file. If None, any whitespace between columns is used as the
        separator.

    externalcommentchar : str
        The character indicating that a line in the external catalog file is to
        be ignored.

    conesearch : list of float

        This is used to specify cone-search parameters. It should be a three
        element list:

        [center_ra_deg, center_decl_deg, search_radius_deg]

    conesearchworkers : int
        The number of parallel workers to launch for the cone-search operation.

    columnfilters : list of str

        This is a list of strings indicating any filters to apply on each column
        in the light curve catalog. All column filters are applied in the
        specified sequence and are combined with a logical AND operator. The
        format of each filter string should be:

        '<lc_catalog column>|<operator>|<operand>'

        where:

        - <lc_catalog column> is a column in the lc_catalog pickle file
        - <operator> is one of: 'lt', 'gt', 'le', 'ge', 'eq', 'ne', which
          correspond to the usual operators: <, >, <=, >=, ==, != respectively.
        - <operand> is a float, int, or string.

    field_fitsfile : str or None
        If this is not None, it should be the path to a FITS image containing
        the objects these light curves are for. If this is provided,
        `make_lclist` will use the WCS information in the FITS itself if
        `field_wcsfrom` is None (or from a WCS header file pointed to by
        `field_wcsfrom`) to obtain x and y pixel coordinates for all of the
        objects in the field. A finder chart will also be made using
        `astrobase.plotbase.fits_finder_chart` using the corresponding
        `field_scale`, `_stretch`, `_colormap`, `_findersize`, `_pltopts`,
        `_grid`, and `_gridcolors` kwargs for that function, reproduced here to
        enable customization of the finder chart plot.

    field_wcsfrom : str or None
        If `wcsfrom` is None, the WCS to transform the RA/Dec to pixel x/y will
        be taken from the FITS header of `fitsfile`. If this is not None, it
        must be a FITS or similar file that contains a WCS header in its first
        extension.

    field_scale : astropy.visualization.Interval object
        `scale` sets the normalization for the FITS pixel values. This is an
        astropy.visualization Interval object.
        See http://docs.astropy.org/en/stable/visualization/normalization.html
        for details on `scale` and `stretch` objects.

    field_stretch : astropy.visualization.Stretch object
        `stretch` sets the stretch function for mapping FITS pixel values to
        output pixel values. This is an astropy.visualization Stretch object.
        See http://docs.astropy.org/en/stable/visualization/normalization.html
        for details on `scale` and `stretch` objects.

    field_colormap : matplotlib Colormap object
        `colormap` is a matplotlib color map object to use for the output image.

    field_findersize : None or tuple of two ints
        If `findersize` is None, the output image size will be set by the NAXIS1
        and NAXIS2 keywords in the input `fitsfile` FITS header. Otherwise,
        `findersize` must be a tuple with the intended x and y size of the image
        in inches (all output images will use a DPI = 100).

    field_pltopts : dict
        `field_pltopts` controls how the overlay points will be plotted. This
        a dict with standard matplotlib marker, etc. kwargs as key-val pairs,
        e.g. 'markersize', 'markerfacecolor', etc. The default options make red
        outline circles at the location of each object in the overlay.

    field_grid : bool
        `grid` sets if a grid will be made on the output image.

    field_gridcolor : str
        `gridcolor` sets the color of the grid lines. This is a usual matplotib
        color spec string.

    field_zoomcontain : bool
        `field_zoomcontain` controls if the finder chart will be zoomed to
        just contain the overlayed points. Everything outside the footprint of
        these points will be discarded.

    copylcsto : str
        If this is provided, it is interpreted as a directory target to copy
        all the light curves that match the specified conditions.

    Returns
    -------

    tuple
        Returns a two elem tuple: (matching_object_lcfiles, matching_objectids)
        if conesearch and/or column filters are used. If `xmatchexternal` is
        also used, a three-elem tuple is returned: (matching_object_lcfiles,
        matching_objectids, extcat_matched_objectids).
This is used to perform cone - search cross - match and column - filter operations on a light curve catalog generated by make_lclist .