text: """Code that evaluates the quality of the SVM products generated by the drizzlepac package.
  
  The JSON files generated here can be converted directly into a Pandas DataFrame
  using the syntax:
  
  >>> import json
  >>> import pandas as pd
  >>> with open("<rootname>_astrometry_resids.json") as jfile:
  >>>     resids = json.load(jfile)
  >>> pdtab = pd.DataFrame(resids)
  
  These DataFrames can then be concatenated using:
  
  >>> allpd = pdtab.concat([pdtab2, pdtab3])
  
  where 'pdtab2' and 'pdtab3' are DataFrames generated from other datasets.  For
  more information on how to merge DataFrames, see
  
  https://pandas.pydata.org/pandas-docs/stable/user_guide/merging.html
  
  Visualization of these Pandas DataFrames with Bokeh can follow the example
  from:
  
  https://programminghistorian.org/en/lessons/visualizing-with-bokeh
  
  """
  
  # Standard library imports
  import argparse
  import collections
  from datetime import datetime
  import glob
  import json
  import math
  import os
  import pdb
  import pickle
  import sys
  import time
  
  # Related third party imports
  from astropy.coordinates import SkyCoord
  from astropy.io import ascii, fits
  from astropy.stats import sigma_clipped_stats
  from astropy.table import Table, Column
  from bokeh.layouts import row, column
  from bokeh.plotting import figure, output_file, save
  from bokeh.models import ColumnDataSource, Label
  from bokeh.models.tools import HoverTool
  from itertools import chain
  import numpy as np
  from photutils import DAOStarFinder
  from scipy import ndimage
  from scipy.spatial import KDTree
  
  # Local application imports
  from drizzlepac import util, wcs_functions
  from drizzlepac.haputils import hla_flag_filter
  from drizzlepac.haputils import catalog_utils
  from drizzlepac.haputils import astrometric_utils as au
  from drizzlepac.haputils import cell_utils
  import drizzlepac.haputils.comparison_utils as cu
  import drizzlepac.haputils.diagnostic_utils as du
  try:
      import drizzlepac.devutils.comparison_tools.compare_sourcelists as csl
  except ImportError:
      print("Unable to import compare_sourcelists! HAP vs. HLA sourcelist comparisons will be skipped.")
  from drizzlepac.haputils import read_hla_catalog
  from stsci.tools import logutil
  from stwcs import wcsutil
  from stwcs.wcsutil import HSTWCS
  
  from drizzlepac.haputils import quality_analysis as qa
  from drizzlepac.haputils.pandas_utils import PandasDFReader
  
  
  __taskname__ = 'svm_quality_analysis'
  
  MSG_DATEFMT = '%Y%j%H%M%S'
  SPLUNK_MSG_FORMAT = '%(asctime)s %(levelname)s src=%(name)s- %(message)s'
  log = logutil.create_logger(__name__, level=logutil.logging.NOTSET, stream=sys.stdout,
                              format=SPLUNK_MSG_FORMAT, datefmt=MSG_DATEFMT)
  # ----------------------------------------------------------------------------------------------------------------------
  
  
  def characterize_gaia_distribution(hap_obj, json_timestamp=None, json_time_since_epoch=None,
                                     log_level=logutil.logging.NOTSET):
      """Statistically describe distribution of GAIA sources in footprint.
  
      Computes and writes the file to a json file:
  
      - Number of GAIA sources
      - X centroid location
      - Y centroid location
      - X offset of centroid from image center
      - Y offset of centroid from image center
      - X standard deviation
      - Y standard deviation
      - minimum closest neighbor distance
      - maximum closest neighbor distance
      - mean closest neighbor distance
      - standard deviation of closest neighbor distances
  
      Parameters
      ----------
      hap_obj : drizzlepac.haputils.Product.FilterProduct
          hap product object to process
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      Returns
      -------
      Nothing
      """
      log.setLevel(log_level)
      log.info('\n\n*****     Begin Quality Analysis Test: characterize_gaia_distribution.     *****\n')
  
      # get table of GAIA sources in footprint
      gaia_table = generate_gaia_catalog(hap_obj, columns_to_remove=['mag', 'objID', 'GaiaID'])
  
      # return without creating a .json if no GAIA sources were found
      if not gaia_table:
          log.error("*** No GAIA sources were found. Characterization of GAIA distribution cannot be computed. "
                    "No json file will be produced.***")
          return
  
      # if log_level is either 'DEBUG' or 'NOTSET', write out GAIA sources to DS9 region file
      if log_level <= logutil.logging.DEBUG:
          reg_file = "{}_gaia_sources.reg".format(hap_obj.drizzle_filename[:-9])
          cname = gaia_table.colnames[0]
          orig_cname = cname
          gaia_table.rename_column(cname, "#{}".format(cname))  # Make reg file DS9-compatible
          gaia_table.write(reg_file, format='ascii.csv')
          log.debug("Wrote GAIA source RA and Dec positions to DS9 region file '{}'".format(reg_file))
          # Reset the column name back
          cname = gaia_table.colnames[0]
          gaia_table.rename_column(cname, orig_cname)
  
      # convert RA, Dec to image X, Y
      outwcs = HSTWCS(hap_obj.drizzle_filename + "[1]")
      x, y = outwcs.all_world2pix(gaia_table['RA'], gaia_table['DEC'], 1)
  
      # compute stats for the distribution
      centroid = [np.mean(x), np.mean(y)]
      centroid_offset = []
      for idx in range(0, 2):
          centroid_offset.append(outwcs.wcs.crpix[idx] - centroid[idx])
      std_dev = [np.std(x), np.std(y)]
  
      # Find straight-line distance to the closest neighbor for each GAIA source
      xys = np.array([x, y])
      xys = xys.reshape(len(x), 2)
      tree = KDTree(xys)
      neighborhood = tree.query(xys, 2)
      min_seps = np.empty([0])
      for sep_pair in neighborhood[0]:
          min_seps = np.append(min_seps, sep_pair[1])
  
      # add statistics to out_dict
      out_dict = collections.OrderedDict()
      out_dict["Number of GAIA sources"] = len(gaia_table)
      axis_list = ["X", "Y"]
      title_list = ["centroid", "offset", "standard deviation"]
      for item_value, item_title in zip([centroid, centroid_offset, std_dev], title_list):
          for axis_item in enumerate(axis_list):
              log.info("{} {} ({}): {}".format(axis_item[1], item_title, "pixels",
                                               item_value[axis_item[0]]))
              out_dict["{} {}".format(axis_item[1], item_title)] = item_value[axis_item[0]]
      min_sep_stats = [min_seps.min(), min_seps.max(), min_seps.mean(), min_seps.std()]
      min_sep_title_list = ["minimum neighbor distance",
                            "maximum neighbor distance",
                            "mean neighbor distance",
                            "standard deviation of neighbor distances"]
      for item_value, item_title in zip(min_sep_stats, min_sep_title_list):
          log.info("{} ({}): {}".format(item_title, "pixels", item_value))
          out_dict[item_title] = item_value
  
      # write catalog to HapDiagnostic-formatted .json file.
      diag_obj = du.HapDiagnostic(log_level=log_level)
      diag_obj.instantiate_from_hap_obj(hap_obj,
                                        data_source="{}.characterize_gaia_distribution".format(__taskname__),
                                        description="A statistical characterization of the distribution of "
                                                    "GAIA sources in image footprint",
                                        timestamp=json_timestamp, time_since_epoch=json_time_since_epoch)
      diag_obj.add_data_item(out_dict, "distribution characterization statistics",
                             descriptions={"Number of GAIA sources": "Number of GAIA sources in image footprint",
                                           "X centroid": "X centroid", "Y centroid": "Y centroid",
                                           "X offset": "X offset of centroid from image center",
                                           "Y offset": "Y offset of centroid from image center",
                                           "X standard deviation": "X standard deviation of mean offset",
                                           "Y standard deviation": "Y standard deviation of mean offset",
                                           "minimum neighbor distance": "distance of closest neighbor",
                                           "maximum neighbor distance": "distance of furthest neighbor",
                                           "mean neighbor distance": "mean distance of neighbors",
                                           "standard deviation of neighbor distances": "standard deviation of mean distance"},
                             units={"Number of GAIA sources": "unitless",
                                    "X centroid": "pixels", "Y centroid": "pixels",
                                    "X offset": "pixels",
                                    "Y offset": "pixels",
                                    "X standard deviation": "pixels",
                                    "Y standard deviation": "pixels",
                                    "minimum neighbor distance": "pixels",
                                    "maximum neighbor distance": "pixels",
                                    "mean neighbor distance": "pixels",
                                    "standard deviation of neighbor distances": "pixels"})
      diag_obj.write_json_file(hap_obj.drizzle_filename[:-9] + "_svm_gaia_distribution_characterization.json",
                               clobber=True)
  
  
  # ----------------------------------------------------------------------------------------------------------------------
  
  def compare_num_sources(drizzle_list, json_timestamp=None, json_time_since_epoch=None,
                          log_level=logutil.logging.NOTSET):
      """Determine the number of viable sources actually listed in SVM output catalogs.
  
      Parameters
      ----------
      drizzle_list: list of strings
          Drizzle files for the Total products which were mined to generate the output catalogs.
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      .. note:: This routine can be run either as a direct call from the hapsequencer.py routine,
      or it can invoked by a simple Python driver (or from within a Python session) by providing
      the names of the previously computed files as lists. The files must exist in the working directory.
      """
      log.setLevel(log_level)
      log.info('\n\n*****     Begin Quality Analysis Test: compare_num_sources.     *****\n')
  
      pnt_suffix = '_point-cat.ecsv'
      seg_suffix = '_segment-cat.ecsv'
  
      # Get a list of the catalog files (*.ecsv) in the current working directory
      catalog_files = []
      basepath = os.getcwd()
      for entry in os.listdir(basepath):
          if os.path.isfile(os.path.join(basepath, entry)):
              if entry.endswith("ecsv"):
                  catalog_files.append(entry)
  
      # Generate a separate JSON file for each detector
      # Drizzle filename example: hst_11665_06_wfc3_ir_total_ib4606_drz.fits
      # The filename is all lower-case by design.
      for drizzle_file in drizzle_list:
          if not os.path.exists(drizzle_file):
              log.warning("[compare_num_sources] Input {} not found. Skipping comparison.".format(drizzle_file))
              continue
          tokens = drizzle_file.split('_')
          detector = tokens[4]
  
          sources_dict = {'detector': detector, 'point': 0, 'segment': 0}
  
          # Construct the output JSON filename
          json_filename = drizzle_file[:-9] + '_svm_num_sources.json'
  
          # Construct catalog names for catalogs that should have been produced
          prefix = '_'.join(tokens[0:-1])
          cat_names = [prefix + pnt_suffix, prefix + seg_suffix]
  
          # If the catalog were actually produced, get the number of sources.
          # A catalog may not be produced because it was not requested, or there
          # was an error.  However, for the purposes of this program, it is OK
          # that no catalog was produced.
          for catalog in cat_names:
              if catalog in catalog_files:
  
                  # if the catalog exists, open it and find the number of sources string
                  num_sources = -1
                  cat_type = ""
                  file = open(catalog, 'r')
                  for line in file:
                      sline = line.strip()
  
                      # All the comments are grouped at the start of the file. When
                      # the first non-comment line is found, there is no need to look further.
                      if not sline.startswith('#'):
                          log.info("Number of sources not reported in Catalog: {}.".format(catalog))
                          break
  
                      # When the matching comment line is found, get the value.
                      if sline.find('Number of sources') != -1:
                          num_sources = sline.split(' ')[-1][0:-1]
                          log.info("Catalog: {} Number of sources: {}.".format(catalog, num_sources))
                          break
  
                  cat_type = 'point' if catalog.find("point") != -1 else 'segment'
                  sources_dict[cat_type] = int(num_sources)
  
          # Set up the diagnostic object and write out the results
          diagnostic_obj = du.HapDiagnostic()
          diagnostic_obj.instantiate_from_fitsfile(drizzle_file,
                                                   data_source="{}.compare_num_sources".format(__taskname__),
                                                   description="Number of sources in Point and Segment "
                                                               "catalogs",
                                                   timestamp=json_timestamp,
                                                   time_since_epoch=json_time_since_epoch)
          diagnostic_obj.add_data_item(sources_dict, 'number_of_sources',
                                       descriptions={'detector': 'Detector in use',
                                                     'point': 'Number of detected sources in Point catalog',
                                                     'segment': 'Number of detected sources in Segmented catalog'},
                                       units={'detector': 'unitless', 'point': 'unitless', 'segment': 'unitless'})
          diagnostic_obj.write_json_file(json_filename)
          log.info("Generated quality statistics (number of sources) as {}.".format(json_filename))
  
          # Clean up
          del diagnostic_obj
  
  # ------------------------------------------------------------------------------------------------------------
  
  
  def compare_ra_dec_crossmatches(hap_obj, json_timestamp=None, json_time_since_epoch=None,
                                  log_level=logutil.logging.NOTSET):
      """Compare the equatorial coordinates of cross-matches sources between the Point and Segment catalogs.
      The results .json file contains the following information:
  
          - image header information
          - cross-match details (input catalog lengths, number of cross-matched sources, coordinate system)
          - catalog containing RA and dec values of cross-matched point catalog sources
          - catalog containing RA and dec values of cross-matched segment catalog sources
          - catalog containing MagAp1 and MagAp2 values of cross-matched point catalog sources
          - catalog containing MagAp1 and MagAp2 values of cross-matched segment catalog sources
          - Statistics describing the on-sky separation of the cross-matched point and segment catalogs
          (non-clipped and sigma-clipped mean, median and standard deviation values)
          - Statistics describing the mean, std, and median for the differences of magnitudes of the
            cross-matched point and segment catalogs
  
      Parameters
      ----------
      hap_obj : drizzlepac.haputils.Product.FilterProduct
          hap filter product object to process
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      Returns
      --------
      nothing.
      """
      log.setLevel(log_level)
      log.info('\n\n*****     Begin Quality Analysis Test: compare_ra_dec_crossmatches.     *****\n')
  
      sl_names = [hap_obj.point_cat_filename, hap_obj.segment_cat_filename]
      img_names = [hap_obj.drizzle_filename, hap_obj.drizzle_filename]
      good_flag_sum = 255  # all bits good
  
      for name in sl_names:
          if name.rstrip == '' or not os.path.exists(name):
              log.warning("[compare_ra_dec_crossmatches] Catalog {} Missing!  No comparison can be made.".format(name))
              return
  
      diag_obj = du.HapDiagnostic(log_level=log_level)
      diag_obj.instantiate_from_hap_obj(hap_obj,
                                        data_source="{}.compare_ra_dec_crossmatches".format(__taskname__),
                                        description="matched point and segment catalog RA and Dec values",
                                        timestamp=json_timestamp,
                                        time_since_epoch=json_time_since_epoch)
      json_results_dict = collections.OrderedDict()
      # add reference and comparison catalog filenames as header elements
      json_results_dict["point catalog filename"] = sl_names[0]
      json_results_dict["segment catalog filename"] = sl_names[1]
  
      # 1: Read in sourcelists files into astropy table or 2-d array so that individual columns from each
      # sourcelist can be easily accessed later in the code.
      point_data, seg_data = cu.slFiles2dataTables(sl_names)
      log.info("Valid point data columns:   {}".format(list(point_data.keys())))
      log.info("Valid segment data columns: {}".format(list(seg_data.keys())))
      log.info("\n")
      log.info("Data columns to be compared:")
      columns_to_compare = list(set(point_data.keys()).intersection(set(seg_data.keys())))
      for listItem in sorted(columns_to_compare):
          log.info(listItem)
      log.info("\n")
      # 2: Run starmatch_hist to get list of matched sources common to both input sourcelists
      sl_lengths = [len(point_data['RA']), len(seg_data['RA'])]
      json_results_dict['point catalog length'] = sl_lengths[0]
      json_results_dict['segment catalog length'] = sl_lengths[1]
      matching_lines_ref, matching_lines_img = cu.getMatchedLists(sl_names, img_names, sl_lengths,
                                                                  log_level=log_level)
      json_results_dict['number of cross-matches'] = len(matching_lines_ref)
  
      # Report number and percentage of the total number of detected ref and comp sources that were matched
      log.info("Cross-matching results")
      log.info(
          "Point sourcelist:  {} of {} total sources cross-matched ({}%)".format(len(matching_lines_ref),
                                                                                 sl_lengths[0],
                                                                                 100.0 *
                                                                                 (float(len(matching_lines_ref))
                                                                                  / float(sl_lengths[0]))))
      log.info(
          "Segment sourcelist: {} of {} total sources cross-matched ({}%)".format(len(matching_lines_img),
                                                                                  sl_lengths[1],
                                                                                  100.0 *
                                                                                  (float(
                                                                                      len(matching_lines_img))
                                                                                   / float(sl_lengths[1]))))
      # return without creating a .json if no cross-matches are found
      if len(matching_lines_ref) == 0 or len(matching_lines_img) == 0:
          log.error("*** No matching sources were found. Comparisons cannot be computed. "
                    "No json file will be produced.***")
          return
      # 2: Create masks to remove missing values or values not considered "good" according to user-specified
      # good bit values
      # 2a: create mask that identifies lines any value from any column is missing
      missing_mask = cu.mask_missing_values(point_data, seg_data, matching_lines_ref, matching_lines_img,
                                            columns_to_compare)
      # 2b: create mask based on flag values
      matched_values = cu.extractMatchedLines("FLAGS", point_data, seg_data, matching_lines_ref,
                                              matching_lines_img)
  
      bitmask = cu.make_flag_mask(matched_values, good_flag_sum, missing_mask)
  
      matched_values_ra = cu.extractMatchedLines("RA", point_data, seg_data, matching_lines_ref,
                                                 matching_lines_img, bitmask=bitmask)
      matched_values_dec = cu.extractMatchedLines("DEC", point_data, seg_data, matching_lines_ref,
                                                  matching_lines_img, bitmask=bitmask)
      matched_values_magap1 = cu.extractMatchedLines("MAGNITUDE1", point_data, seg_data, matching_lines_ref,
                                                     matching_lines_img, bitmask=bitmask)
      matched_values_magap2 = cu.extractMatchedLines("MAGNITUDE2", point_data, seg_data, matching_lines_ref,
                                                     matching_lines_img, bitmask=bitmask)
  
      if matched_values_ra.shape[1] > 0 and matched_values_ra.shape[1] == matched_values_dec.shape[1]:
          # get coordinate system type from fits headers
  
          point_frame = fits.getval(img_names[0], "radesys", ext=('sci', 1)).lower()
          seg_frame = fits.getval(img_names[1], "radesys", ext=('sci', 1)).lower()
          # Add 'ref_frame' and 'comp_frame" values to header so that will SkyCoord() execute OK
          json_results_dict["point frame"] = point_frame
          json_results_dict["segment frame"] = seg_frame
  
          # convert reference and comparison RA/Dec values into SkyCoord objects
          matched_values_point = SkyCoord(matched_values_ra[0, :], matched_values_dec[0, :], frame=point_frame,
                                          unit="deg")
          matched_values_seg = SkyCoord(matched_values_ra[1, :], matched_values_dec[1, :], frame=seg_frame,
                                        unit="deg")
          # convert to ICRS coord system
          if point_frame != "icrs":
              matched_values_point = matched_values_point.icrs
          if seg_frame != "icrs":
              matched_values_seg = matched_values_seg.icrs
  
          # compute on-sky separations in arcseconds
          sep = matched_values_seg.separation(matched_values_point).arcsec
  
          # Compute and store statistics on separations
          sep_stat_dict = collections.OrderedDict()
          sep_stat_dict["Non-clipped min"] = np.min(sep)
          sep_stat_dict["Non-clipped max"] = np.max(sep)
          sep_stat_dict["Non-clipped mean"] = np.mean(sep)
          sep_stat_dict["Non-clipped median"] = np.median(sep)
          sep_stat_dict["Non-clipped standard deviation"] = np.std(sep)
          sigma = 3
          maxiters = 3
          clipped_stats = sigma_clipped_stats(sep, sigma=sigma, maxiters=maxiters)
          sep_stat_dict["{}x{} sigma-clipped mean".format(maxiters, sigma)] = clipped_stats[0]
          sep_stat_dict["{}x{} sigma-clipped median".format(maxiters, sigma)] = clipped_stats[1]
          sep_stat_dict["{}x{} sigma-clipped standard deviation".format(maxiters, sigma)] = clipped_stats[2]
  
          #
          # Compute statistics on the photometry differences
          #
  
          # Compute the differences (Point - Segment)
          delta_phot_magap1 = np.subtract(matched_values_magap1[0], matched_values_magap1[1])
          delta_phot_magap2 = np.subtract(matched_values_magap2[0], matched_values_magap2[1])
  
          # Compute some basic statistics: mean difference and standard deviation, and median difference
          phot_stat_dict = collections.OrderedDict()
          phot_stat_dict["mean_dmagap1"] = np.mean(delta_phot_magap1)
          phot_stat_dict["std_dmagap1"] = np.std(delta_phot_magap1)
          phot_stat_dict["median_dmagap1"] = np.median(delta_phot_magap1)
          phot_stat_dict["mean_dmagap2"] = np.mean(delta_phot_magap2)
          phot_stat_dict["std_dmagap2"] = np.std(delta_phot_magap2)
          phot_stat_dict["median_dmagap2"] = np.median(delta_phot_magap2)
  
          # Create output catalogs for json file
          out_cat_point = Table([matched_values_ra[0], matched_values_dec[0], matched_values_magap1[0],
                                matched_values_magap2[0]], names=("Right ascension", "Declination",
                                "MagAp1", "MagAp2"))
          out_cat_seg = Table([matched_values_ra[1], matched_values_dec[1], sep, matched_values_magap1[1],
                              matched_values_magap2[1]], names=("Right ascension", "Declination",
                              "Separation",
                              "MagAp1", "MagAp2"))
          for table_item in [out_cat_point, out_cat_seg]:
              for col_name in ["Right ascension", "Declination"]:
                  table_item[col_name].unit = "degrees"  # Add correct units
              for col_name in ["MagAp1", "MagAp2"]:
                  table_item[col_name].unit = "ABMag"  # Add correct units
          out_cat_seg['Separation'].unit = "arcseconds"
  
          # add various data items to diag_obj
          diag_obj.add_data_item(json_results_dict, "Cross-match details",
                                 descriptions={"point catalog filename": "ECSV point catalog filename",
                                               "segment catalog filename": "ECSV segment catalog filename",
                                               "point catalog length": "Number of entries in point catalog",
                                               "segment catalog length": "Number of entries in segment catalog",
                                               "number of cross-matches": "Number of cross-matches between point and segment catalogs",
                                               "point frame": "Coordinate reference frame",
                                               "segment frame": "Coordinate reference frame"},
                                 units={"point catalog filename": "unitless",
                                        "segment catalog filename": "unitless",
                                        "point catalog length": "unitless",
                                        "segment catalog length": "unitless",
                                        "number of cross-matches": "unitless",
                                        "point frame": "unitless",
                                        "segment frame": "unitless"})
          diag_obj.add_data_item(out_cat_point, "Cross-matched point catalog",
                                 descriptions={"Right ascension": "ICRS Right ascension",
                                               "Declination": "ICRS Declination",
                                               "MagAp1": "Magnitude Aperture 1",
                                               "MagAp2": "Magnitude Aperture 2"},
                                 units={"Right ascension": "degrees", "Declination": "degrees",
                                        "MagAp1": "ABMag", "MagAp2": "ABMag"})
          diag_obj.add_data_item(out_cat_seg, "Cross-matched segment catalog",
                                 descriptions={"Right ascension": "ICRS Right ascension",
                                               "Declination": "ICRS Declination",
                                               "Separation": "Segment minus Point on-sky coordinate separation",
                                               "MagAp1": "Magnitude Aperture 1",
                                               "MagAp2": "Magnitude Aperture 2"},
                                 units={"Right ascension": "degrees", "Declination": "degrees",
                                        "Separation": "arcseconds",
                                        "MagAp1": "ABMag", "MagAp2": "ABMag"})
          diag_obj.add_data_item(sep_stat_dict, "Segment - point on-sky separation statistics",
                                 descriptions={"Non-clipped min": "Non-clipped min difference",
                                               "Non-clipped max": "Non-clipped max difference",
                                               "Non-clipped mean": "Non-clipped mean difference",
                                               "Non-clipped median": "Non-clipped median difference",
                                               "Non-clipped standard deviation": "Non-clipped standard deviation of differences",
                                               "3x3 sigma-clipped mean": "3x3 sigma-clipped mean difference",
                                               "3x3 sigma-clipped median": "3x3 sigma-clipped median difference",
                                               "3x3 sigma-clipped standard deviation": "3x3 sigma-clipped standard deviation of differences"},
                                 units={"Non-clipped min": "arcseconds", "Non-clipped max": "arcseconds",
                                        "Non-clipped mean": "arcseconds", "Non-clipped median": "arcseconds",
                                        "Non-clipped standard deviation": "arcseconds",
                                        "3x3 sigma-clipped mean": "arcseconds", "3x3 sigma-clipped median": "arcseconds",
                                        "3x3 sigma-clipped standard deviation": "arcseconds"})
          diag_obj.add_data_item(phot_stat_dict, "Delta_Photometry",
                                 descriptions={'mean_dmagap1': 'dMagAp1_Mean_Differences(Point-Segment)',
                                               'std_dmagap1': 'dMagAp1_StdDev_of_Mean_Differences',
                                               'median_dmagap1': 'dMagAp1_Median_Differences(Point-Segment)',
                                               'mean_dmagap2': 'dMagAp2_Mean_Differences(Point-Segment)',
                                               'std_dmagap2': 'dMagAp2_StdDev_of_Mean_Differences',
                                               'median_dmagap2': 'dMagAp2_Median_Differences(Point-Segment)'},
                                 units={'mean_dmagap1': 'ABMag',
                                        'std_dmagap1': 'ABMag',
                                        'median_dmagap1': 'ABMag',
                                        'mean_dmagap2': 'ABMag',
                                        'std_dmagap2': 'ABMag',
                                        'median_dmagap2': 'ABMag'})
  
          # write everything out to the json file
          json_filename = hap_obj.drizzle_filename[:-9]+"_svm_point_segment_crossmatch.json"
          diag_obj.write_json_file(json_filename, clobber=True)
      else:
          log.warning("Point vs. segment catalog cross match test could not be performed.")
  
  # ------------------------------------------------------------------------------------------------------------
  
  
  def compare_interfilter_crossmatches(total_obj_list, json_timestamp=None, json_time_since_epoch=None,
                                       log_level=logutil.logging.NOTSET):
      """Compare X, Y positions of cross matched sources from different filter-level products
  
      Parameters
      ----------
      total_obj_list : list
          list of drizzlepac.haputils.Product.TotalProduct objects to process
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      Returns
      -------
      Nothing
      """
      # Initiate logging!
      log.setLevel(log_level)
  
      log.info('\n\n*****     Begin Quality Analysis Test: compare_interfilter_crossmatches.     *****\n')
      # Check to make sure there's at last 2 filter-level products. If not, return.
      num_filter_prods = 0
      for total_obj in total_obj_list:
          num_filter_prods += len(total_obj.fdp_list)
      if num_filter_prods == 1:
          log.error("At least 2 filter-level products are necessary to perform analysis of inter-filter "
                    "cross matched sources. Only 1 filter-level product was found. ")
          return
      elif num_filter_prods == 0:
          log.error("At least 2 filter-level products are necessary to perform analysis of inter-filter "
                    "cross matched sources. No filter-level products were found. ")
          return
      else:
          log.info("Found {} filter-level products for use in analysis of inter-filter cross matched source"
                   " positions".format(num_filter_prods))
  
          filtobj_dict = {}
          xmatch_ref_imgname = None
          max_sources = 0
          ctr = 1
          for total_obj in total_obj_list:
              for filt_obj in total_obj.fdp_list:
                  log.info("{} {}: {} {}".format(">"*20, ctr, filt_obj.drizzle_filename, "<"*20))
                  filtobj_dict[filt_obj.drizzle_filename] = find_hap_point_sources(filt_obj, log_level=log_level)
                  n_sources = len(filtobj_dict[filt_obj.drizzle_filename]['sources'])
                  log.info("Identified {} sources in {}".format(n_sources, filt_obj.drizzle_filename))
                  if n_sources > max_sources:
                      max_sources = n_sources
                      xmatch_ref_imgname = filt_obj.drizzle_filename
                      xmatch_ref_catname = filtobj_dict[filt_obj.drizzle_filename]['cat_name']
                  ctr += 1
                  log.info("")
          log.info("Crossmatch reference image {} contains {} sources.".format(xmatch_ref_imgname, max_sources))
  
          # Perform coord transform and write temp cat files for cross match
          temp_cat_file_list = []
          log.info("")
          for imgname in filtobj_dict.keys():
              filtobj_dict[imgname] = transform_coords(filtobj_dict[imgname],
                                                       xmatch_ref_imgname,
                                                       log_level=log_level)
  
              filtobj_dict[imgname]["sources"].write(filtobj_dict[imgname]['cat_name'], format="ascii.ecsv")
              log.info("Wrote temporary source catalog {}".format(filtobj_dict[imgname]['cat_name']))
  
              # write out ds9 region files if log level is 'debug'
              if log_level == logutil.logging.DEBUG:
                  out_reg_stuff = {"xyorig": ['xcentroid', 'ycentroid'],
                                   "rd": ['ra', 'dec'],
                                   "xyref": ['xcentroid_ref', 'ycentroid_ref']}
                  for reg_type in out_reg_stuff.keys():
                      reg_table = filtobj_dict[imgname]["sources"].copy()
                      reg_table.keep_columns(out_reg_stuff[reg_type])
                      reg_table.rename_column(out_reg_stuff[reg_type][0], "#"+out_reg_stuff[reg_type][0])
                      reg_filename = filtobj_dict[imgname]['cat_name'].replace(".ecsv",
                                                                               "_{}_all.reg".format(reg_type))
                      reg_table.write(reg_filename, format='ascii.csv')
                      log.debug("wrote region file {}".format(reg_filename))
                      del reg_filename
  
          # Perform cross-match based on X, Y coords
          for imgname in filtobj_dict.keys():
              if imgname != xmatch_ref_imgname:
                  log.info(" ")
                  xmatch_comp_imgname = imgname
                  xmatch_comp_catname = filtobj_dict[imgname]['cat_name']
                  filtername_ref = filtobj_dict[xmatch_ref_imgname]['filt_obj'].filters
                  filtername_comp = filtobj_dict[imgname]['filt_obj'].filters
  
                  # Perform crossmatching; get lists of crossmatched sources in the reference and comparison
                  log.info("{} Crossmatching {} -> {} {}".format(">" * 20, xmatch_comp_catname,
                                                                 xmatch_ref_catname, "<" * 20))
                  sl_names = [xmatch_ref_catname, xmatch_comp_catname]
                  img_names = [xmatch_ref_imgname, xmatch_comp_imgname]
                  sl_lengths = [max_sources, len(filtobj_dict[xmatch_comp_imgname]["sources"])]
                  matching_lines_ref, matching_lines_comp = cu.getMatchedLists(sl_names, img_names, sl_lengths,
                                                                               log_level=log_level)
  
                  # Report number and percentage of the total number of detected ref and comp sources that were
                  # matched
  
                  xmresults = []
                  xmresults.append("Reference sourcelist:  {} of {} total reference sources ({}%) cross-matched.".format(len(matching_lines_ref),
                                                                                                                         sl_lengths[0],
                                                                                                                         100.0 * (float(len(matching_lines_ref)) / float(sl_lengths[0]))))
                  xmresults.append("Comparison sourcelist: {} of {} total comparison sources ({}%) cross-matched.".format(len(matching_lines_comp),
                                                                                                                          sl_lengths[1],
                                                                                                                          100.0 * (float(len(matching_lines_comp)) / float(sl_lengths[1]))))
                  # display crossmatch results
                  padding = math.ceil((max([len(xmresults[0]), len(xmresults[1])]) - 18) / 2)
                  log.info("{}Crossmatch results".format(" "*padding))
                  for item in xmresults:
                      log.info(item)
                  log.info("")
                  if len(matching_lines_ref) > 0:
                      # instantiate diagnostic object to store test results for eventual .json file output
                      diag_obj = du.HapDiagnostic(log_level=log_level)
                      diag_obj.instantiate_from_hap_obj(filtobj_dict[xmatch_comp_imgname]['filt_obj'],
                                                        data_source="{}.compare_interfilter_crossmatches".format(__taskname__),
                                                        description="Interfilter cross-matched comparison and reference catalog X and Y values",
                                                        timestamp=json_timestamp,
                                                        time_since_epoch=json_time_since_epoch)
                      json_results_dict = collections.OrderedDict()
                      json_results_dict["reference image name"] = xmatch_ref_imgname
                      json_results_dict["comparison image name"] = xmatch_comp_imgname
                      json_results_dict['reference catalog length'] = sl_lengths[0]
                      json_results_dict['comparison catalog length'] = sl_lengths[1]
                      json_results_dict['number of cross-matches'] = len(matching_lines_ref)
                      json_results_dict['percent of all identified reference sources crossmatched'] = 100.0 * (
                              float(len(matching_lines_ref)) / float(sl_lengths[0]))
                      json_results_dict['percent of all identified comparison sources crossmatched'] = 100.0 * (
                                  float(len(matching_lines_comp)) / float(sl_lengths[1]))
                      json_results_dict['reference image platescale'] = filtobj_dict[xmatch_ref_imgname]['filt_obj'].meta_wcs.pscale
  
                      # store cross-match details
                      diag_obj.add_data_item(json_results_dict, "Interfilter cross-match details",
                                             descriptions={
                                                 "reference image name": "Crossmatch reference image name",
                                                 "comparison image name": "Crossmatch comparison image name",
                                                 "reference catalog length": "Number of entries in point catalog",
                                                 "comparison catalog length": "Number of entries in segment catalog",
                                                 "number of cross-matches": "Number of cross-matches between point and segment catalogs",
                                                 "percent of all identified reference sources crossmatched": "percent of all identified reference sources crossmatched",
                                                 "percent of all identified comparison sources crossmatched": "percent of all identified comparison sources crossmatched",
                                                 "reference image platescale": "Platescale of the crossmatch reference image"},
                                             units={"reference image name": "unitless",
                                                    "comparison image name": "unitless",
                                                    "reference catalog length": "unitless",
                                                    "comparison catalog length": "unitless",
                                                    "number of cross-matches": "unitless",
                                                    "percent of all identified reference sources crossmatched": "unitless",
                                                    "percent of all identified comparison sources crossmatched": "unitless",
                                                    "reference image platescale": "arcseconds/pixel"})
  
                      # Generate tables containing just "xcentroid_ref" and "ycentroid_ref" columns with only
                      # the cross-matched reference sources
                      matched_ref_coords = filtobj_dict[xmatch_ref_imgname]["sources"].copy()
                      matched_ref_coords.keep_columns(['xcentroid_ref', 'ycentroid_ref'])
  
                      matched_ref_coords = matched_ref_coords[matching_lines_ref]
  
                      # store reference matched sources catalog
                      diag_obj.add_data_item(matched_ref_coords, "Interfilter cross-matched reference catalog",
                                             descriptions={"xcentroid_ref": "xcentroid_ref",
                                                           "ycentroid_ref": "ycentroid_ref"},
                                             units={"xcentroid_ref": "pixels", "ycentroid_ref": "pixels"})
  
                      # write out ds9 region files if log level is 'debug'
                      if log_level == logutil.logging.DEBUG:
                          reg_filename = "{}_{}_ref_matches.reg".format(filtername_comp, filtername_ref)
                          matched_ref_coords.write(reg_filename, format='ascii.csv')
                          log.debug("wrote region file {}".format(reg_filename))
                          del reg_filename
  
                      # Generate tables containing just "xcentroid_ref" and "ycentroid_ref" columns with only
                      # the cross-matched comparison sources
                      matched_comp_coords = filtobj_dict[imgname]["sources"].copy()
                      matched_comp_coords.keep_columns(['xcentroid_ref', 'ycentroid_ref'])
                      matched_comp_coords = matched_comp_coords[matching_lines_comp]
  
                      # store comparison matched sources catalog
                      diag_obj.add_data_item(matched_comp_coords,
                                             "Interfilter cross-matched comparison catalog",
                                             descriptions={"xcentroid_ref": "xcentroid_ref",
                                                           "ycentroid_ref": "ycentroid_ref"},
                                             units={"xcentroid_ref": "pixels", "ycentroid_ref": "pixels"})
  
                      # write out ds9 region file if log level is 'debug'
                      if log_level == logutil.logging.DEBUG:
                          reg_filename = "{}_{}_comp_matches.reg".format(filtername_comp, filtername_ref)
                          matched_comp_coords.write(reg_filename, format='ascii.csv')
                          log.debug("wrote region file {}".format(reg_filename))
                          del reg_filename
  
                      # compute statistics
                      xy_separations_table = Table()
                      for colname in ["xcentroid_ref", "ycentroid_ref"]:
                          sep = matched_comp_coords[colname] - matched_ref_coords[colname]
  
                          # Add column of comp-ref differences to table
                          xy_separations_table['delta_{}'.format(colname)] = sep
  
                          # Compute and store statistics on separations
                          sep_stat_dict = collections.OrderedDict()
                          sep_stat_dict["Non-clipped min"] = np.min(sep)
                          sep_stat_dict["Non-clipped max"] = np.max(sep)
                          sep_stat_dict["Non-clipped mean"] = np.mean(sep)
                          sep_stat_dict["Non-clipped median"] = np.median(sep)
                          sep_stat_dict["Non-clipped standard deviation"] = np.std(sep)
                          sigma = 3
                          maxiters = 3
                          clipped_stats = sigma_clipped_stats(sep, sigma=sigma, maxiters=maxiters)
                          sep_stat_dict["{}x{} sigma-clipped mean".format(maxiters, sigma)] = clipped_stats[0]
                          sep_stat_dict["{}x{} sigma-clipped median".format(maxiters, sigma)] = clipped_stats[1]
                          sep_stat_dict["{}x{} sigma-clipped standard deviation".format(maxiters, sigma)] = clipped_stats[2]
                          sep_stat_dict["Reference image platescale"] = filtobj_dict[xmatch_ref_imgname][
                              'filt_obj'].meta_wcs.pscale
  
                          # Store statistics as new data section
                          diag_obj.add_data_item(sep_stat_dict, "Interfilter cross-matched {} comparison - reference separation statistics".format(colname),
                                                 descriptions={"Non-clipped min": "Non-clipped min difference",
                                                               "Non-clipped max": "Non-clipped max difference",
                                                               "Non-clipped mean": "Non-clipped mean difference",
                                                               "Non-clipped median": "Non-clipped median difference",
                                                               "Non-clipped standard deviation": "Non-clipped standard deviation of differences",
                                                               "3x3 sigma-clipped mean": "3x3 sigma-clipped mean difference",
                                                               "3x3 sigma-clipped median": "3x3 sigma-clipped median difference",
                                                               "3x3 sigma-clipped standard deviation": "3x3 sigma-clipped standard deviation of differences",
                                                               "Reference image platescale": "Platescale of the crossmatch reference image"},
                                                 units={"Non-clipped min": "pixels",
                                                        "Non-clipped max": "pixels",
                                                        "Non-clipped mean": "pixels",
                                                        "Non-clipped median": "pixels",
                                                        "Non-clipped standard deviation": "pixels",
                                                        "3x3 sigma-clipped mean": "pixels",
                                                        "3x3 sigma-clipped median": "pixels",
                                                        "3x3 sigma-clipped standard deviation": "pixels",
                                                        "Reference image platescale": "arcseconds/pixel"})
  
                          # display stats
                          padding = math.ceil((74 - 61) / 2)
                          log.info("{}Interfilter cross-matched {} comparison - reference separation statistics".format(" "*padding, colname))
                          max_str_length = 0
                          for stat_key in sep_stat_dict.keys():
                              if len(stat_key) > max_str_length:
                                  max_str_length = len(stat_key)
                          for stat_key in sep_stat_dict.keys():
                              padding = " " * (max_str_length - len(stat_key))
                              log.info("{}{}: {} {}".format(padding, stat_key, sep_stat_dict[stat_key],
                                                            diag_obj.out_dict['data']["Interfilter cross-matched {} comparison - reference separation statistics".format(colname)]['units'][stat_key]))
                          log.info("")
  
                      # store separations table
                      diag_obj.add_data_item(xy_separations_table,
                                             "Interfilter cross-matched comparison - reference separations",
                                             descriptions={"delta_x": "delta_x",
                                                           "delta_y": "delta_y"},
                                             units={"delta_x": "pixels", "delta_y": "pixels"})
  
                      # write everything out to the json file
                      json_filename = filtobj_dict[xmatch_comp_imgname]['filt_obj'].drizzle_filename[
                                      :-9] + "_svm_interfilter_crossmatch.json"
                      diag_obj.write_json_file(json_filename, clobber=True)
                  else:
                      filt_names = "{} - {}".format(filtername_comp, filtername_ref)
                      log.warning("{} interfilter cross match test could not be performed.".format(filt_names))
  
      # Housekeeping. Delete the *_point-cat-fxm.ecsv files created for cross-matching, and the
      # filtobj_dict dictionary
      log.info("")
      for imgname in filtobj_dict.keys():
          log.info("removing temporary catalog file {}".format(filtobj_dict[imgname]['cat_name']))
          os.remove(filtobj_dict[imgname]['cat_name'])
      del filtobj_dict
  
  
  # ------------------------------------------------------------------------------------------------------------
  
  
  def transform_coords(filtobj_subdict, xmatch_ref_imgname, log_level=logutil.logging.NOTSET):
      """Transform comparison image frame of reference x, y coords to RA and dec, then back to x, y coords in
      the cross-match reference image's frame of reference
  
      Parameters
      ----------
      filtobj_subdict : dict
          dictionary containing a drizzlepac.haputils.Product.FilterProduct and corresponding source catalog
          generated earlier in the run by **find_hap_point_sources**.
  
      xmatch_ref_imgname : str
          name of the reference image whose WCS info will be used to convert RA and dec values into x, y coords
          in a common frame
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      Returns
      -------
      filtobj_dict : dict
          the input **filtobj_dict** dictionary with the 'sources' table updated to include freshly computed ra,
          dec, x_centroid_ref, and y_centroid_ref columns
      """
      # Initiate logging!
      log.setLevel(log_level)
  
      # 1: stack up xcentroid and ycentorid columns from sources table
      xy_centroid_values = np.stack((filtobj_subdict['sources']['xcentroid'],
                                     filtobj_subdict['sources']['ycentroid']), axis=1)
  
      # 2: perform coordinate transforms.
      origin = 0
      fits_exten = "[1]"
      imgname = filtobj_subdict['filt_obj'].drizzle_filename
  
      # 2a: perform xcentroid, ycentroid -> ra, dec transform
      ra_dec_values = hla_flag_filter.xytord(xy_centroid_values, imgname, fits_exten, origin=origin)
  
      # 2b: perform ra, dec -> x_centroid_ref, y_centroid_ref transform
      xy_centroid_ref_values = hla_flag_filter.rdtoxy(ra_dec_values, xmatch_ref_imgname, fits_exten,
                                                      origin=origin)
  
      # 3: add new columns to filtobj_subdict['sources'] table
      title_data_dict = collections.OrderedDict()
      title_data_dict["ra"] = ra_dec_values[:, 0]
      title_data_dict["dec"] = ra_dec_values[:, 1]
      title_data_dict["xcentroid_ref"] = xy_centroid_ref_values[:, 0]
      title_data_dict["ycentroid_ref"] = xy_centroid_ref_values[:, 1]
      col_ctr = 3
      for col_name in title_data_dict.keys():
          col_to_add = Column(name=col_name, data=title_data_dict[col_name], dtype=np.float64)
          filtobj_subdict['sources'].add_column(col_to_add, index=col_ctr)
          col_ctr += 1
  
      return filtobj_subdict
  
  
  def find_gaia_sources(hap_obj, json_timestamp=None, json_time_since_epoch=None,
                        log_level=logutil.logging.NOTSET):
      """Creates a catalog of all GAIA sources in the footprint of a specified HAP final product image, and
      stores the GAIA object catalog as a hap diagnostic json file. The catalog contains RA, Dec and magnitude
      of each identified source. The catalog is sorted in descending order by brightness.
  
      Parameters
      ----------
      hap_obj : drizzlepac.haputils.Product.TotalProduct, drizzlepac.haputils.Product.FilterProduct, or
          drizzlepac.haputils.Product.ExposureProduct, depending on input.
          hap product object to process
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      Returns
      -------
      Nothing.
      """
      log.setLevel(log_level)
      log.info('\n\n*****     Begin Quality Analysis Test: find_gaia_sources.     *****\n')
  
      gaia_table = generate_gaia_catalog(hap_obj, columns_to_remove=['objID', 'GaiaID'])
  
      if not gaia_table:
          log.error("*** No GAIA sources were found. Finding of GAIA sources cannot be computed. "
                    "No json file will be produced.***")
          return
  
      # write catalog to HapDiagnostic-formatted .json file.
      diag_obj = du.HapDiagnostic(log_level=log_level)
      diag_obj.instantiate_from_hap_obj(hap_obj,
                                        data_source="{}.find_gaia_sources".format(__taskname__),
                                        description="A table of GAIA sources in image footprint",
                                        timestamp=json_timestamp,
                                        time_since_epoch=json_time_since_epoch)
      diag_obj.add_data_item(gaia_table, "GAIA sources", descriptions={"RA": "Right Ascension", "DEC": "Declination", "mag": "AB Magnitude"}, units={"RA": "degrees", "DEC": "degrees", "mag": "unitless"})  # write catalog of identified GAIA sources
      diag_obj.add_data_item({"Number of GAIA sources": len(gaia_table)}, "Number of GAIA sources", descriptions={"Number of GAIA sources": 'Number of GAIA sources in image footprint'}, units={"Number of GAIA sources": "unitless"})  # write the number of GAIA sources
      diag_obj.write_json_file(hap_obj.drizzle_filename[:-9]+"_svm_gaia_sources.json", clobber=True)
  
      # Clean up
      del diag_obj
      del gaia_table
  
  # ----------------------------------------------------------------------------------------------------------------------
  
  
  def find_hap_point_sources(filt_obj, log_level=logutil.logging.NOTSET):
      """Identifies point sources in HAP imagery products and returns a dictionary containg **filt_obj** and
      a catalog of identified sources.
  
      Parameters
      ----------
      filt_obj : drizzlepac.haputils.Product.FilterProduct
          HAP filter product to process
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      Returns
      -------
      A three-element dictionary containing **filt_obj**, the source catalog and the temporary catalog filename
      keyed 'filt_obj', 'sources' and 'cat_name', respectively.
      """
      # Initiate logging!
      log.setLevel(log_level)
  
      # Need to create footprint mask of number_of_images_per_pixel
      footprint = cell_utils.SkyFootprint(filt_obj.meta_wcs)
      exp_list = [e.full_filename for e in filt_obj.edp_list]
      footprint.build(exp_list)
      # Initialize image object
      img_obj = catalog_utils.CatalogImage(filt_obj.drizzle_filename, footprint.total_mask, log_level)
      img_obj.compute_background(filt_obj.configobj_pars.get_pars("catalog generation")['bkg_box_size'],
                                 filt_obj.configobj_pars.get_pars("catalog generation")['bkg_filter_size'])
      img_obj.build_kernel(filt_obj.configobj_pars.get_pars("catalog generation")['bkg_box_size'],
                           filt_obj.configobj_pars.get_pars("catalog generation")['bkg_filter_size'],
                           filt_obj.configobj_pars.get_pars("catalog generation")['dao']['TWEAK_FWHMPSF'])
  
      # Perform background subtraction
      image = img_obj.data.copy()
      image -= img_obj.bkg_background_ra
  
      # create mask to reject any sources located less than 10 pixels from a image/chip edge
      wht_image = img_obj.data.copy()
      binary_inverted_wht = np.where(wht_image == 0, 1, 0)
      exclusion_mask = ndimage.binary_dilation(binary_inverted_wht, iterations=10)
  
      # Identify sources
      nsigma = filt_obj.configobj_pars.get_pars("alignment")["generate_source_catalogs"]["nsigma"]
      log.info("DAOStarFinder(fwhm={}, threshold={}*{})".format(img_obj.kernel_fwhm,
                                                                nsigma, img_obj.bkg_rms_median))
      daofind = DAOStarFinder(fwhm=img_obj.kernel_fwhm, threshold=nsigma * img_obj.bkg_rms_median)
      sources = daofind(image, mask=exclusion_mask)
      cat_name = filt_obj.product_basename + "_point-cat-fxm.ecsv"
  
      return {"filt_obj": filt_obj, "sources": sources, "cat_name": cat_name}
  
  
  # ----------------------------------------------------------------------------------------------------------------------
  
  
  def generate_gaia_catalog(hap_obj, columns_to_remove=None):
      """Uses astrometric_utils.create_astrometric_catalog() to create a catalog of all GAIA sources in the
      image footprint. This catalog contains right ascension, declination, and magnitude values, and is sorted
      in descending order by brightness.
  
      Parameters
      ----------
      hap_obj : drizzlepac.haputils.Product.TotalProduct, drizzlepac.haputils.Product.FilterProduct, or
          drizzlepac.haputils.Product.ExposureProduct, depending on input.
          hap product object to process
  
      columns_to_remove : list
          list of columns to remove from the table of GAIA sources returned by this function
  
      Returns
      -------
      gaia_table : astropy table
          table containing right ascension, declination, and magnitude of all GAIA sources identified in the
          image footprint, sorted in descending order by brightness.
      """
      # Gather list of input flc/flt images
      img_list = []
      log.debug("GAIA catalog will be created using the following input images:")
      # Create a list of the input flc.fits/flt.fits that were drizzled to create the final HAP product being
      # processed here. edp_item.info and hap_obj.info are both structured as follows:
      # <proposal id>_<visit #>_<instrument>_<detector>_<input filename>_<filter>_<drizzled product
      # image filetype>
      # Example: '10265_01_acs_wfc_j92c01b9q_flc.fits_f606w_drc'
      # what is being extracted here is just the input filename, which in this case is 'j92c01b9q_flc.fits'.
      if hasattr(hap_obj, "edp_list"):  # for total and filter product objects
          if not os.path.exists(hap_obj.drizzle_filename) or len(hap_obj.edp_list) == 0:
              # No valid products to evaluate
              log.warning('[generate_gaia_catalog] {} was not created.  Skipping.'.format(hap_obj.drizzle_filename))
              return None
  
          for edp_item in hap_obj.edp_list:
              parse_info = edp_item.info.split("_")
              imgname = "{}_{}".format(parse_info[4], parse_info[5])
              log.debug(imgname)
              img_list.append(imgname)
      else:  # For single-exposure product objects
          parse_info = hap_obj.info.split("_")
          imgname = "{}_{}".format(parse_info[4], parse_info[5])
          log.debug(imgname)
          img_list.append(imgname)
  
      # generate catalog of GAIA sources
      gaia_table = au.create_astrometric_catalog(img_list, gaia_only=True, use_footprint=True)
  
      if len(gaia_table) > 0:
          # trim off specified columns, but
          #    only if the specified columns already exist in the table
          #
          if columns_to_remove:
              existing_cols = [col for col in columns_to_remove if col in gaia_table.colnames]
              gaia_table.remove_columns(existing_cols)
  
          # remove sources outside image footprint
          outwcs = wcsutil.HSTWCS(hap_obj.drizzle_filename, ext=1)
          x, y = outwcs.all_world2pix(gaia_table['RA'], gaia_table['DEC'], 1)
          imghdu = fits.open(hap_obj.drizzle_filename)
          in_img_data = imghdu['WHT'].data.copy()
          in_img_data = np.where(in_img_data == 0, np.nan, in_img_data)
          mask = au.within_footprint(in_img_data, outwcs, x, y)
          gaia_table = gaia_table[mask]
  
      # Report results to log
      if len(gaia_table) == 0:
          log.warning("No GAIA sources were found!")
      elif len(gaia_table) == 1:
          log.info("1 GAIA source was found.")
      else:
          log.info("{} GAIA sources were found.".format(len(gaia_table)))
      return gaia_table
  
  # ----------------------------------------------------------------------------------------------------------------------
  
  
  def compare_photometry(drizzle_list, json_timestamp=None, json_time_since_epoch=None,
                         log_level=logutil.logging.NOTSET):
      """Compare photometry measurements for sources cross matched between the Point and Segment catalogs.
  
      DEPRECATED
  
      Report the magnitudes, as well as the mean difference, standard deviation of the mean, and median
      differences between the Point and Segment catalogs.
  
      Parameters
      ----------
      drizzle_list: list of strings
          Drizzle files for the Filter products which were mined to generate the output catalogs.
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      .. note:: This routine can be run either as a direct call from the hapsequencer.py routine,
      or it can invoked by a simple Python driver (or from within a Python session) by providing
      the names of the previously computed files as lists. The files must exist in the working directory.
      """
      log.setLevel(log_level)
      log.info('\n\n*****     Begin Quality Analysis Test: compare_photometry.     *****\n')
  
      pnt_suffix = '_point-cat.ecsv'
      seg_suffix = '_segment-cat.ecsv'
  
      good_flag_sum = 255
  
      phot_column_names = ["MagAp1", "MagAp2"]
      error_column_names = ["MagErrAp1", "MagErrAp2"]
  
      # Generate a separate JSON file for each detector and filter product
      # Drizzle filename example: hst_11665_06_wfc3_ir_f110w_ib4606_drz.fits.
      # The "product" in this context is a filter name.
      # The filename is all lower-case by design.
      for drizzle_file in drizzle_list:
          if not os.path.exists(drizzle_file):
              log.warning("[compare_photometry] Input {} not found. Skipping comparison.".format(drizzle_file))
              return  # So calling routine can continue to next test
  
          tokens = drizzle_file.split('_')
          detector = tokens[4]
          filter_name = tokens[5]
          ipppss = tokens[6]
  
          # Set up the diagnostic object
          diagnostic_obj = du.HapDiagnostic()
          diagnostic_obj.instantiate_from_fitsfile(drizzle_file,
                                                   data_source="{}.compare_photometry".format(__taskname__),
                                                   description="Photometry differences in Point and "
                                                               "Segment catalogs",
                                                   timestamp=json_timestamp,
                                                   time_since_epoch=json_time_since_epoch)
          summary_dict = {'detector': detector, 'filter_name': filter_name}
  
          # Construct the output JSON filename
          json_filename = '_'.join([ipppss, detector, 'svm', filter_name, 'photometry.json'])
  
          # Construct catalog names for catalogs that should have been produced
          # For any drizzled product, only two catalogs can be produced at most (point and segment).
          prefix = '_'.join(tokens[0:-1])
          cat_names = [prefix + pnt_suffix, prefix + seg_suffix]
  
          # Check that both catalogs exist
          for catalog in cat_names:
              does_exist = os.path.isfile(catalog)
              if not does_exist:
                  log.warning("Catalog {} does not exist.  Both the Point and Segment catalogs must exist "
                              "for comparison.".format(catalog))
                  log.warning("Program skipping comparison of catalogs associated "
                              "with {}.\n".format(drizzle_file))
                  return  # So calling routine can continue to next test
  
          # If the catalogs were actually produced, then get the data.
          tab_point_measurements = ascii.read(cat_names[0])
          tab_seg_measurements = ascii.read(cat_names[1])
  
          # Unfortunately the Point and Segment catalogs use different names for the X and Y values
          # Point: ([X|Y]-Center)  Segment: ([X|Y]-Centroid. Reset the coordinate columns to be only X or Y.
          tab_point_measurements.rename_column('X-Center', 'X')
          tab_point_measurements.rename_column('Y-Center', 'Y')
          tab_seg_measurements.rename_column('X-Centroid', 'X')
          tab_seg_measurements.rename_column('Y-Centroid', 'Y')
          cat_lengths = [len(tab_point_measurements), len(tab_seg_measurements)]
  
          # Determine the column names common to both catalogs as a list
          common_columns = list(set(tab_point_measurements.colnames).intersection(
              set(tab_seg_measurements.colnames)))
  
          # Use the utilities in devutils to match the sources in the two lists - get
          # the indices of the matches.
          matches_point_to_seg, matches_seg_to_point = cu.getMatchedLists(cat_names,
                                                                          [drizzle_file,
                                                                           drizzle_file],
                                                                          cat_lengths,
                                                                          log_level=log_level)
  
          # Move on to the next comparison without creating a .json if no cross-matches are found
          if len(matches_point_to_seg) == 0 or len(matches_seg_to_point) == 0:
              log.warning("Catalog {} and Catalog {} had no matching sources.".format(cat_names[0],
                                                                                      cat_names[1]))
              log.warning("Program skipping comparison of catalog indices associated "
                          "with {}. No JSON file will be produced.\n".format(drizzle_file))
              continue
  
          # There are nan values present in the catalogs - create a mask which identifies these rows
          # which are missing valid data
          missing_values_mask = cu.mask_missing_values(tab_point_measurements, tab_seg_measurements,
                                                       matches_point_to_seg, matches_seg_to_point,
                                                       common_columns)
  
          # Extract the Flag column from the two catalogs and get an ndarray (2, length)
          flag_matching = cu.extractMatchedLines('Flags', tab_point_measurements, tab_seg_measurements,
                                                 matches_point_to_seg, matches_seg_to_point)
  
          # Generate a mask to accommodate the missing, as well as the "flagged" entries
          flag_values_mask = cu.make_flag_mask(flag_matching, good_flag_sum, missing_values_mask)
  
          # Extract the columns of interest from the two catalogs for each desired measurement
          # and get an ndarray (2, length)
          # array([[21.512, ..., 2.944], [21.6 , ..., 22.98]],
          #       [[21.872, ..., 2.844], [21.2 , ..., 22.8]])
          for index, phot_column_name in enumerate(phot_column_names):
              matching_phot_rows = cu.extractMatchedLines(phot_column_name, tab_point_measurements,
                                                          tab_seg_measurements, matches_point_to_seg,
                                                          matches_seg_to_point, bitmask=flag_values_mask)
  
              # Compute the differences (Point - Segment)
              delta_phot = np.subtract(matching_phot_rows[0], matching_phot_rows[1])
  
              # Compute some basic statistics: mean difference and standard deviation, median difference,
              median_delta_phot = np.median(delta_phot)
              mean_delta_phot = np.mean(delta_phot)
              std_delta_phot = np.std(delta_phot)
  
              # NEED A BETTER WAY TO ASSOCIATE THE ERRORS WITH THE MEASUREMENTS
              # Compute the corresponding error of the differences
              matching_error_rows = cu.extractMatchedLines(error_column_names[index],
                                                           tab_point_measurements, tab_seg_measurements,
                                                           matches_point_to_seg, matches_seg_to_point,
                                                           bitmask=flag_values_mask)
  
              # Compute the error of the delta value (square root of the sum of the squares)
              result_error = np.sqrt(np.add(np.square(matching_error_rows[0]),
                                            np.square(matching_error_rows[1])))
  
              stat_key = 'Delta_' + phot_column_name
              stat_dict = {stat_key: {'Mean': mean_delta_phot, 'StdDev': std_delta_phot,
                           'Median': median_delta_phot}}
  
              # Write out the results
              diagnostic_obj.add_data_item(stat_dict,
                                           'Statistics_' + phot_column_name,
                                           descriptions={stat_key + '.Mean': phot_column_name + '_Mean_Differences(Point-Segment)',
                                                         stat_key + '.StdDev': phot_column_name + '_StdDev_of_Mean_Differences',
                                                         stat_key + '.Median': phot_column_name + '_Median_Differences(Point-Segment)'},
                                           units={stat_key + '.Mean': 'ABMag',
                                                  stat_key + '.StdDev': 'ABMag',
                                                  stat_key + '.Median': 'ABMag'})
  
          diagnostic_obj.write_json_file(json_filename)
          log.info("Generated photometry comparison for Point - Segment matches "
                   "sources {}.".format(json_filename))
  
          # Clean up
          del diagnostic_obj
  
  # ----------------------------------------------------------------------------------------------------------------------
  
  
  def report_wcs(total_product_list, json_timestamp=None, json_time_since_epoch=None,
                 log_level=logutil.logging.NOTSET):
      """Report the WCS information for each exposure of a total data product.
  
      Parameters
      ----------
      total_product_list: list of HAP TotalProduct objects, one object per instrument detector
      (drizzlepac.haputils.Product.TotalProduct)
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and
          written to the .log file.  Default value is 'NOTSET'.
      """
      log.setLevel(log_level)
      log.info('\n\n*****     Begin Quality Analysis Test: report_wcs.     *****\n')
      aposteriori_list = ['FIT', 'SVM', 'MVM']
  
      # Generate a separate JSON file for each ExposureProduct object
      for total_product in total_product_list:
          detector = total_product.detector
  
          # Loop over all the individual exposures in the list which comprise the total product
          for edp_object in total_product.edp_list:
              # Set up the diagnostic object
              diagnostic_obj = du.HapDiagnostic()
              diagnostic_obj.instantiate_from_hap_obj(edp_object,
                                                      data_source="{}.report_wcs".format(__taskname__),
                                                      description="WCS information",
                                                      timestamp=json_timestamp,
                                                      time_since_epoch=json_time_since_epoch)
  
              # Construct the output JSON filename
              json_filename = '_'.join([edp_object.product_basename, 'svm_wcs.json'])
  
              # For exposures with multiple science extensions (multiple chips),
              # generate a combined WCS
              num_sci_ext, extname = util.count_sci_extensions(edp_object.full_filename)
              extname_list = []
              for x in range(num_sci_ext):
                  extname_list.append((extname, x+1))
  
              metawcs = wcs_functions.make_mosaic_wcs(edp_object.full_filename)
  
              # Get information from the active WCS
              active_wcs_dict = {'primary_wcsname': metawcs.wcs.name,
                                 'crpix1': metawcs.wcs.crpix[0],
                                 'crpix2': metawcs.wcs.crpix[1],
                                 'crval1': metawcs.wcs.crval[0],
                                 'crval2': metawcs.wcs.crval[1],
                                 'scale': metawcs.pscale,
                                 'orientation': metawcs.orientat,
                                 'exposure': edp_object.exposure_name}
  
              diagnostic_obj.add_data_item(active_wcs_dict, 'PrimaryWCS',
                                           descriptions={'primary_wcsname': 'Active WCS',
                                                         'crpix1': 'X coord of reference pixel',
                                                         'crpix2': 'Y coord of reference pixel',
                                                         'crval1': 'RA of reference pixel',
                                                         'crval2': 'Dec of reference pixel',
                                                         'scale': 'Plate scale',
                                                         'orientation': 'Position angle of Image Y axis (East of North)',
                                                         'exposure': 'Exposure name'},
                                           units={'primary_wcsname': 'unitless',
                                                  'crpix1': 'pixels',
                                                  'crpix2': 'pixels',
                                                  'crval1': 'degrees',
                                                  'crval2': 'degrees',
                                                  'scale': 'arcseconds/pixel',
                                                  'orientation': 'degrees',
                                                  'exposure': 'unitless'})
  
              # Determine the possible alternate WCS solutions in the header
              dict_of_wcskeys_names = wcsutil.altwcs.wcsnames(edp_object.full_filename, ext=1)
  
              # Ignore the OPUS ("O") WCS, as well as the duplicate of the active WCS
              dict_of_wcskeys_names.pop('O')
  
              reverse_dict = {}
              for key, value in dict_of_wcskeys_names.items():
                  reverse_dict.setdefault(value, set()).add(key)
              keys_with_dups = set(chain.from_iterable(values for key, values in reverse_dict.items() if len(values) > 1))
  
              # Make a list of the keys which contain duplicate values...
              if keys_with_dups:
                  list_keys = list(keys_with_dups)
                  # ...ignore the primary key as it is important, and...
                  #MDD
                  #list_keys.remove(' ')
                  # ...remove the duplicates.
                  for popkey in list_keys:
                      dict_of_wcskeys_names.pop(popkey)
  
              # The remaining dictionary items all need to share the same IDC base
              # solution as the "active" solution in order to make a consistent comparison -
              # remove any outliers.
              loc = metawcs.wcs.name.find("_")
              root_idc = metawcs.wcs.name[0:loc]
  
              bad_match_key = []
              for key, value in dict_of_wcskeys_names.items():
                  if root_idc in value:
                      continue
                  else:
                      bad_match_key.append(key)
  
              for bad_key in bad_match_key:
                  dict_of_wcskeys_names.pop(bad_key)
  
              log.info("Remaining WCS keys and names based on same reference root {}".format(dict_of_wcskeys_names))
  
              # If there is anything left to compare, then do it.
              # Want to gather the WCS for the default (IDC_rootname), apriori (any name without
              # FIT, SVM, or MVM), and aposteriori (any name containing FIT, SVM, or MVM) solutions.
              if len(dict_of_wcskeys_names) > 1:
  
                  # Activate an alternate WCS in order to gather its information.
                  # First copy the original primary WCS to an alternate (in case there was
                  # not already a duplicate). Use key 'Z'.  *** FIX MDD Should check for Z in use.
                  wcsutil.altwcs.archiveWCS(edp_object.full_filename, ext=extname_list, wcskey='Z')
  
                  icnt = 0
                  # Restore an alternate to be the primary WCS
                  # Examples of WCS name formats:
                  # Default WCS name: IDC_xxxxxxxxx
                  # A priori WCS name: IDC_xxxxxxxxx-HSC30
                  # A posteriori WCS name: IDC_xxxxxxxxx-FIT_REL_GAIADR1
                  for key, value in dict_of_wcskeys_names.items():
                      if key != ' ':
                          alt_key = key
                          alt_wcs_name = value.strip()
  
                          if '-' not in alt_wcs_name and alt_wcs_name.startswith('IDC_'):
                              alt_name = 'AlternateWCS_default'
                              delta_name = 'DeltaWCS_default'
                          elif '-' in alt_wcs_name and 'FIT' not in alt_wcs_name:
                              alt_name = 'AlternateWCS_apriori'
                              delta_name = 'DeltaWCS_apriori'
                          elif list(filter(alt_wcs_name.find, aposteriori_list)) != []:
                              alt_name = 'AlternateWCS_aposteriori'
                              delta_name = 'DeltaWCS_aposteriori'
                          else:
                              log.info('WCS name {} is unknown.  Moving on to next WCS solution.\n'.format(alt_wcs_name))
                              continue
  
                          wcsutil.altwcs.restoreWCS(edp_object.full_filename, ext=extname_list, wcskey=key)
  
                          # Create a metawcs for this alternate WCS
                          alt_metawcs = wcs_functions.make_mosaic_wcs(edp_object.full_filename)
  
                          # Get information from the alternate/active WCS
                          alt_wcs_dict = {'alt_wcsname': alt_wcs_name,
                                          'crpix1': alt_metawcs.wcs.crpix[0],
                                          'crpix2': alt_metawcs.wcs.crpix[1],
                                          'crval1': alt_metawcs.wcs.crval[0],
                                          'crval2': alt_metawcs.wcs.crval[1],
                                          'scale': alt_metawcs.pscale,
                                          'orientation': alt_metawcs.orientat,
                                          'exposure': edp_object.exposure_name}
  
                          diagnostic_obj.add_data_item(alt_wcs_dict, alt_name,
                                                       descriptions={'alt_wcsname': 'Alternate WCS',
                                                                     'crpix1': 'X coord of reference pixel',
                                                                     'crpix2': 'Y coord of reference pixel',
                                                                     'crval1': 'RA of reference pixel',
                                                                     'crval2': 'Dec of reference pixel',
                                                                     'scale': 'Plate scale',
                                                                     'orientation': 'Position angle of Image Y axis (East of North)',
                                                                     'exposure': 'Exposure name'},
                                                       units={'alt_wcsname': 'unitless',
                                                              'crpix1': 'pixels',
                                                              'crpix2': 'pixels',
                                                              'crval1': 'degrees',
                                                              'crval2': 'degrees',
                                                              'scale': 'arcseconds/pixel',
                                                              'orientation': 'degrees',
                                                              'exposure': 'unitless'})
  
                          delta_wcs_name = metawcs.wcs.name + '_minus_' + alt_wcs_name
                          diff_wcs_dict = {'delta_wcsname': delta_wcs_name,
                                           'd_crpix1': active_wcs_dict['crpix1'] - alt_wcs_dict['crpix1'],
                                           'd_crpix2': active_wcs_dict['crpix2'] - alt_wcs_dict['crpix2'],
                                           'd_crval1': active_wcs_dict['crval1'] - alt_wcs_dict['crval1'],
                                           'd_crval2': active_wcs_dict['crval2'] - alt_wcs_dict['crval2'],
                                           'd_scale': active_wcs_dict['scale'] - alt_wcs_dict['scale'],
                                           'd_orientation': active_wcs_dict['orientation'] - alt_wcs_dict['orientation'],
                                           'exposure': edp_object.exposure_name}
  
                          diagnostic_obj.add_data_item(diff_wcs_dict, delta_name,
                                                       descriptions={'delta_wcsname': 'Active-Alternate WCS',
                                                                     'd_crpix1': 'delta_X of reference pixel',
                                                                     'd_crpix2': 'delta_Y of reference pixel',
                                                                     'd_crval1': 'delta_RA of reference pixel',
                                                                     'd_crval2': 'delta_Dec of reference pixel',
                                                                     'd_scale': 'delta_Plate scale',
                                                                     'd_orientation': 'delta_Position angle of Image Y axis (East of North)',
                                                                     'exposure': 'Exposure name'},
                                                       units={'delta_wcsname': 'unitless',
                                                              'd_crpix1': 'pixels',
                                                              'd_crpix2': 'pixels',
                                                              'd_crval1': 'degrees',
                                                              'd_crval2': 'degrees',
                                                              'd_scale': 'arcseconds/pixel',
                                                              'd_orientation': 'degrees',
                                                              'exposure': 'unitless'})
  
                          # Delete the original alternate WCS...
                          wcsutil.altwcs.deleteWCS(edp_object.full_filename, ext=extname_list, wcskey=alt_key)
  
                          # ... and archive the current primary with its original key
                          wcsutil.altwcs.archiveWCS(edp_object.full_filename, ext=extname_list, wcskey=alt_key)
  
                          icnt += 1
  
                      else:
                          continue
  
                  # When comparisons are done between the original primary WCS and all
                  # the alternates, restore the original primary WCS with key Z.
                  wcsutil.altwcs.restoreWCS(edp_object.full_filename, ext=extname_list, wcskey='Z')
  
                  # Delete the extra copy of the primary
                  wcsutil.altwcs.deleteWCS(edp_object.full_filename, ext=extname_list, wcskey='Z')
  
              else:
                  log.info("This dataset only has the Primary and OPUS WCS values")
  
              diagnostic_obj.write_json_file(json_filename)
  
              # Clean up
              del diagnostic_obj
  
      # This routine does not return any values
  
  
  # ----------------------------------------------------------------------------------------------------------------------
  
  
  def run_hla_sourcelist_comparison(total_list, diagnostic_mode=False, json_timestamp=None,
                                    json_time_since_epoch=None, log_level=logutil.logging.NOTSET):
      """ This subroutine automates execution of drizzlepac/devutils/comparison_tools/compare_sourcelists.py to
      compare HAP-generated filter catalogs with their HLA classic counterparts.
  
      NOTE: In order for this subroutine to run, the following environment variables need to be set:
      - HLA_CLASSIC_BASEPATH
      - HLA_BUILD_VER
  
      Alternatively, if the HLA classic path is unavailable, The comparison can be run using locally stored HLA
      classic files. The relevant HLA classic imagery and sourcelist files must be placed in a subdirectory of
      the current working directory called 'hla_classic'.
  
      Parameters
      ----------
      total_list: list
          List of TotalProduct objects, one object per instrument/detector combination is
          a visit.  The TotalProduct objects are comprised of FilterProduct and ExposureProduct
          objects.
  
      diagnostic_mode : Boolean, optional.
          create intermediate diagnostic files? Default value is False.
  
      json_timestamp: str, optional
          Universal .json file generation date and time (local timezone) that will be used in the instantiation
          of the HapDiagnostic object. Format: MM/DD/YYYYTHH:MM:SS (Example: 05/04/2020T13:46:35). If not
          specified, default value is logical 'None'
  
      json_time_since_epoch : float
          Universal .json file generation time that will be used in the instantiation of the HapDiagnostic
          object. Format: Time (in seconds) elapsed since January 1, 1970, 00:00:00 (UTC). If not specified,
          default value is logical 'None'
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file.  Default value is 'NOTSET'.
  
      RETURNS
      -------
      Nothing.
      """
      log.setLevel(log_level)
      log.info('\n\n*****     Begin Quality Analysis Test: run_hla_sourcelist_comparison.     *****\n')
      # get HLA classic path details from envroment variables
      hla_classic_basepath = os.getenv('HLA_CLASSIC_BASEPATH')
      hla_build_ver = os.getenv("HLA_BUILD_VER")
      for tot_obj in total_list:
          combo_comp_pdf_list = []
          if hla_classic_basepath and hla_build_ver and os.path.exists(hla_classic_basepath):
              hla_cassic_basepath = os.path.join(hla_classic_basepath, tot_obj.instrument, hla_build_ver)
              hla_classic_path = os.path.join(hla_cassic_basepath,
                                              tot_obj.prop_id,
                                              tot_obj.prop_id + "_" + tot_obj.obset_id)  # Generate path to HLA classic products
          elif os.path.exists(os.path.join(os.getcwd(), "hla_classic")):  # For local testing
              hla_classic_basepath = os.path.join(os.getcwd(), "hla_classic")
              hla_classic_path = hla_classic_basepath
          else:
              return  # bail out if HLA classic path can't be found.
          for filt_obj in tot_obj.fdp_list:
              hap_imgname = filt_obj.drizzle_filename
              hla_imgname = glob.glob("{}/{}{}_dr*.fits".format(hla_classic_path,
                                                                filt_obj.basename,
                                                                filt_obj.filters))[0]
              if not os.path.exists(hap_imgname) or not os.path.exists(hla_imgname):  # Skip filter if one or both of the images can't be found
                  continue
              for hap_sourcelist_name in [filt_obj.point_cat_filename, filt_obj.segment_cat_filename]:
                  if hap_sourcelist_name.endswith("point-cat.ecsv"):
                      hla_classic_cat_type = "dao"
                      plotfile_prefix = filt_obj.product_basename + "_point"
                  else:
                      hla_classic_cat_type = "sex"
                      plotfile_prefix = filt_obj.product_basename + "_segment"
                  if hla_classic_basepath and hla_build_ver and os.path.exists(hla_classic_basepath):
                      hla_sourcelist_name = "{}/logs/{}{}_{}phot.txt".format(hla_classic_path,
                                                                             filt_obj.basename,
                                                                             filt_obj.filters,
                                                                             hla_classic_cat_type)
                  else:
                      hla_sourcelist_name = "{}/{}{}_{}phot.txt".format(hla_classic_path, filt_obj.basename,
                                                                        filt_obj.filters, hla_classic_cat_type)
                  if not os.path.exists(hap_sourcelist_name) or not os.path.exists(hla_sourcelist_name):  # Skip catalog type if one or both of the catalogs can't be found
                      continue
  
                  # convert HLA Classic RA and Dec values to HAP reference frame so the RA and Dec comparisons are correct
                  updated_hla_sourcelist_name = correct_hla_classic_ra_dec(hla_sourcelist_name,
                                                                           hap_imgname,
                                                                           hla_classic_cat_type,
                                                                           log_level)
                  log.info("HAP image:                   {}".format(os.path.basename(hap_imgname)))
                  log.info("HLA Classic image:           {}".format(os.path.basename(hla_imgname)))
                  log.info("HAP catalog:                 {}".format(os.path.basename(hap_sourcelist_name)))
                  log.info("HLA Classic catalog:         {}".format(os.path.basename(updated_hla_sourcelist_name)))
                  # once all file exist checks are passed, execute sourcelist comparison
                  return_status = csl.comparesourcelists(slNames=[updated_hla_sourcelist_name,
                                                                  hap_sourcelist_name],
                                                         imgNames=[hla_imgname, hap_imgname],
                                                         good_flag_sum=255,
                                                         plotGen="file",
                                                         plotfile_prefix=plotfile_prefix,
                                                         output_json_filename=hap_sourcelist_name.replace(".ecsv",
                                                                                                          "_svm_compare_sourcelists.json"),
                                                         verbose=True,
                                                         log_level=log_level,
                                                         json_timestamp=json_timestamp,
                                                         json_time_since_epoch=json_time_since_epoch,
                                                         debugMode=diagnostic_mode)
                  combo_comp_pdf_filename = "{}_comparison_plots.pdf".format(plotfile_prefix)
                  if os.path.exists(combo_comp_pdf_filename):
                      combo_comp_pdf_list.append(combo_comp_pdf_filename)
          if len(combo_comp_pdf_list) > 0:  # combine all plots generated by compare_sourcelists.py for this total object into a single pdf file
              total_combo_comp_pdf_filename = "{}_svm_comparison_plots.pdf".format(tot_obj.drizzle_filename[:-9].replace("_total", ""))
              cu.pdf_merger(total_combo_comp_pdf_filename, combo_comp_pdf_list)
              log.info("Sourcelist comparison plots saved to file {}.".format(total_combo_comp_pdf_filename))
  
  
  # ----------------------------------------------------------------------------------------------------------------------
  
  def correct_hla_classic_ra_dec(orig_hla_classic_sl_name, hap_imgname, cattype, log_level):
      """
      This subroutine runs Rick White's read_hla_catalog script to convert the RA and Dec values from a HLA Classic
      sourcelist into same reference frame used by the HAP sourcelists. Additionally, the new RA and Dec values
      are then transformed to produce X and Y coordinates that are in the HAP image frame of reference. A new
      version of the input file with the converted X and Y and RA and Dec values is written to the current
      working directory named <INPUT SOURCELIST NAME>_corrected.txt.
  
      Parameters
      ----------
      orig_hla_classic_sl_name : string
          name of the HLA Classic sourcelist whose RA and Dec values will be converted.
  
      hap_imgname : string
          name of HAP image. The WCS info from this image will be used to transform the updated RA and DEC
          values to new X and Y values in the this image's frame of reference
  
      cattype : string
          HLA Classic catalog type. Either 'sex' (source extractor) or 'dao' (DAOphot).
  
      log_level : int
          The desired level of verboseness in the log statements displayed on the screen and written to the .log file.
  
      Returns
      -------
      mod_sl_name : string
          Name of the new version of the input file with the converted RA and Dec values
      """
      try:
          mod_sl_name = os.path.basename(orig_hla_classic_sl_name)
  
          # Execute read_hla_catalog.read_hla_catalog() to convert RA and Dec values
          dataset = mod_sl_name.replace("_{}phot.txt".format(cattype), "")
          modcat = read_hla_catalog.read_hla_catalog(dataset, cattype=cattype, applyomega=True, multiwave=False,
                                                     verbose=True, trim=False, log_level=log_level)
          # sort catalog with updated RA, DEC values so that ordering is the same as the uncorrected table and everything maps correclty.
          if cattype == "dao":
              sortcoltitle = "ID"
              x_coltitle = "X-Center"
              y_coltitle = "Y-Center"
          if cattype == "sex":
              sortcoltitle = "NUMBER"
              x_coltitle = "X-IMAGE"
              y_coltitle = "Y-IMAGE"
  
          modcat.sort(sortcoltitle)
  
          # Identify RA and Dec column names in the new catalog table object
          for ra_col_title in ["ra", "RA", "ALPHA_J2000", "alpha_j2000"]:
              if ra_col_title in modcat.colnames:
                  true_ra_col_title = ra_col_title
                  log.debug("RA Col_name: {}".format(true_ra_col_title))
                  break
          for dec_col_title in ["dec", "DEC", "Dec", "DELTA_J2000", "delta_j2000"]:
              if dec_col_title in modcat.colnames:
                  true_dec_col_title = dec_col_title
                  log.debug("DEC Col_name: {}".format(true_dec_col_title))
                  break
  
          # transform new RA and Dec values into X and Y values in the HAP reference frame
          ra_dec_values = np.stack((modcat[true_ra_col_title], modcat[true_dec_col_title]), axis=1)
          new_xy_values = hla_flag_filter.rdtoxy(ra_dec_values, hap_imgname, '[sci,1]')
  
          # get HLA Classic sourcelist data, replace existing RA and Dec column data with the converted RA and Dec column data
          cat = Table.read(orig_hla_classic_sl_name, format='ascii')
          cat['RA'] = modcat[true_ra_col_title]
          cat['DEC'] = modcat[true_dec_col_title]
  
          # update existing X and Y values with new X and Y values transformed from new RA and Dec values.
          cat[x_coltitle] = new_xy_values[:, 0]
          cat[y_coltitle] = new_xy_values[:, 1]
  
          # Write updated version of HLA Classic catalog to current working directory
          mod_sl_name = mod_sl_name.replace(".txt", "_corrected.txt")
          log.info("Corrected version of HLA Classic file {} with new X, Y and RA, Dec values written to {}.".format(os.path.basename(orig_hla_classic_sl_name), mod_sl_name))
          cat.write(mod_sl_name, format="ascii.csv")
          return mod_sl_name
  
      except Exception:
          log.warning("There was a problem converting the RA and Dec values. Using original uncorrected HLA Classic sourcelist instead.")
          log.warning("Comparisons may be of questionable quality.")
          return orig_hla_classic_sl_name
  
  
  # ----------------------------------------------------------------------------------------------------------------------
  def run_quality_analysis(total_obj_list, run_compare_num_sources=True, run_find_gaia_sources=True,
                           run_compare_hla_sourcelists=True, run_compare_ra_dec_crossmatches=True,
                           run_characterize_gaia_distribution=True, run_compare_photometry=True,
                           run_compare_interfilter_crossmatches=True, run_report_wcs=True,
                           log_level=logutil.logging.NOTSET):
      """Run the quality analysis functions
  
      Parameters
      ----------
      total_obj_list : list
          List of one or more HAP drizzlepac.haputils.Product.TotalProduct object(s) to process
  
      run_compare_num_sources : bool, optional
          Run 'compare_num_sources' test? Default value is True.
  
      run_find_gaia_sources : bool, optional
          Run 'find_gaia_sources' test? Default value is True.
  
      run_compare_hla_sourcelists : bool, optional
          Run 'run_compare_sourcelists' test? Default value is True.
  
      run_compare_ra_dec_crossmatches : bool, optional
          Run 'compare_ra_dec_crossmatches' test? Default value is True.
  
      run_characterize_gaia_distribution : bool, optional
          Run 'characterize_gaia_distribution' test? Default value is True.
  
      run_compare_photometry : bool, optional
          Run 'compare_photometry' test? Default value is True.
  
      run_compare_interfilter_crossmatches : bool, optional
          Run 'compare_filter_cross_match' test? Default value is True.
  
      run_report_wcs : bool, optional
          Run 'report_wcs' test? Devault value is True.
  
      log_level : int, optional
          The desired level of verboseness in the log statements displayed on the screen and written to the
          .log file. Default value is 'NOTSET'.
  
      Returns
      -------
      Nothing.
      """
      log.setLevel(log_level)
  
      # generate a timestamp values that will be used to make creation time, creation date and epoch values
      # common to each json file
      json_timestamp = datetime.now().strftime("%m/%d/%YT%H:%M:%S")
      json_time_since_epoch = time.time()
  
      # Determine number of sources in Point and Segment catalogs
      if run_compare_num_sources:
          try:
              total_drizzle_list = []
              for total_obj in total_obj_list:
                  total_drizzle_list.append(total_obj.drizzle_filename)
              compare_num_sources(total_drizzle_list, json_timestamp=json_timestamp,
                                  json_time_since_epoch=json_time_since_epoch, log_level=log_level)
          except Exception:
              log.warning("HAP Point vs. HAP Segment sourcelist length comparison (compare_num_sources) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
      # Identify the number of GAIA sources in final product footprints
      if run_find_gaia_sources:
          try:
              for total_obj in total_obj_list:
                  find_gaia_sources(total_obj, json_timestamp=json_timestamp,
                                    json_time_since_epoch=json_time_since_epoch, log_level=log_level)
                  for filter_obj in total_obj.fdp_list:
                      find_gaia_sources(filter_obj, json_timestamp=json_timestamp,
                                        json_time_since_epoch=json_time_since_epoch, log_level=log_level)
                      for exp_obj in filter_obj.edp_list:
                          find_gaia_sources(exp_obj, json_timestamp=json_timestamp,
                                            json_time_since_epoch=json_time_since_epoch, log_level=log_level)
          except Exception:
              log.warning("GAIA sources count (find_gaia_sources) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
      # Compare HAP sourcelists to their HLA Classic counterparts
      if run_compare_hla_sourcelists:
          try:
              if log_level == logutil.logging.DEBUG:
                  diag_mode = True
              else:
                  diag_mode = False
              run_hla_sourcelist_comparison(total_obj_list,
                                            diagnostic_mode=diag_mode,
                                            json_timestamp=json_timestamp,
                                            json_time_since_epoch=json_time_since_epoch,
                                            log_level=log_level)
          except Exception:
              log.warning("HAP vs. HLA sourcelist comparison (compare_sourcelists) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
      # Get point/segment cross-match RA/Dec statistics
      if run_compare_ra_dec_crossmatches:
          try:
              for total_obj in total_obj_list:
                  for filter_obj in total_obj.fdp_list:
                      compare_ra_dec_crossmatches(filter_obj, json_timestamp=json_timestamp,
                                                  json_time_since_epoch=json_time_since_epoch, log_level=log_level)
          except Exception:
              log.warning("HAP Point vs. HAP Segment sourcelist RA/Dec comparison (compare_ra_dec_crossmatches) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
      # Statistically characterize GAIA distribution
      if run_characterize_gaia_distribution:
          try:
              for total_obj in total_obj_list:
                  for filter_obj in total_obj.fdp_list:
                      characterize_gaia_distribution(filter_obj, json_timestamp=json_timestamp,
                                                     json_time_since_epoch=json_time_since_epoch,
                                                     log_level=log_level)
          except Exception:
              log.warning("GAIA source distribution characterization (characterize_gaia_distribution) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
      # Photometry of cross-matched sources in Point and Segment catalogs for Filter products
      if run_compare_photometry:
          try:
              tot_len = len(total_obj_list)
              filter_drizzle_list = []
              temp_list = []
              for tot in total_obj_list:
                  temp_list = [x.drizzle_filename for x in tot.fdp_list]
                  filter_drizzle_list.extend(temp_list)
              compare_photometry(filter_drizzle_list, json_timestamp=json_timestamp,
                                 json_time_since_epoch=json_time_since_epoch, log_level=log_level)
          except Exception:
              log.warning("HAP Point vs. HAP Segment sourcelist photometry comparison (compare_photometry) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
      # Compare inter-filter cross matched HAP sources
      if run_compare_interfilter_crossmatches:
          try:
              compare_interfilter_crossmatches(total_obj_list, json_timestamp=json_timestamp,
                                               json_time_since_epoch=json_time_since_epoch, log_level=log_level)
          except Exception:
              log.warning("HAP Point sourcelist interfilter comparison (compare_interfilter_crossmatches) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
      # Report WCS info
      if run_report_wcs:
          try:
              report_wcs(total_obj_list, json_timestamp=json_timestamp, json_time_since_epoch=json_time_since_epoch,
                         log_level=log_level)
          except Exception:
              log.warning("WCS reporting (report_wcs) encountered a problem.")
              log.exception("message")
              log.warning("Continuing to next test...")
  
  
  # ============================================================================================================
  
  
  if __name__ == "__main__":
      # process command-line inputs with argparse
      parser = argparse.ArgumentParser(description='Perform quality assessments of the SVM products generated '
                                                   'by the drizzlepac package. NOTE: if no QA switches (-cgd, '
                                                   '-cns, -cp, -cxm, or -fgs) are specified, ALL QA steps will '
                                                   'be executed.')
      parser.add_argument('input_filename', help='_total_list.pickle file that holds vital information about '
                                                 'the processing run')
      parser.add_argument('-cgd', '--run_characterize_gaia_distribution', required=False, action='store_true',
                          help="Statistically describe distribution of GAIA sources in footprint.")
      parser.add_argument('-cns', '--run_compare_num_sources', required=False, action='store_true',
                          help='Determine the number of viable sources actually listed in SVM output catalogs.')
      parser.add_argument('-cp', '--run_compare_photometry', required=False, action='store_true',
                          help="Compare photometry measurements for sources cross matched between the Point "
                               "and Segment catalogs.")
      parser.add_argument('-cxm', '--run_compare_ra_dec_crossmatches', required=False, action='store_true',
                          help="Compare RA/Dec position measurements for sources cross matched between the "
                               "Point and Segment catalogs.")
      parser.add_argument('-fgs', '--run_find_gaia_sources', required=False, action='store_true',
                          help="Determine the number of GAIA sources in the footprint of a specified HAP final "
                               "product image")
      parser.add_argument('-fxm', '--run_compare_interfilter_crossmatches', required=False, action='store_true',
                          help="Compare positions of cross matched sources from different filter-level "
                               "products")
      parser.add_argument('-hla', '--run_compare_hla_sourcelists', required=False, action='store_true',
                          help="Compare HAP sourcelists to their HLA classic counterparts")
      parser.add_argument('-wcs', '--run_report_wcs', required=False, action='store_true',
                          help="Report the WCS information for each exposure of a total data product")
      parser.add_argument('-l', '--log_level', required=False, default='info',
                          choices=['critical', 'error', 'warning', 'info', 'debug'],
                          help='The desired level of verboseness in the log statements displayed on the screen '
                               'and written to the .log file. The level of verboseness from left to right, and '
                               'includes all log statements with a log_level left of the specified level. '
                               'Specifying "critical" will only record/display "critical" log statements, and '
                               'specifying "error" will record/display both "error" and "critical" log '
                               'statements, and so on.')
      user_args = parser.parse_args()
  
      # set up logging
      log_dict = {"critical": logutil.logging.CRITICAL,
                  "error": logutil.logging.ERROR,
                  "warning": logutil.logging.WARNING,
                  "info": logutil.logging.INFO,
                  "debug": logutil.logging.DEBUG}
      log_level = log_dict[user_args.log_level]
      log.setLevel(log_level)
  
      # verify that input file exists
      if not os.path.exists(user_args.input_filename):
          err_msg = "File {} doesn't exist.".format(user_args.input_filename)
          log.critical(err_msg)
          raise Exception(err_msg)
  
      #  check that at least one QA switch is turned on
      all_qa_steps_off = True
      max_step_str_length = 0
      for kv_pair in user_args._get_kwargs():
          if kv_pair[0] not in ['input_filename', 'run_all', 'log_level']:
              if len(kv_pair[0])-4 > max_step_str_length:
                  max_step_str_length = len(kv_pair[0])-4
              if kv_pair[1]:
                  all_qa_steps_off = False
  
      # if no QA steps are explicitly turned on in the command-line call, run ALL the QA steps
      if all_qa_steps_off:
          log.info("No specific QA switches were turned on. All QA steps will be executed.")
          user_args.run_characterize_gaia_distribution = True
          user_args.run_compare_num_sources = True
          user_args.run_compare_photometry = True
          user_args.run_compare_hla_sourcelists = True
          user_args.run_compare_ra_dec_crossmatches = True
          user_args.run_compare_interfilter_crossmatches = True
          user_args.run_find_gaia_sources = True
          user_args.run_report_wcs = True
  
      # display status summary indicating which QA steps are turned on and which steps are turned off
      toplinestring = "-"*(int(max_step_str_length/2)-6)
      log.info("{}QA step run status{}".format(toplinestring, toplinestring))
      for kv_pair in user_args._get_kwargs():
          if kv_pair[0] not in ['input_filename', 'run_all', 'log_level']:
              if kv_pair[1]:
                  run_status = "ON"
              else:
                  run_status = "off"
              log.info("{}{}   {}".format(kv_pair[0][4:], " "*(max_step_str_length-(len(kv_pair[0])-4)),
                                          run_status))
      log.info("-"*(max_step_str_length+6))
  
      # execute specified tests
      filehandler = open(user_args.input_filename, 'rb')
      total_obj_list = pickle.load(filehandler)
      run_quality_analysis(total_obj_list,
                           run_compare_num_sources=user_args.run_compare_num_sources,
                           run_find_gaia_sources=user_args.run_find_gaia_sources,
                           run_compare_ra_dec_crossmatches=user_args.run_compare_ra_dec_crossmatches,
                           run_characterize_gaia_distribution=user_args.run_characterize_gaia_distribution,
                           run_compare_hla_sourcelists=user_args.run_compare_hla_sourcelists,
                           run_compare_photometry=user_args.run_compare_photometry,
                           run_compare_interfilter_crossmatches=user_args.run_compare_interfilter_crossmatches,
                           run_report_wcs=user_args.run_report_wcs,
                           log_level=log_level)
  
  
  # -----------------------------------------------------------------------------
  #
  # Generate plots for these results
  #
  # -----------------------------------------------------------------------------
  HOVER_COLUMNS = ['gen_info.instrument',
                   'gen_info.detector',
                   'gen_info.filter',
                   'gen_info.imgname',
                   'header.DATE-OBS',
                   'header.RA_TARG',
                   'header.DEC_TARG',
                   'header.GYROMODE',
                   'header.EXPTIME']
  
  FIGURE_TOOLS = 'pan,wheel_zoom,box_zoom,zoom_in,zoom_out,xbox_select,reset,save'
  
  
  def build_svm_plots(data_source, output_basename=''):
      """Create all the plots for the results generated by these comparisons
  
      Parameters
      ----------
      data_source : str
          Filename for master data file which contains all the results.  This will
          typically be an HSF5 file generated by the JSON harvester.
  
      """
      if output_basename == '':
          output_basename = "svm_qa"
      else:
          output_basename = "{}_svm_qa".format(output_basename)
  
      # Generate plots for point-segment catalog cross-match comparisons
      xmatch_col_names = HOVER_COLUMNS + ['Cross-match_details.number_of_cross-matches',
                                          'Cross-match_details.point_catalog_filename',
                                          'Cross-match_details.point_catalog_length',
                                          'Cross-match_details.point_frame',
                                          'Cross-match_details.segment_catalog_filename',
                                          'Cross-match_details.segment_catalog_length',
                                          'Cross-match_details.segment_frame',
                                          'Cross-matched_point_catalog.Declination',
                                          'Cross-matched_point_catalog.Right ascension',
                                          'Cross-matched_segment_catalog.Declination',
                                          'Cross-matched_segment_catalog.Right ascension',
                                          'Segment_-_point_on-sky_separation_statistics.3x3_sigma-clipped_mean',
                                          'Segment_-_point_on-sky_separation_statistics.3x3_sigma-clipped_median',
                                          'Segment_-_point_on-sky_separation_statistics.3x3_sigma-clipped_standard_deviation',
                                          'Segment_-_point_on-sky_separation_statistics.Non-clipped_max',
                                          'Segment_-_point_on-sky_separation_statistics.Non-clipped_mean',
                                          'Segment_-_point_on-sky_separation_statistics.Non-clipped_median',
                                          'Segment_-_point_on-sky_separation_statistics.Non-clipped_min',
                                          'Segment_-_point_on-sky_separation_statistics.Non-clipped_standard_deviation']
  
      xmatch_cols = get_pandas_data(data_source, xmatch_col_names)
  
      xmatch_plots_name = build_crossmatch_plots(xmatch_cols, xmatch_col_names, output_basename=output_basename)
  
  
  # -----------------------------------------------------------------------------
  # Functions for generating each data plot
  #
  
  def build_crossmatch_plots(xmatchCDS, data_cols, output_basename='svm_qa'):
      """
      Generate the cross-match statistics plots for the comparison between the
      point catalog and the segment catalog.
  
      Parameters
      ----------
      xmatchCDS : Pandas ColumnDataSource
          This object contains all the columns relevant to the cross-match plots.
  
      data_cols : list
          The list of column names for the columns read in to the `xmatchCDS` object.
  
      output_basename : str
          String to use as the start of the filename for the output plot pages.
  
      Returns
      -------
      output : str
          Name of HTML file where the plot was saved.
  
      """
      output_basename = "{}_crossmatch_comparison".format(output_basename)
  
      if not output_basename.endswith('.html'):
          output = output_basename + '.html'
      else:
          output = output_basename
      # Set the output file immediately as advised by Bokeh.
      output_file(output)
  
      num_hover_cols = len(HOVER_COLUMNS)
  
      colormap = [qa.DETECTOR_LEGEND[x] for x in xmatchCDS.data[data_cols[1]]]
      xmatchCDS.data['colormap'] = colormap
      inst_det = ["{}/{}".format(i, d) for (i, d) in zip(xmatchCDS.data[data_cols[0]],
                                                         xmatchCDS.data[data_cols[1]])]
      xmatchCDS.data[qa.INSTRUMENT_COLUMN] = inst_det
  
      plot_list = []
  
      hist0, edges0 = np.histogram(xmatchCDS.data[data_cols[num_hover_cols]], bins=50)
      title0 = 'Number of Point-to-Segment Cross-matched sources'
      p0 = [plot_histogram(title0, hist0, edges0, y_start=0, fill_color='navy',
                           background_fill_color='#fafafa', xlabel='Number of Cross-matched sources',
                           ylabel='Number of products')]
      plot_list += p0
  
      hist1, edges1 = np.histogram(xmatchCDS.data[data_cols[num_hover_cols + 11]], bins=50)
      title1 = 'Mean Separation (Sigma-clipped) of Point-to-Segment Cross-matched sources'
      p1 = [plot_histogram(title1, hist1, edges1, y_start=0,
                           fill_color='navy', background_fill_color='#fafafa',
                           xlabel='Mean Separation of Cross-matched sources (arcseconds)',
                           ylabel='Number of products')]
      plot_list += p1
  
      hist2, edges2 = np.histogram(xmatchCDS.data[data_cols[num_hover_cols + 12]], bins=50)
      title2 = 'Median Separation (Sigma-clipped) of Point-to-Segment Cross-matched sources'
      p2 = [plot_histogram(title2, hist2, edges2, y_start=0,
                           fill_color='navy', background_fill_color='#fafafa',
                           xlabel='Median Separation of Cross-matched sources (arcseconds)',
                           ylabel='Number of products')]
      plot_list += p2
  
      hist3, edges3 = np.histogram(xmatchCDS.data[data_cols[num_hover_cols + 13]], bins=50)
      title3 = 'Standard-deviation (sigma-clipped) of Separation of Point-to-Segment Cross-matched sources'
      p3 = [plot_histogram(title3, hist3, edges3, y_start=0,
                           fill_color='navy', background_fill_color='#fafafa',
                           xlabel='STD(Separation) of Cross-matched sources (arcseconds)',
                           ylabel='Number of products')]
      plot_list += p3
  
      # Save the plot to an HTML file
      save(column(plot_list))
  
      return output
  
  # -----------------------------------------------------------------------------
  # Utility functions for plotting
  #
  
  # MDD Deprecated - only used by build_crossmatch_plots in this module which
  # have not been migrated.  The get_pandas_data in svm_quality_graphics.py
  # should be used on migrated code.  NOTE: The returned value  from the
  # migrated version of get_pandas_data is a DataFrame
  # and not a ColumnDataSource.   I suggest you do not convert the DataFrame
  # into a ColumnDataSource until you are ready to plot as this is really
  # what the ColumnDataSource is there to support (IMHO).
  
  
  def get_pandas_data(data_source, data_columns):
  
      """Load the harvested data, stored in a CSV file, into local arrays.
  
      Parameters
      ==========
      data_source : str
          Name of the file containing the Pandas Dataframe created by the harvester.
  
      data_columns : list
          List of column names which should be extracted from the master file
          `data_source` for use in a plot.
  
      Returns
      =======
      phot_data : Pandas ColumnDataSource
      Dataframe which is a subset of the input Pandas dataframe written out as
      a CSV file.  The subset dataframe consists of only the requested columns
      and rows where all of the requested columns did not contain NaNs.
  
      """
  
      # Instantiate a Pandas Dataframe Reader (lazy instantiation)
      # df_handle = PandasDFReader_CSV("svm_qa_dataframe.csv")
      df_handle = PandasDFReader(data_source, log_level=logutil.logging.NOTSET)
  
      # In this particular case, the names of the desired columns do not
      # have to be further manipulated, for example, to add dataset specific
      # names.
      #
      # Get the relevant column data, eliminating all rows which have NaNs
      # in any of the relevant columns.
      if data_source.endswith(".h5"):
          data_cols = df_handle.get_columns_HDF5(data_columns)
      else:
          data_cols = df_handle.get_columns_CSV(data_columns)
  
      # Setup the source of the data to be plotted so the axis variables can be
      # referenced by column name in the Pandas dataframe
      dataDF = ColumnDataSource(data_cols)
      num_of_datasets = len(dataDF.data['index'])
      print('Number of datasets: {}'.format(num_of_datasets))
  
      return dataDF
  
  
  # MDD Deprecated - only used by build_crossmatch_plots which have not been migrated.
  # This is only still here if you want to run the graphics in this routine to see
  # what the code does.  This routine should NOT be migrated.
  def plot_histogram(title, hist, edges, y_start=0,
                     fill_color='navy', background_fill_color='#fafafa',
                     xlabel='', ylabel=''):
      p = figure(title=title, tools=FIGURE_TOOLS,
                 background_fill_color=background_fill_color)
      p.quad(top=hist, bottom=0, left=edges[:-1], right=edges[1:],
             fill_color=fill_color, line_color="white", alpha=0.5)
  
      p.y_range.start = y_start
      p.legend.location = "center_right"
      p.xaxis.axis_label = xlabel
      p.yaxis.axis_label = ylabel
      p.grid.grid_line_color = "white"
      return p
  