diff --git a/drizzlepac/devutils/search_skyframes.py b/drizzlepac/devutils/search_skyframes.py index d652c0a6e..c0b52d841 100644 --- a/drizzlepac/devutils/search_skyframes.py +++ b/drizzlepac/devutils/search_skyframes.py @@ -101,7 +101,7 @@ def augment_results(results): # populate dateobs_list, path_list for idx in results.index: rootname = results.exposure[idx] - imgname = make_poller_files.locate_fitspath_from_rootname(rootname) + imgname = make_poller_files.locate_fitsfile(rootname) dateobs_list.append(fits.getval(imgname, "DATE-OBS")) path_list.append(imgname) @@ -305,7 +305,7 @@ def make_footprint_fits_file(skycell_name, img_list, footprint_imgname): parser.add_argument('-f', '--spec', required=False, default="None", help='Filter name(s) to search for. To search for ACS observations that use two ' 'spectral elements, enter the names of both spectral elements in any order ' - 'seperated by a dash. Example two-spectral element input: f606w-pol60v') + 'separated by a dash. Example two-spectral element input: f606w-pol60v') parser.add_argument('-m', '--master_observations_file', required=False, default=os.getenv("ALL_EXP_FILE"), help='Name of the master observations .csv file containing comma-separated columns ' diff --git a/drizzlepac/hapmultisequencer.py b/drizzlepac/hapmultisequencer.py index ce361f970..ac1deaf15 100644 --- a/drizzlepac/hapmultisequencer.py +++ b/drizzlepac/hapmultisequencer.py @@ -231,9 +231,10 @@ def create_drizzle_products(total_obj_list, custom_limits=None): # ---------------------------------------------------------------------------------------------------------------------- -def run_mvm_processing(input_filename, diagnostic_mode=False, use_defaults_configs=True, - input_custom_pars_file=None, output_custom_pars_file=None, phot_mode="both", - custom_limits=None, output_file_prefix=None, log_level=logutil.logging.INFO): +def run_mvm_processing(input_filename, skip_gaia_alignment=False, diagnostic_mode=False, + use_defaults_configs=True, input_custom_pars_file=None, output_custom_pars_file=None, + phot_mode="both", custom_limits=None, output_file_prefix=None, + log_level=logutil.logging.INFO): """Run the HST Advanced Products (HAP) generation code. This routine is the sequencer or controller which invokes the high-level functionality to process the multi-visit data. @@ -244,6 +245,10 @@ def run_mvm_processing(input_filename, diagnostic_mode=False, use_defaults_confi The 'poller file' where each line contains information regarding an exposures considered part of the multi-visit. + skip_gaia_alignment : bool, optional + Skip alignment of all input images to known Gaia/HSC sources in the input image footprint? If set to + 'True', the existing input image alignment solution will be used instead. The default is False. + diagnostic_mode : bool, optional Allows printing of additional diagnostic information to the log. Also, can turn on creation and use of pickled information. @@ -355,11 +360,15 @@ def run_mvm_processing(input_filename, diagnostic_mode=False, use_defaults_confi log.info("The configuration parameters have been read and applied to the drizzle objects.") # TODO: This is the place where updated WCS info is migrated from drizzlepac params to filter objects - - reference_catalog = run_align_to_gaia(total_obj_list, custom_limits=custom_limits, - log_level=log_level, diagnostic_mode=diagnostic_mode) - if reference_catalog: - product_list += [reference_catalog] + if skip_gaia_alignment: + log.info("Gaia alignment step skipped. Existing input image alignment solution will be used instead.") + else: + reference_catalog = run_align_to_gaia(total_obj_list, + custom_limits=custom_limits, + log_level=log_level, + diagnostic_mode=diagnostic_mode) + if reference_catalog: + product_list += [reference_catalog] # Run AstroDrizzle to produce drizzle-combined products log.info("\n{}: Create drizzled imagery products.".format(str(datetime.datetime.now()))) diff --git a/drizzlepac/hapsequencer.py b/drizzlepac/hapsequencer.py index 7ada49890..c8b98d8d1 100755 --- a/drizzlepac/hapsequencer.py +++ b/drizzlepac/hapsequencer.py @@ -162,46 +162,6 @@ def create_catalog_products(total_obj_list, log_level, diagnostic_mode=False, ph # images and some of the measurements can be appended to the total catalog total_product_catalogs.identify(mask=total_product_obj.mask) - # Determine how to continue if "aperture" or "segment" fails to find sources for this total - # detection product - take into account the initial setting of phot_mode. - # If no sources were found by either the point or segmentation algorithms, go on to - # the next total detection product (detector) in the visit with the initially requested - # phot_mode. If the point or segmentation algorithms found sources, need to continue - # processing for that (those) algorithm(s) only. - - # When both algorithms have been requested... - if input_phot_mode == 'both': - # If no sources found with either algorithm, skip to the next total detection product - if total_product_catalogs.catalogs['aperture'].sources is None and total_product_catalogs.catalogs['segment'].sources is None: - log.info("No sources found with Segmentation or Point algorithms for TDP {} - skip to next TDP".format(total_product_obj.drizzle_filename)) - del total_product_catalogs.catalogs['aperture'] - del total_product_catalogs.catalogs['segment'] - continue - - # Only point algorithm found sources, continue to the filter catalogs for just point - if total_product_catalogs.catalogs['aperture'].sources is not None and total_product_catalogs.catalogs['segment'].sources is None: - log.info("Sources only found with Point algorithm for TDP {} - phot_mode set only to POINT for this TDP".format(total_product_obj.drizzle_filename)) - phot_mode = 'aperture' - del total_product_catalogs.catalogs['segment'] - - # Only segment algorithm found sources, continue to the filter catalogs for just segmentation - if total_product_catalogs.catalogs['aperture'].sources is None and total_product_catalogs.catalogs['segment'].sources is not None: - log.info("Sources only found with Segmentation algorithm for TDP {} - phot_mode set only to SEGMENT for this TDP".format(total_product_obj.drizzle_filename)) - phot_mode = 'segment' - del total_product_catalogs.catalogs['aperture'] - - # Only requested the point algorithm - elif input_phot_mode == 'aperture': - if total_product_catalogs.catalogs['aperture'].sources is None: - del total_product_catalogs.catalogs['aperture'] - continue - - # Only requested the segmentation algorithm - elif input_phot_mode == 'segment': - if total_product_catalogs.catalogs['segment'].sources is None: - del total_product_catalogs.catalogs['segment'] - continue - # Build dictionary of total_product_catalogs.catalogs[*].sources to use for # filter photometric catalog generation sources_dict = {} @@ -243,7 +203,6 @@ def create_catalog_products(total_obj_list, log_level, diagnostic_mode=False, ph # a filter "subset" table which will be combined with the total detection table. filter_name = filter_product_obj.filters filter_product_catalogs.measure(filter_name) - log.info("Flagging sources in filter product catalog") filter_product_catalogs = run_sourcelist_flagging(filter_product_obj, filter_product_catalogs, @@ -355,30 +314,33 @@ def create_catalog_products(total_obj_list, log_level, diagnostic_mode=False, ph # rate of cosmic-ray contamination for the total detection product reject_catalogs = total_product_catalogs.verify_crthresh(n1_exposure_time) - if not reject_catalogs or diagnostic_mode: - for filter_product_obj in total_product_obj.fdp_list: - filter_product_catalogs = filter_catalogs[filter_product_obj.drizzle_filename] - - # Now write the catalogs out for this filter product - log.info("Writing out filter product catalog") - # Write out photometric (filter) catalog(s) - filter_product_catalogs.write(reject_catalogs) + if diagnostic_mode: + # If diagnostic mode, we want to inspect the original full source catalogs + reject_catalogs = False - # append filter product catalogs to list - if phot_mode in ['aperture', 'both']: - product_list.append(filter_product_obj.point_cat_filename) - if phot_mode in ['segment', 'both']: - product_list.append(filter_product_obj.segment_cat_filename) + for filter_product_obj in total_product_obj.fdp_list: + filter_product_catalogs = filter_catalogs[filter_product_obj.drizzle_filename] - log.info("Writing out total product catalog") - # write out list(s) of identified sources - total_product_catalogs.write(reject_catalogs) + # Now write the catalogs out for this filter product + log.info("Writing out filter product catalog") + # Write out photometric (filter) catalog(s) + filter_product_catalogs.write(reject_catalogs) - # append total product catalogs to manifest list + # append filter product catalogs to list if phot_mode in ['aperture', 'both']: - product_list.append(total_product_obj.point_cat_filename) + product_list.append(filter_product_obj.point_cat_filename) if phot_mode in ['segment', 'both']: - product_list.append(total_product_obj.segment_cat_filename) + product_list.append(filter_product_obj.segment_cat_filename) + + log.info("Writing out total product catalog") + # write out list(s) of identified sources + total_product_catalogs.write(reject_catalogs) + + # append total product catalogs to manifest list + if phot_mode in ['aperture', 'both']: + product_list.append(total_product_obj.point_cat_filename) + if phot_mode in ['segment', 'both']: + product_list.append(total_product_obj.segment_cat_filename) return product_list @@ -871,22 +833,26 @@ def run_sourcelist_flagging(filter_product_obj, filter_product_catalogs, log_lev pickle_out.close() log.info("Wrote hla_flag_filter param pickle file {} ".format(out_pickle_filename)) # TODO: REMOVE ABOVE CODE ONCE FLAGGING PARAMS ARE OPTIMIZED + if catalog_data is not None and len(catalog_data) > 0: + source_cat = hla_flag_filter.run_source_list_flagging(drizzled_image, + flt_list, + param_dict, + exptime, + plate_scale, + median_sky, + catalog_name, + catalog_data, + cat_type, + drz_root_dir, + filter_product_obj.hla_flag_msk, + ci_lookup_file_path, + output_custom_pars_file, + log_level, + diagnostic_mode) + else: + source_cat = catalog_data - filter_product_catalogs.catalogs[cat_type].source_cat = hla_flag_filter.run_source_list_flagging(drizzled_image, - flt_list, - param_dict, - exptime, - plate_scale, - median_sky, - catalog_name, - catalog_data, - cat_type, - drz_root_dir, - filter_product_obj.hla_flag_msk, - ci_lookup_file_path, - output_custom_pars_file, - log_level, - diagnostic_mode) + filter_product_catalogs.catalogs[cat_type].source_cat = source_cat return filter_product_catalogs diff --git a/drizzlepac/haputils/_detection_utils.py b/drizzlepac/haputils/_detection_utils.py new file mode 100644 index 000000000..d87227737 --- /dev/null +++ b/drizzlepac/haputils/_detection_utils.py @@ -0,0 +1,647 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +This module implements common utility functions and classes for the star +finders. This is a copy of photutils private functions, which have been +refactored. Use of these tools can likely be removed after drizzlepac +requires photutils >= 1.2.0. +""" + +import math +import warnings + +from astropy.convolution import Kernel2D +from astropy.stats import gaussian_fwhm_to_sigma +from astropy.units import Quantity +from astropy.utils import lazyproperty +from astropy.utils.exceptions import AstropyUserWarning +import numpy as np +from photutils.detection import find_peaks +from photutils.utils.exceptions import NoDetectionsWarning + + +def _filter_data(data, kernel, mode='constant', fill_value=0.0, + check_normalization=False): + """ + Convolve a 2D image with a 2D kernel. + + The kernel may either be a 2D `~numpy.ndarray` or a + `~astropy.convolution.Kernel2D` object. + + Parameters + ---------- + data : array_like + The 2D array of the image. + + kernel : array-like (2D) or `~astropy.convolution.Kernel2D` + The 2D kernel used to filter the input ``data``. Filtering the + ``data`` will smooth the noise and maximize detectability of + objects with a shape similar to the kernel. + + mode : {'constant', 'reflect', 'nearest', 'mirror', 'wrap'}, optional + The ``mode`` determines how the array borders are handled. For + the ``'constant'`` mode, values outside the array borders are + set to ``fill_value``. The default is ``'constant'``. + + fill_value : scalar, optional + Value to fill data values beyond the array borders if ``mode`` + is ``'constant'``. The default is ``0.0``. + + check_normalization : bool, optional + If `True` then a warning will be issued if the kernel is not + normalized to 1. + """ + from scipy import ndimage + + if kernel is not None: + if isinstance(kernel, Kernel2D): + kernel_array = kernel.array + else: + kernel_array = kernel + + if check_normalization: + if not np.allclose(np.sum(kernel_array), 1.0): + warnings.warn('The kernel is not normalized.', + AstropyUserWarning) + + # scipy.ndimage.convolve currently strips units, but be explicit + # in case that behavior changes + unit = None + if isinstance(data, Quantity): + unit = data.unit + data = data.value + + # NOTE: astropy.convolution.convolve fails with zero-sum + # kernels (used in findstars) (cf. astropy #1647) + # NOTE: if data is int and kernel is float, ndimage.convolve + # will return an int image - here we make the data float so + # that a float image is always returned + result = ndimage.convolve(data.astype(float), kernel_array, + mode=mode, cval=fill_value) + + if unit is not None: + result = result * unit # can't use *= with older astropy + + return result + else: + return data + + +class _StarFinderKernel: + """ + Class to calculate a 2D Gaussian density enhancement kernel. + + The kernel has negative wings and sums to zero. It is used by both + `DAOStarFinder` and `IRAFStarFinder`. + + Parameters + ---------- + fwhm : float + The full-width half-maximum (FWHM) of the major axis of the + Gaussian kernel in units of pixels. + + ratio : float, optional + The ratio of the minor and major axis standard deviations of the + Gaussian kernel. ``ratio`` must be strictly positive and less + than or equal to 1.0. The default is 1.0 (i.e., a circular + Gaussian kernel). + + theta : float, optional + The position angle (in degrees) of the major axis of the + Gaussian kernel, measured counter-clockwise from the positive x + axis. + + sigma_radius : float, optional + The truncation radius of the Gaussian kernel in units of sigma + (standard deviation) [``1 sigma = FWHM / + 2.0*sqrt(2.0*log(2.0))``]. The default is 1.5. + + normalize_zerosum : bool, optional + Whether to normalize the Gaussian kernel to have zero sum, The + default is `True`, which generates a density-enhancement kernel. + + Notes + ----- + The class attributes include the dimensions of the elliptical kernel + and the coefficients of a 2D elliptical Gaussian function expressed + as: + + ``f(x,y) = A * exp(-g(x,y))`` + + where + + ``g(x,y) = a*(x-x0)**2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)**2`` + + References + ---------- + .. [1] https://en.wikipedia.org/wiki/Gaussian_function + """ + + def __init__(self, fwhm, ratio=1.0, theta=0.0, sigma_radius=1.5, + normalize_zerosum=True): + + if fwhm < 0: + raise ValueError('fwhm must be positive.') + + if ratio <= 0 or ratio > 1: + raise ValueError('ratio must be positive and less or equal ' + 'than 1.') + + if sigma_radius <= 0: + raise ValueError('sigma_radius must be positive.') + + self.fwhm = fwhm + self.ratio = ratio + self.theta = theta + self.sigma_radius = sigma_radius + self.xsigma = self.fwhm * gaussian_fwhm_to_sigma + self.ysigma = self.xsigma * self.ratio + + theta_radians = np.deg2rad(self.theta) + cost = np.cos(theta_radians) + sint = np.sin(theta_radians) + xsigma2 = self.xsigma**2 + ysigma2 = self.ysigma**2 + + self.a = (cost**2 / (2.0 * xsigma2)) + (sint**2 / (2.0 * ysigma2)) + # CCW + self.b = 0.5 * cost * sint * ((1.0 / xsigma2) - (1.0 / ysigma2)) + self.c = (sint**2 / (2.0 * xsigma2)) + (cost**2 / (2.0 * ysigma2)) + + # find the extent of an ellipse with radius = sigma_radius*sigma; + # solve for the horizontal and vertical tangents of an ellipse + # defined by g(x,y) = f + self.f = self.sigma_radius**2 / 2.0 + denom = (self.a * self.c) - self.b**2 + + # nx and ny are always odd + self.nx = 2 * int(max(2, math.sqrt(self.c * self.f / denom))) + 1 + self.ny = 2 * int(max(2, math.sqrt(self.a * self.f / denom))) + 1 + + self.xc = self.xradius = self.nx // 2 + self.yc = self.yradius = self.ny // 2 + + # define the kernel on a 2D grid + yy, xx = np.mgrid[0:self.ny, 0:self.nx] + self.circular_radius = np.sqrt((xx - self.xc)**2 + (yy - self.yc)**2) + self.elliptical_radius = (self.a * (xx - self.xc)**2 + + 2.0 * self.b * (xx - self.xc) + * (yy - self.yc) + + self.c * (yy - self.yc)**2) + + self.mask = np.where( + (self.elliptical_radius <= self.f) + | (self.circular_radius <= 2.0), 1, 0).astype(int) + self.npixels = self.mask.sum() + + # NOTE: the central (peak) pixel of gaussian_kernel has a value of 1. + self.gaussian_kernel_unmasked = np.exp(-self.elliptical_radius) + self.gaussian_kernel = self.gaussian_kernel_unmasked * self.mask + + # denom = variance * npixels + denom = ((self.gaussian_kernel**2).sum() + - (self.gaussian_kernel.sum()**2 / self.npixels)) + self.relerr = 1.0 / np.sqrt(denom) + + # normalize the kernel to zero sum + if normalize_zerosum: + self.data = ((self.gaussian_kernel + - (self.gaussian_kernel.sum() / self.npixels)) + / denom) * self.mask + else: + self.data = self.gaussian_kernel + + self.shape = self.data.shape + + +class _StarCutout: + """ + Class to hold a 2D image cutout of a single star for the star finder + classes. + + Parameters + ---------- + data : 2D array_like + The cutout 2D image from the input unconvolved 2D image. + + convdata : 2D array_like + The cutout 2D image from the convolved 2D image. + + slices : tuple of two slices + A tuple of two slices representing the minimal box of the cutout + from the original image. + + xpeak, ypeak : float + The (x, y) pixel coordinates of the peak pixel. + + kernel : `_StarFinderKernel` + The convolution kernel. The shape of the kernel must match that + of the input ``data``. + + threshold_eff : float + The absolute image value above which to select sources. This + threshold should be the threshold value input to the star finder + class multiplied by the kernel relerr. + """ + + def __init__(self, data, convdata, slices, xpeak, ypeak, kernel, + threshold_eff): + + self.data = data + self.convdata = convdata + self.slices = slices + self.xpeak = xpeak + self.ypeak = ypeak + self.kernel = kernel + self.threshold_eff = threshold_eff + + self.shape = data.shape + self.nx = self.shape[1] # always odd + self.ny = self.shape[0] # always odd + self.cutout_xcenter = int(self.nx // 2) + self.cutout_ycenter = int(self.ny // 2) + + self.xorigin = self.slices[1].start # in original image + self.yorigin = self.slices[0].start # in original image + + self.mask = kernel.mask # kernel mask + self.npixels = kernel.npixels # unmasked pixels + self.data_masked = self.data * self.mask + + +def _find_stars(data, kernel, threshold_eff, min_separation=None, + mask=None, exclude_border=False): + """ + Find stars in an image. + + Parameters + ---------- + data : 2D array_like + The 2D array of the image. + + kernel : `_StarFinderKernel` + The convolution kernel. + + threshold_eff : float + The absolute image value above which to select sources. This + threshold should be the threshold input to the star finder class + multiplied by the kernel relerr. + + min_separation : float, optional + The minimum separation for detected objects in pixels. + + mask : 2D bool array, optional + A boolean mask with the same shape as ``data``, where a `True` + value indicates the corresponding element of ``data`` is masked. + Masked pixels are ignored when searching for stars. + + exclude_border : bool, optional + Set to `True` to exclude sources found within half the size of + the convolution kernel from the image borders. The default is + `False`, which is the mode used by IRAF's `DAOFIND`_ and + `starfind`_ tasks. + + Returns + ------- + objects : list of `_StarCutout` + A list of `_StarCutout` objects containing the image cutout for + each source. + + .. _DAOFIND: https://iraf.net/irafhelp.php?val=daofind + + .. _starfind: https://iraf.net/irafhelp.php?val=starfind + """ + convolved_data = _filter_data(data, kernel.data, mode='constant', + fill_value=0.0, check_normalization=False) + + # define a local footprint for the peak finder + if min_separation is None: # daofind + footprint = kernel.mask.astype(bool) + else: + # define a circular footprint + idx = np.arange(-min_separation, min_separation + 1) + xx, yy = np.meshgrid(idx, idx) + footprint = np.array((xx**2 + yy**2) <= min_separation**2, dtype=int) + + # pad the data and convolved image by the kernel x/y radius to allow + # for detections near the edges + if not exclude_border: + ypad = kernel.yradius + xpad = kernel.xradius + pad = ((ypad, ypad), (xpad, xpad)) + pad_mode = 'constant' + data = np.pad(data, pad, mode=pad_mode, constant_values=[0.]) + if mask is not None: + mask = np.pad(mask, pad, mode=pad_mode, constant_values=[0.]) + convolved_data = np.pad(convolved_data, pad, mode=pad_mode, + constant_values=[0.]) + + # find local peaks in the convolved data + with warnings.catch_warnings(): + # suppress any NoDetectionsWarning from find_peaks + warnings.filterwarnings('ignore', category=NoDetectionsWarning) + tbl = find_peaks(convolved_data, threshold_eff, footprint=footprint, + mask=mask) + + if tbl is None: + return None + + coords = np.transpose([tbl['y_peak'], tbl['x_peak']]) + + star_cutouts = [] + for (ypeak, xpeak) in coords: + # now extract the object from the data, centered on the peak + # pixel in the convolved image, with the same size as the kernel + x0 = xpeak - kernel.xradius + x1 = xpeak + kernel.xradius + 1 + y0 = ypeak - kernel.yradius + y1 = ypeak + kernel.yradius + 1 + + if x0 < 0 or x1 > data.shape[1]: + continue # pragma: no cover + if y0 < 0 or y1 > data.shape[0]: + continue # pragma: no cover + + slices = (slice(y0, y1), slice(x0, x1)) + data_cutout = data[slices] + convdata_cutout = convolved_data[slices] + + # correct pixel values for the previous image padding + if not exclude_border: + x0 -= kernel.xradius + x1 -= kernel.xradius + y0 -= kernel.yradius + y1 -= kernel.yradius + xpeak -= kernel.xradius + ypeak -= kernel.yradius + slices = (slice(y0, y1), slice(x0, x1)) + + star_cutouts.append(_StarCutout(data_cutout, convdata_cutout, slices, + xpeak, ypeak, kernel, threshold_eff)) + + return star_cutouts + + +class _DAOFindProperties: + """ + Class to calculate the properties of each detected star, as defined + by `DAOFIND`_. + + Parameters + ---------- + star_cutout : `_StarCutout` + A `_StarCutout` object containing the image cutout for the star. + + kernel : `_StarFinderKernel` + The convolution kernel. The shape of the kernel must match that + of the input ``star_cutout``. + + sky : float, optional + The local sky level around the source. ``sky`` is used only to + calculate the source peak value, flux, and magnitude. The + default is 0. + + .. _DAOFIND: https://iraf.net/irafhelp.php?val=daofind + """ + + def __init__(self, star_cutout, kernel, sky=0.): + if not isinstance(star_cutout, _StarCutout): + raise ValueError('data must be an _StarCutout object') + + if star_cutout.data.shape != kernel.shape: + raise ValueError('cutout and kernel must have the same shape') + + self.cutout = star_cutout + self.kernel = kernel + self.sky = sky # DAOFIND has no sky input -> same as sky=0. + + self.data = star_cutout.data + self.data_masked = star_cutout.data_masked + self.npixels = star_cutout.npixels # unmasked pixels + self.nx = star_cutout.nx + self.ny = star_cutout.ny + self.xcenter = star_cutout.cutout_xcenter + self.ycenter = star_cutout.cutout_ycenter + + @lazyproperty + def data_peak(self): + return self.data[self.ycenter, self.xcenter] + + @lazyproperty + def conv_peak(self): + return self.cutout.convdata[self.ycenter, self.xcenter] + + @lazyproperty + def roundness1(self): + # set the central (peak) pixel to zero + cutout_conv = self.cutout.convdata.copy() + cutout_conv[self.ycenter, self.xcenter] = 0.0 # for sum4 + + # calculate the four roundness quadrants. + # the cutout size always matches the kernel size, which have odd + # dimensions. + # quad1 = bottom right + # quad2 = bottom left + # quad3 = top left + # quad4 = top right + # 3 3 4 4 4 + # 3 3 4 4 4 + # 3 3 x 1 1 + # 2 2 2 1 1 + # 2 2 2 1 1 + quad1 = cutout_conv[0:self.ycenter + 1, self.xcenter + 1:] + quad2 = cutout_conv[0:self.ycenter, 0:self.xcenter + 1] + quad3 = cutout_conv[self.ycenter:, 0:self.xcenter] + quad4 = cutout_conv[self.ycenter + 1:, self.xcenter:] + + sum2 = -quad1.sum() + quad2.sum() - quad3.sum() + quad4.sum() + if sum2 == 0: + return 0. + + sum4 = np.abs(cutout_conv).sum() + if sum4 <= 0: + return None + + return 2.0 * sum2 / sum4 + + @lazyproperty + def sharpness(self): + npixels = self.npixels - 1 # exclude the peak pixel + data_mean = (np.sum(self.data_masked) - self.data_peak) / npixels + + return (self.data_peak - data_mean) / self.conv_peak + + def daofind_marginal_fit(self, axis=0): + """ + Fit 1D Gaussians, defined from the marginal x/y kernel + distributions, to the marginal x/y distributions of the original + (unconvolved) image. + + These fits are used calculate the star centroid and roundness + ("GROUND") properties. + + Parameters + ---------- + axis : {0, 1}, optional + The axis for which the marginal fit is performed: + + * 0: for the x axis + * 1: for the y axis + + Returns + ------- + dx : float + The fractional shift in x or y (depending on ``axis`` value) + of the image centroid relative to the maximum pixel. + + hx : float + The height of the best-fitting Gaussian to the marginal x or + y (depending on ``axis`` value) distribution of the + unconvolved source data. + """ + + # define triangular weighting functions along each axis, peaked + # in the middle and equal to one at the edge + x = self.xcenter - np.abs(np.arange(self.nx) - self.xcenter) + 1 + y = self.ycenter - np.abs(np.arange(self.ny) - self.ycenter) + 1 + xwt, ywt = np.meshgrid(x, y) + + if axis == 0: # marginal distributions along x axis + wt = xwt[0] # 1D + wts = ywt # 2D + size = self.nx + center = self.xcenter + sigma = self.kernel.xsigma + dxx = center - np.arange(size) + elif axis == 1: # marginal distributions along y axis + wt = np.transpose(ywt)[0] # 1D + wts = xwt # 2D + size = self.ny + center = self.ycenter + sigma = self.kernel.ysigma + dxx = np.arange(size) - center + + # compute marginal sums for given axis + wt_sum = np.sum(wt) + dx = center - np.arange(size) + + # weighted marginal sums + kern_sum_1d = np.sum(self.kernel.gaussian_kernel_unmasked * wts, + axis=axis) + kern_sum = np.sum(kern_sum_1d * wt) + kern2_sum = np.sum(kern_sum_1d**2 * wt) + + dkern_dx = kern_sum_1d * dx + dkern_dx_sum = np.sum(dkern_dx * wt) + dkern_dx2_sum = np.sum(dkern_dx**2 * wt) + kern_dkern_dx_sum = np.sum(kern_sum_1d * dkern_dx * wt) + + data_sum_1d = np.sum(self.data * wts, axis=axis) + data_sum = np.sum(data_sum_1d * wt) + data_kern_sum = np.sum(data_sum_1d * kern_sum_1d * wt) + data_dkern_dx_sum = np.sum(data_sum_1d * dkern_dx * wt) + data_dx_sum = np.sum(data_sum_1d * dxx * wt) + + # perform linear least-squares fit (where data = sky + hx*kernel) + # to find the amplitude (hx) + # reject the star if the fit amplitude is not positive + hx_numer = data_kern_sum - (data_sum * kern_sum) / wt_sum + if hx_numer <= 0.: + return np.nan, np.nan + + hx_denom = kern2_sum - (kern_sum**2 / wt_sum) + if hx_denom <= 0.: + return np.nan, np.nan + + # compute fit amplitude + hx = hx_numer / hx_denom + # sky = (data_sum - (hx * kern_sum)) / wt_sum + + # compute centroid shift + dx = ((kern_dkern_dx_sum + - (data_dkern_dx_sum - dkern_dx_sum * data_sum)) + / (hx * dkern_dx2_sum / sigma**2)) + + hsize = size / 2. + if abs(dx) > hsize: + if data_sum == 0.: + dx = 0.0 + else: + dx = data_dx_sum / data_sum + if abs(dx) > hsize: + dx = 0.0 + + return dx, hx + + @lazyproperty + def dx_hx(self): + return self.daofind_marginal_fit(axis=0) + + @lazyproperty + def dy_hy(self): + return self.daofind_marginal_fit(axis=1) + + @lazyproperty + def dx(self): + return self.dx_hx[0] + + @lazyproperty + def dy(self): + return self.dy_hy[0] + + @lazyproperty + def xcentroid(self): + return self.cutout.xpeak + self.dx + + @lazyproperty + def ycentroid(self): + return self.cutout.ypeak + self.dy + + @lazyproperty + def hx(self): + return self.dx_hx[1] + + @lazyproperty + def hy(self): + return self.dy_hy[1] + + @lazyproperty + def roundness2(self): + """ + The star roundness. + + This roundness parameter represents the ratio of the difference + in the height of the best fitting Gaussian function in x minus + the best fitting Gaussian function in y, divided by the average + of the best fitting Gaussian functions in x and y. A circular + source will have a zero roundness. A source extended in x or y + will have a negative or positive roundness, respectively. + """ + + if np.isnan(self.hx) or np.isnan(self.hy): + return np.nan + else: + return 2.0 * (self.hx - self.hy) / (self.hx + self.hy) + + @lazyproperty + def peak(self): + return self.data_peak - self.sky + + @lazyproperty + def npix(self): + """ + The total number of pixels in the rectangular cutout image. + """ + + return self.data.size + + @lazyproperty + def flux(self): + return ((self.conv_peak / self.cutout.threshold_eff) + - (self.sky * self.npix)) + + @lazyproperty + def mag(self): + if self.flux <= 0: + return np.nan + else: + return -2.5 * np.log10(self.flux) diff --git a/drizzlepac/haputils/align_utils.py b/drizzlepac/haputils/align_utils.py index 740484c33..95683bd19 100755 --- a/drizzlepac/haputils/align_utils.py +++ b/drizzlepac/haputils/align_utils.py @@ -169,6 +169,25 @@ def __init__(self, input_list, clobber=False, dqname='DQ', rms_estimator=self.alignment_pars['rms_estimator'], threshold_flag=self.alignment_pars['threshold']) catimg.build_kernel(fwhmpsf) + catimg.crclean = self.alignment_pars['classify'] + log.info("CATIMG.CRCLEAN: {}".format(catimg.crclean)) + if catimg.crclean: + log.info("Recomputing BACKGROUND after applying single-image CR clean") + for chip in range(1, catimg.num_sci + 1): + sciarr = amutils.crclean_image(catimg.imghdu[("SCI", chip)].data, catimg.threshold[chip], + catimg.kernel, catimg.kernel_fwhm, background=catimg.bkg[chip]) + catimg.imghdu[("SCI", chip)].data = sciarr + # recompute now that the CRs have been removed (set to 0) from the science arrays + catimg.compute_background(box_size=self.alignment_pars['box_size'], + win_size=self.alignment_pars['win_size'], + nsigma=self.alignment_pars['nsigma'], + bkg_estimator=self.alignment_pars['bkg_estimator'], + rms_estimator=self.alignment_pars['rms_estimator'], + threshold_flag=self.alignment_pars['threshold']) + log.info("Finished computing revised BACKROUND") + catimg.build_kernel(fwhmpsf) + log.info("Finished determining revised kernel") + # Use FWHM from first good exposure as default for remainder of exposures if not default_fwhm_set and catimg.kernel is not None: fwhmpsf = catimg.fwhmpsf @@ -666,7 +685,8 @@ def find_alignment_sources(self, output=True, dqname='DQ', crclean=False, dqmask = self.build_dqmask(chip=chip) sciarr = self.imghdu[("SCI", chip)].data.copy() # TODO: replace detector_pars with dict from OO Config class - extract_pars = {'classify': alignment_pars['classify'], + # Turning off 'classify' since same CRs are being removed before segmentation now + extract_pars = {'classify': False, # alignment_pars['classify'], 'centering_mode': alignment_pars['centering_mode'], 'nlargest': alignment_pars['MAX_SOURCES_PER_CHIP'], 'deblend': alignment_pars['deblend']} @@ -738,10 +758,12 @@ def match_relative_fit(imglist, reference_catalog, **fit_pars): else: fitgeom = 'rscale' + rel_fitgeom = 'rscale' + common_pars = fit_pars['pars'] del fit_pars['pars'] - nclip = None if fitgeom == 'shift' else 3 + nclip = 1 if fitgeom == 'rscale' else 0 # Only perform sigma-clipping for 'rscale' # 0: Specify matching algorithm to use match = tweakwcs.TPMatch(**fit_pars) @@ -757,10 +779,10 @@ def match_relative_fit(imglist, reference_catalog, **fit_pars): # limits. match_relcat = tweakwcs.align_wcs(imglist, None, match=match, - minobj=common_pars['minobj'][fitgeom], + minobj=common_pars['minobj'][rel_fitgeom], expand_refcat=True, - fitgeom=fitgeom, - nclip=nclip) + fitgeom=rel_fitgeom, + nclip=1) # Implement a consistency check even before trying absolute alignment # If relative alignment in question, no use in aligning to GAIA if not check_consistency(imglist): @@ -782,7 +804,7 @@ def match_relative_fit(imglist, reference_catalog, **fit_pars): msg = "Image {} --".format(i.meta['name']) msg += "\n SHIFT:({:9.4f},{:9.4f}) NMATCHES: {} ".format(off[0], off[1], nmatches) msg += "\n ROT:{:9.4f} SCALE:{:9.4f}".format(rot, scale) - msg += "\n Using fitgeom = '{}'".format(fitgeom) + msg += "\n Using fitgeom = '{}'".format(rel_fitgeom) log.info(msg) # This logic enables performing only relative fitting and skipping fitting to GAIA @@ -853,7 +875,7 @@ def match_default_fit(imglist, reference_catalog, **fit_pars): common_pars = fit_pars['pars'] del fit_pars['pars'] - nclip = None if fitgeom == 'shift' else 3 + nclip = 1 if fitgeom == 'rscale' else 0 # Only perform sigma-clipping for 'rscale' log.info("{} (match_default_fit) Cross matching and fitting " "{}".format("-" * 20, "-" * 27)) @@ -917,7 +939,7 @@ def match_2dhist_fit(imglist, reference_catalog, **fit_pars): common_pars = fit_pars['pars'] del fit_pars['pars'] - nclip = None if fitgeom == 'shift' else 3 + nclip = 1 if fitgeom == 'rscale' else 0 # Only perform sigma-clipping for 'rscale' log.info("{} (match_2dhist_fit) Cross matching and fitting " "{}".format("-" * 20, "-" * 28)) @@ -961,14 +983,15 @@ def check_consistency(imglist, rot_tolerance=0.1, shift_tolerance=1.0): rots[i] = finfo['proper_rot'] nmatches[i] = finfo['nmatches'] else: - # Set default as value larger than we can use - nmatches[i] = 1000000 + # 'status' == 'FAILED': Set default as negative value + nmatches[i] = -1 + is_consistent = False # We should only need to check for consistency when less than 5 # matches were used for the fit, leading to a higher potential for # a singular solution or mis-identified cross-matches. if nmatches.min() > 4: - return imglist + return is_consistent # compute deltas to look for outliers delta_rots = rots[1:] - rots[:-1] diff --git a/drizzlepac/haputils/analyze.py b/drizzlepac/haputils/analyze.py index 49ca315de..fb458c5e7 100644 --- a/drizzlepac/haputils/analyze.py +++ b/drizzlepac/haputils/analyze.py @@ -1,10 +1,15 @@ -""" Utility to analyze an input dataset and determine whether the dataset can be aligned +""" Utilities to analyze an input and determine whether the input is viable for a given process -The function analyze_data opens an input list containing FLT and/or FLC FITS filenames -in order to access the primary header data. Based upon the values of specific +The fundamental function analyze_data opens an input list containing FLT and/or FLC FITS +filenames in order to access the primary header data. Based upon the values of specific FITS keywords, the function determines whether or not each file within this dataset can or should be reconciled against an astrometric catalog and, for multiple images, used to create a mosaic. + +The functions, mvm_analyze_wrapper and analyze_wrapper, are thin wrappers around analyze_data +to accommodate special case uses. The mvm_analyze_wrapper takes a filename as input and +returns a boolean indicator. The analyze_wrapper has the same function signature as +analyze_data, but it returns a list instead of an astropy table. """ import math import sys @@ -20,10 +25,10 @@ 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, +log = logutil.create_logger(__name__, level=logutil.logging.DEBUG, stream=sys.stdout, format=SPLUNK_MSG_FORMAT, datefmt=MSG_DATEFMT) -__all__ = ['analyze_data'] +__all__ = ['analyze_data', 'analyze_wrapper', 'mvm_analyze_wrapper'] # Define global default keyword names for these fields OBSKEY = 'OBSTYPE' @@ -51,7 +56,41 @@ class Messages(Enum): OK, WARN, NOPROC = 1, -1, -2 -def analyze_wrapper(input_file_list, log_level=logutil.logging.NOTSET): +def mvm_analyze_wrapper(input_filename, log_level=logutil.logging.DEBUG): + """ + Thin wrapper for the analyze_data function to return a viability indicator regarding a image for MVM processing. + + Parameters + ========== + input_filename : string + Full fileName of data to be analyzed for viability to be processed as an MVM constituent. + + Returns + ======= + use_for_mvm : boolean + Boolean which indicates whether the input image should be used for MVM processing - + True: use for MVM, False: do NOT use for MVM + """ + # Set logging level to user-specified level + log.setLevel(log_level) + + # Invoke the low-level analyze_data routine with type = "MVM" + filtered_table = analyze_data([input_filename], type = "MVM") + + # There is only one row in this output table + use_for_mvm = False + if filtered_table['doProcess'] == 0: + use_for_mvm = False + log.warning("Image, {}, cannot be used for MVM processing. Issue: {}.\n". \ + format(input_filename, filtered_table['processMsg'][0])) + else: + use_for_mvm = True + log.info("Image, {}, will be used for MVM processing.".format(input_filename)) + + return use_for_mvm + + +def analyze_wrapper(input_file_list, log_level=logutil.logging.DEBUG): """ Thin wrapper for the analyze_data function to return a list of viable images. @@ -70,7 +109,7 @@ def analyze_wrapper(input_file_list, log_level=logutil.logging.NOTSET): This routine returns a list containing only viable images instead of a table which provides information, as well as a doProcess bool, regarding each image. """ - # set logging level to user-specified level + # Set logging level to user-specified level log.setLevel(log_level) process_list = [] @@ -89,7 +128,7 @@ def analyze_wrapper(input_file_list, log_level=logutil.logging.NOTSET): return process_list -def analyze_data(input_file_list, log_level=logutil.logging.NOTSET): +def analyze_data(input_file_list, log_level=logutil.logging.DEBUG, type=""): """ Determine if images within the dataset can be aligned @@ -104,6 +143,12 @@ def analyze_data(input_file_list, log_level=logutil.logging.NOTSET): The desired level of verboseness in the log statements displayed on the screen and written to the .log file. Default value is 20, or 'info'. + type : string + String indicating whether this file is for MVM or some other processing. + If type == "MVM", then Grism/Prism data is ignored. If type == "" (default) or any + other string, the Grism/Prism data is considered available for processing unless there is + some other issue (i.e., exposure time of zero). + Returns ======= output_table : object @@ -134,7 +179,7 @@ def analyze_data(input_file_list, log_level=logutil.logging.NOTSET): Please be aware of the FITS keyword value NONE vs the Python None. """ - # set logging level to user-specified level + # Set logging level to user-specified level log.setLevel(log_level) acs_filt_name_list = [FILKEY1, FILKEY2] @@ -303,12 +348,20 @@ def analyze_data(input_file_list, log_level=logutil.logging.NOTSET): # for all the good data. split_sfilter = sfilter.upper().split('_') for item in split_sfilter: - if item.startswith(('G', 'PR')) and not is_zero: + # This is the only circumstance when Grism/Prism data WILL be processed. + if item.startswith(('G', 'PR')) and not is_zero and type.upper() != "MVM": no_proc_key = None no_proc_value = None - elif item.startswith(('G', 'PR')) and is_zero: - no_proc_value += " and EXPTIME = 0.0" - + log.info("The Grism/Prism data, {}, will be processed.".format(input_file)) + # Grism/Prism WILL NOT be processed primarily if MVM processing or with an exposure time of zero. + elif item.startswith(('G', 'PR')): + if type.upper() == "MVM": + no_proc_value += ", Grism/Prism data and MVM processing" + log.warning("The Grism/Prism data {} with MVM processing will be ignored.".format(input_file)) + elif is_zero: + no_proc_value += ", Grism/Prism data and EXPTIME = 0.0" + log.warning("The Grism/Prism data {} with zero exposure time will be ignored.".format(input_file)) + if item.startswith(('BLOCK')): no_proc_key = FILKEY no_proc_value = sfilter diff --git a/drizzlepac/haputils/astrometric_utils.py b/drizzlepac/haputils/astrometric_utils.py index 503cdbf38..7168525c3 100644 --- a/drizzlepac/haputils/astrometric_utils.py +++ b/drizzlepac/haputils/astrometric_utils.py @@ -20,6 +20,7 @@ import sys import time from distutils.version import LooseVersion +import copy import numpy as np import scipy.stats as st @@ -295,6 +296,7 @@ def create_astrometric_catalog(inputs, catalog="GAIAedr3", output="ref_cat.ecsv" ref_table.meta['catalog'] = catalog ref_table.meta['gaia_only'] = gaia_only ref_table.meta['epoch'] = epoch + ref_table.meta['converted'] = False # Convert common set of columns into standardized column names convert_astrometric_table(ref_table, catalog) @@ -1001,6 +1003,7 @@ def extract_sources(img, dqmask=None, fwhm=3.0, kernel=None, photmode=None, # If classify is turned on, it should modify the segmentation map dqmap = None if classify: + log.info('Removing suspected cosmic-rays from source catalog') cat = SourceCatalog(imgarr, segm) # Remove likely cosmic-rays based on central_moments classification bad_srcs = np.where(classify_sources(cat, fwhm) == 0)[0] + 1 @@ -1157,6 +1160,56 @@ def extract_sources(img, dqmask=None, fwhm=3.0, kernel=None, photmode=None, return tbl, segm, dqmap +def crclean_image(imgarr, segment_threshold, kernel, fwhm, + source_box=7, deblend=False, background=0.0): + """ Apply moments-based single-image CR clean algorithm to image + + Returns + ------- + imgarr : ndarray + Input array with pixels flagged as CRs set to 0.0 + """ + # Perform basic image segmentation to look for sources in image + segm = detect_sources(imgarr, segment_threshold, npixels=source_box, + filter_kernel=kernel, connectivity=4) + + log.debug("Detected {} sources of all types ".format(segm.nlabels)) + # photutils >= 0.7: segm=None; photutils < 0.7: segm.nlabels=0 + if segm is None or segm.nlabels == 0: + log.info("No sources to identify as cosmic-rays!") + return imgarr + + if deblend: + segm = deblend_sources(imgarr, segm, npixels=5, + filter_kernel=kernel, nlevels=32, + contrast=0.01) + + # Measure segment sources to compute centralized moments for each source + cat = SourceCatalog(imgarr, segm) + # Remove likely cosmic-rays based on central_moments classification + bad_srcs = classify_sources(cat, fwhm) + # Get the label IDs for sources flagged as CRs, IDs are 1-based not 0-based + bad_ids = np.where(bad_srcs == 0)[0] + 1 + log.debug("Recognized {} sources as cosmic-rays".format(len(bad_ids))) + + # Instead of ignoring them by removing them from the catalog + # zero all these sources out from the image itself. + # Create a mask from the segmentation image with only the flagged CRs + cr_segm = copy.deepcopy(segm) + cr_segm.keep_labels(bad_ids) + cr_mask = np.where(cr_segm.data > 0) + # Use the mask to replace all pixels associated with the CRs + # with values specified as the background + if not isinstance(background, np.ndarray): + background = np.ones_like(imgarr, dtype=np.float32) * background + + imgarr[cr_mask] = background[cr_mask] + del cr_segm + + log.debug("Finshed zeroing out cosmic-rays") + + return imgarr + def classify_sources(catalog, fwhm, sources=None): """ Convert moments_central attribute for source catalog into star/cr flag. @@ -1180,29 +1233,37 @@ def classify_sources(catalog, fwhm, sources=None): An ndarray where a value of 1 indicates a likely valid, non-cosmic-ray source, and a value of 0 indicates a likely cosmic-ray. """ - moments = catalog.moments_central - semiminor_axis = catalog.semiminor_sigma - elon = catalog.elongation + # All these return ndarray objects + moments = catalog.moments_central # (num_sources, 4, 4) + semiminor_axis = catalog.semiminor_sigma.to_value() # (num_sources,) + elon = catalog.elongation.to_value() # (num_sources,) + # Determine how many sources we have. if sources is None: sources = (0, len(moments)) num_sources = sources[1] - sources[0] srctype = np.zeros((num_sources,), np.int32) - for src in range(sources[0], sources[1]): - # Protect against spurious detections - src_x = catalog[src].xcentroid - src_y = catalog[src].ycentroid - if np.isnan(src_x) or np.isnan(src_y): - continue - # This identifies moment of maximum value - x, y = np.where(moments[src] == moments[src].max()) - valid_src = (x[0] > 1) and (y[0] > 1) - # These look for CR streaks (not delta CRs) - valid_width = semiminor_axis[src].value < (0.75 * fwhm) # skinny source - valid_elon = elon[src].value > 2 # long source - valid_streak = valid_width and valid_elon # long and skinny... - # If either a delta CR or a CR streak are identified, remove it - if valid_src and not valid_streak: - srctype[src] = 1 + + # Find which moment is the peak moment for all sources + mx = np.array([np.where(moments[src] == moments[src].max())[0][0] for src in range(num_sources)], np.int32) + my = np.array([np.where(moments[src] == moments[src].max())[1][0] for src in range(num_sources)], np.int32) + # Now avoid processing sources which had NaN as either/both of their X,Y positions (rare, but does happen) + src_x = catalog.xcentroid # (num_sources, ) + src_y = catalog.ycentroid # (num_sources, ) + valid_xy = ~np.bitwise_or(np.isnan(src_x), np.isnan(src_y)) + + # Define descriptors of CRs + # First: look for sources where the moment in X and Y is the first moment + valid_src = np.bitwise_and(mx > 1, my > 1)[valid_xy] + # Second: look for sources where the width of the source is less than the PSF FWHM + valid_width = (semiminor_axis < (0.75 * fwhm))[valid_xy] + # Third: identify sources which are elongated at least 2:1 + valid_elon = (elon > 2)[valid_xy] + # Now combine descriptors into a single value + valid_streak = ~np.bitwise_and(valid_width, valid_elon) + src_cr = np.bitwise_and(valid_src, valid_streak) + # Flag identified sources as CRs with a value of 1 + srctype[src_cr] = 1 + return srctype diff --git a/drizzlepac/haputils/astroquery_utils.py b/drizzlepac/haputils/astroquery_utils.py index fcb76ec03..c4f9ebb4e 100644 --- a/drizzlepac/haputils/astroquery_utils.py +++ b/drizzlepac/haputils/astroquery_utils.py @@ -47,6 +47,17 @@ def retrieve_observation(obsid, suffix=['FLC'], archive=False, clobber=False, clobber : Boolean, optional Download and Overwrite existing files? Default is "False". + product_type : str, optional + Specify what type of product you want from the archive, either 'pipeline' + or 'HAP' or 'both' (default). By default, all versions of the products + processed for the requested datasets will be returned. This would include: + - pipeline : files processed by `runastrodriz` to include the latest + distortion calibrations and the best possible alignment to GAIA + with `ipppssoot_fl[tc].fits` filenames for FLT/FLC files. + - HAP : files processed as a single visit and aligned (as possible) to GAIA + with `hst_______fl[tc].fits` + filenames. + Returns ------- local_files : list diff --git a/drizzlepac/haputils/catalog_utils.py b/drizzlepac/haputils/catalog_utils.py index b1665a69f..299aeb234 100755 --- a/drizzlepac/haputils/catalog_utils.py +++ b/drizzlepac/haputils/catalog_utils.py @@ -577,7 +577,6 @@ def __init__(self, fitsfile, param_dict, param_dict_qc, num_images_mask, log_lev # The syntax here is EXTREMELY cludgy, but until a more compact way to do this is found, # it will have to do... self.catalogs = {} - if 'segment' in self.types: self.catalogs['segment'] = HAPSegmentCatalog(self.image, self.param_dict, self.param_dict_qc, self.diagnostic_mode, tp_sources=tp_sources) @@ -678,7 +677,6 @@ def write(self, reject_catalogs, **pars): Supported types of catalogs include: 'aperture', 'segment'. """ # Make sure we at least have a default 2D background computed - for catalog in self.catalogs.values(): if catalog.source_cat is None: catalog.source_cat = catalog.sources @@ -1026,6 +1024,7 @@ def identify_sources(self, **pars): if not sources: log.warning("No point sources were found in Total Detection Product, {}.".format(self.imgname)) log.warning("Processing for point source catalogs for this product is ending.") + self._define_empty_table() return log.info("Measured {} sources in {}".format(len(sources), self.image.imgname)) @@ -1069,10 +1068,56 @@ def identify_sources(self, **pars): self.sources = self.tp_sources['aperture']['sources'] # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + def _define_empty_table(self): + """Create basic empty table based on total product table to signify no valid source were detected""" + final_col_format = {"xcentroid": "10.3f", "ycentroid": "10.3f", + "RA": "13.7f", "DEC": "13.7f", + "id": "7d", "Flags": "5d"} + + final_col_descrip = {"xcentroid": "Pixel Coordinate", "ycentroid": "Pixel Coordinate", + "RA": "Sky coordinate at epoch of observation", + "DEC": "Sky coordinate at epoch of observation", + "id": "Catalog Object Identification Number", + "Flags": "Numeric encoding for conditions on detected sources"} + + final_col_units = {"xcentroid": "pixels", "ycentroid": "pixels", "RA": "degrees", "DEC": "degrees", + "id": "unitless", "Flags": "unitless"} + + final_colnames = [k for k in final_col_format.keys()] + + # Initialize empty table with desired column names, descriptions and units + empty_table = Table(names=final_colnames, + descriptions=final_col_descrip, + units=final_col_units) + + # Add formatting for each column + for fcf_key in final_col_format.keys(): + empty_table[fcf_key].format = final_col_format[fcf_key] + + self.sources = empty_table + def measure_sources(self, filter_name): """Perform aperture photometry on identified sources """ + if len(self.sources) == 0: + # Report configuration values to log + log.info("{}".format("=" * 80)) + log.info("") + log.info("No point sources identified for photometry for") + log.info("image name: {}".format(self.imgname)) + log.info("Generating empty point-source catalog.") + log.info("") + # define this attribute for use by the .write method + self.source_cat = self.sources + + self.subset_filter_source_cat = Table(names=["ID", "MagAp2", "CI", "Flags"]) + self.subset_filter_source_cat.rename_column("MagAp2", "MagAP2_" + filter_name) + self.subset_filter_source_cat.rename_column("CI", "CI_" + filter_name) + self.subset_filter_source_cat.rename_column("Flags", "Flags_" + filter_name) + + return + log.info("Performing aperture photometry on identified point-sources") # Open and background subtract image image = self.image.data.copy() @@ -1201,13 +1246,19 @@ def write_catalog(self, reject_catalogs): Nothing! """ - if not reject_catalogs: - # Write out catalog to ecsv file - self.source_cat = self.annotate_table(self.source_cat, self.param_dict_qc, proc_type="aperture", product=self.image.ghd_product) - # self.source_cat.meta['comments'] = \ - # ["NOTE: The X and Y coordinates in this table are 0-indexed (i.e. the origin is (0,0))."] - self.source_cat.write(self.sourcelist_filename, format=self.catalog_format) - log.info("Wrote catalog file '{}' containing {} sources".format(self.sourcelist_filename, len(self.source_cat))) + # Insure catalog has all necessary metadata + self.source_cat = self.annotate_table(self.source_cat, self.param_dict_qc, proc_type="aperture", + product=self.image.ghd_product) + if reject_catalogs: + # We still want to write out empty files + # This will delete all rows from the existing table + self.source_cat.remove_rows(slice(0, None)) + + # Write out catalog to ecsv file + # self.source_cat.meta['comments'] = \ + # ["NOTE: The X and Y coordinates in this table are 0-indexed (i.e. the origin is (0,0))."] + self.source_cat.write(self.sourcelist_filename, format=self.catalog_format) + log.info("Wrote catalog file '{}' containing {} sources".format(self.sourcelist_filename, len(self.source_cat))) # Write out region file if in diagnostic_mode. if self.diagnostic_mode: @@ -1447,6 +1498,7 @@ def identify_sources(self, **pars): check_big_island_only=False, rw2d_biggest_source=self._rw2d_biggest_source, rw2d_source_fraction=self._rw2d_source_fraction) + segm_img_orig = copy.deepcopy(g_segm_img) # If the science field via the segmentation map is deemed crowded or has big sources/islands, compute the # RickerWavelet2DKernel and call detect_and_eval_segments() again. Still use the custom fwhm as it @@ -1599,6 +1651,7 @@ def identify_sources(self, **pars): log.warning("The Round 2 of segmentation images still contain big sources/islands or a\n" "large source fraction of segments.") log.warning("The segmentation algorithm is unable to continue and no segmentation catalog will be produced.") + self._define_empty_table(rw_segm_img) del g_segm_img del rw_segm_img return @@ -1619,6 +1672,7 @@ def identify_sources(self, **pars): # No segments were detected in the total data product - no further processing done for this TDP, # but processing of another TDP should proceed. elif not rw_segm_img: + self._define_empty_table(rw_segm_img) return # The first round custom/Gaussian segmentation image is good, continue with the processing @@ -1630,15 +1684,16 @@ def identify_sources(self, **pars): # No segments were detected in the total data product - no further processing done for this TDP, # but processing of another TDP should proceed. elif not g_segm_img: + self._define_empty_table(g_segm_img) return # Deblend the segmentation image ncount += 1 - self.deblend_segments(segm_img, - imgarr, - ncount, - filter_kernel=self.kernel, - source_box=self._size_source_box) + segm_img = self.deblend_segments(segm_img, + imgarr, + ncount, + filter_kernel=self.kernel, + source_box=self._size_source_box) # The total product catalog consists of at least the X/Y and RA/Dec coordinates for the detected # sources in the total drizzled image. All the actual measurements are done on the filtered drizzled @@ -1745,7 +1800,7 @@ def detect_and_eval_segments(self, imgarr, kernel, ncount, size_source_box, nsig # a final message for this particular total detection product and return. if segm_img is None: log.warning("End processing for the segmentation catalog due to no sources detected with the current kernel.") - log.warning("No segmentation catalog will be produced for this total detection product, {}.".format(self.imgname)) + log.warning("An empty catalog will be produced for this total detection product, {}.".format(self.imgname)) is_big_crowded = True big_island = 1.0 source_fraction = 1.0 @@ -1787,14 +1842,13 @@ def compute_threshold(self, nsigma, bkg_mean, bkg_rms): """ log.info("Computing the threshold value used for source detection.") - if not self.tp_masks: threshold = bkg_mean + (nsigma * bkg_rms) else: threshold = np.zeros_like(self.tp_masks[0]['rel_weight']) log.info("Using WHT masks as a scale on the RMS to compute threshold detection limit.") for wht_mask in self.tp_masks: - threshold_rms = bkg_rms * np.sqrt(wht_mask['mask'] / wht_mask['rel_weight'].max()) + threshold_rms = bkg_rms * np.sqrt(wht_mask['scale'] * wht_mask['mask'] / wht_mask['rel_weight'].max()) threshold_rms_median = np.nanmedian(threshold_rms[threshold_rms > 0]) threshold_item = bkg_mean + (nsigma * threshold_rms_median) threshold += threshold_item @@ -1898,7 +1952,8 @@ def deblend_segments(self, segm_img, imgarr, ncount, filter_kernel=None, source_ log.info("Deblending segments in total image product.") # Note: SExtractor has "connectivity=8" which is the default for detect_sources(). - + # Initialize return value in case of failure in deblending + segm_deblended_img = segm_img try: # Deblending is a combination of multi-thresholding and watershed # segmentation. Sextractor uses a multi-thresholding technique. @@ -1909,21 +1964,21 @@ def deblend_segments(self, segm_img, imgarr, ncount, filter_kernel=None, source_ npixels=source_box, filter_kernel=filter_kernel, nlevels=self._nlevels, - contrast=self._contrast) + contrast=self._contrast, + labels=segm_img.big_segments) if self.diagnostic_mode: + log.info("Deblended {} out of {} segments".format(len(segm_img.big_segments), segm_img.nlabels)) outname = self.imgname.replace(".fits", "_segment_deblended" + str(ncount) + ".fits") fits.PrimaryHDU(data=segm_deblended_img.data).writeto(outname) - # The deblending was successful, so just copy the deblended sources back to the sources attribute. - segm_img = copy.deepcopy(segm_deblended_img) - del segm_deblended_img - except Exception as x_cept: log.warning("Deblending the segments in image {} was not successful: {}.".format(self.imgname, x_cept)) log.warning("Processing can continue with the non-deblended segments, but the user should\n" "check the output catalog for issues.") + # The deblending was successful, so just return the deblended SegmentationImage to calling routine. + return segm_deblended_img # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - def measure_sources(self, filter_name): @@ -1939,6 +1994,24 @@ def measure_sources(self, filter_name): ------- """ + if self.sources is None or self.sources.nlabels == 0: + # Report configuration values to log + log.info("{}".format("=" * 80)) + log.info("") + log.info("No segmentation sources identified for photometry for") + log.info("image name: {}".format(self.imgname)) + log.info("Generating empty segment catalog.") + log.info("") + self._define_empty_table(None) + + # Capture specified filter columns in order to append to the total detection table + self.subset_filter_source_cat = Table(names=["ID", "MagAp2", "CI", "Flags"]) + self.subset_filter_source_cat.rename_column("MagAp2", "MagAP2_" + filter_name) + self.subset_filter_source_cat.rename_column("CI", "CI_" + filter_name) + self.subset_filter_source_cat.rename_column("Flags", "Flags_" + filter_name) + + return + # Get filter-level science data imgarr = copy.deepcopy(self.image.data) @@ -2312,6 +2385,37 @@ def _define_filter_table(self, filter_table): return(final_filter_table) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + def _define_empty_table(self, segm_img): + """Create basic empty table based on total_table format to signify no valid sources were found""" + + final_col_unit = {"X-Centroid": "pixels", "Y-Centroid": "pixels", + "RA": "degrees", "DEC": "degrees", "Flags": "unitless"} + final_col_format = {"ID": "7d", + "X-Centroid": "10.3f", "Y-Centroid": "10.3f", + "RA": "13.7f", "DEC": "13.7f", + "Flags": "5d"} + final_col_descrip = {"ID": "Catalog Object Identification Number", + "X-Centroid": "Pixel Coordinate", + "Y-Centroid": "Pixel Coordinate", + "RA": "Sky coordinate at epoch of observation", + "DEC": "Sky coordinate at epoch of observation", + "Flags": "Numeric encoding for conditions on detected sources"} + + final_colnames = [k for k in final_col_format.keys()] + # Initialize empty table with desired column names, descriptions and units + empty_table = Table(names=final_colnames, + descriptions=final_col_descrip, + units=final_col_unit) + + # Add formatting for each column + for fcf_key in final_col_format.keys(): + empty_table[fcf_key].format = final_col_format[fcf_key] + + self.source_cat = empty_table + self.sources = copy.deepcopy(segm_img) + if self.sources: + self.sources.nlabels = 0 # Insure nlabels is set to 0 to indicate no valid sources + def _define_total_table(self, updated_table): """Set the overall format for the total detection output catalog. @@ -2434,15 +2538,22 @@ def _evaluate_segmentation_image(self, segm_img, image_data, big_island_only=Fal log.info("Segmentation image is blank.") return is_poor_quality, biggest_source, source_fraction + # If the segmentation image is not blank, start out assuming it is good. + is_poor_quality = False + # segm_img is a SegmentationImage, nbins must be at least 1 or segm_img == None nbins = segm_img.max_label log.info("Number of sources from segmentation map: %d", nbins) - # narray = np.bincount(segm_img) - # n = narray[1:] n, binedges = np.histogram(segm_img.data, range=(1, nbins)) real_pixels = (image_data != 0).sum() + # Compute which segments are larger than the kernel. + deb_limit = self.kernel.size + log.debug("Deblending limit set at: {} pixels".format(deb_limit)) + # add as attribute to SegmentationImage for use later + segm_img.big_segments = np.where(segm_img.areas >= deb_limit)[0] + 1 # Segment labels are 1-based + biggest_source = n.max()/float(real_pixels) log.info("Biggest_source: %f", biggest_source) if biggest_source > max_biggest_source: @@ -2479,11 +2590,16 @@ def write_catalog(self, reject_catalogs): Nothing """ - if not reject_catalogs: - # Write out catalog to ecsv file - self.source_cat = self.annotate_table(self.source_cat, self.param_dict_qc, proc_type="segment", product=self.image.ghd_product) - self.source_cat.write(self.sourcelist_filename, format=self.catalog_format) - log.info("Wrote catalog file '{}' containing {} sources".format(self.sourcelist_filename, len(self.source_cat))) + self.source_cat = self.annotate_table(self.source_cat, self.param_dict_qc, proc_type="segment", + product=self.image.ghd_product) + if reject_catalogs: + # We still want to write out empty files + # This will delete all rows from the existing table + self.source_cat.remove_rows(slice(0, None)) + + # Write out catalog to ecsv file + self.source_cat.write(self.sourcelist_filename, format=self.catalog_format) + log.info("Wrote catalog file '{}' containing {} sources".format(self.sourcelist_filename, len(self.source_cat))) # For debugging purposes only, create a "regions" files to use for ds9 overlay of the segm_img. # Create the image regions file here in case there is a failure. This diagnostic portion of the diff --git a/drizzlepac/haputils/deconvolve_utils.py b/drizzlepac/haputils/deconvolve_utils.py index f335588a4..e77622d60 100644 --- a/drizzlepac/haputils/deconvolve_utils.py +++ b/drizzlepac/haputils/deconvolve_utils.py @@ -5,28 +5,20 @@ import sys import shutil import warnings -from distutils.version import LooseVersion import numpy as np import skimage from astropy.io import fits as fits -import photutils -if LooseVersion(photutils.__version__) < '1.1.0': - from photutils.detection.findstars import (_DAOFindProperties, _StarCutout, - _StarFinderKernel, _find_stars) -else: - from photutils.detection._utils import (_StarCutout, _StarFinderKernel, - _find_stars) - from photutils.detection.daofinder import _DAOFindProperties - from photutils.detection import StarFinderBase, find_peaks from photutils.utils import NoDetectionsWarning from stsci.tools import fileutil as fu from . import astrometric_utils as amutils +from ._detection_utils import (_DAOFindProperties, _StarCutout, + _StarFinderKernel, _find_stars) from .. import astrodrizzle # diff --git a/drizzlepac/haputils/hapcut_utils.py b/drizzlepac/haputils/hapcut_utils.py index 0fd4dad86..61ffe6c0a 100644 --- a/drizzlepac/haputils/hapcut_utils.py +++ b/drizzlepac/haputils/hapcut_utils.py @@ -1,4 +1,3 @@ -from astropy import units as u """The module is a high-level interface to astrocut for use with HAP SVM and MVM files.""" from astrocut import fits_cut @@ -32,7 +31,8 @@ def mvm_id_filenames(sky_coord, cutout_size, log_level=logutil.logging.INFO): """ This function retrieves a table of MVM drizzled image filenames with additional information from the archive. The user can then further cull the table to use as - input to obtain a list of files from the archive. + input to obtain a list of files from the archive. This function will return filter-level, + as well as exposure-level products. Parameters ---------- @@ -102,6 +102,8 @@ def mvm_id_filenames(sky_coord, cutout_size, log_level=logutil.logging.INFO): obs_collection="HST") query_table = vstack([acs_query_table, wfc3_query_table]) + del acs_query_table + del wfc3_query_table # Catch the case where no files are found which satisfied the Query if not query_table: @@ -274,7 +276,7 @@ def mvm_retrieve_files(products, archive=False, clobber=False, log_level=logutil return(local_files) -def make_the_cut(input_files, sky_coord, cutout_size, single_outfile=True, output_dir=".", log_level=logutil.logging.INFO, verbose=False): +def make_the_cut(input_files, sky_coord, cutout_size, output_dir=".", log_level=logutil.logging.INFO, verbose=False): """ This function makes the actual cut in the input MVM drizzled filter-level images. As such it is a high-level interface for the `˜astrocut.cutouts.fits_cut` functionality. @@ -296,16 +298,6 @@ def make_the_cut(input_files, sky_coord, cutout_size, single_outfile=True, outpu in ``cutout_size`` are assumed to be in units of arcseconds. `~astropy.units.Quantity` objects must be in angular units. - single_outfile : bool - Default True. If True, return all cutouts in a single fits file with one cutout per extension, - if False return cutouts in individual fits files. Since both the SCI and WHT extensions of the - input files are cutout, individual fits files will contain two image extensions, a SCI followed - by the WHT. Correspondingly, the single fits file will contain a SCI and then a WHT extension - for each input image. If returing a single file the filename will have the form: - ____.fits. If returning multiple files - each will be named: ____.fits. - Output filename of the form: hst_skycell_cutout-p-ra<##>d<####>-dec<##>d<####>_detector_filter.fits - output_dir : str Default value '.'. The directory to save the cutout file(s) to. @@ -319,9 +311,19 @@ def make_the_cut(input_files, sky_coord, cutout_size, single_outfile=True, outpu Returns ------- - response : str or list - If single_outfile is True returns the single output filepath. Otherwise returns a list of all - the output filepaths. + response : list + Returns a list of all the output filenames. + + + Note: For each input file designated for a cutout, there will be a corresponding output file. + Since both the SCI and WHT extensions of the input files are actually cut, individual fits files + will contain two image extensions, a SCI followed by the WHT. Each filter-level output + filename will be of the form: + hst_cutout_skycell-p-ra<##>d<####>-dec<##>d<####>_instrument_detector_filter[_platescale].fits + Each exposure-level filename will be of the form: + hst_cutout_skycell-p-ra<##>d<####>-dec<##>d<####>_instrument_detector_filter[_platescale]-ipppssoo.fits + + where platescale is not present representing the default of "fine" or has the value of "coarse". """ @@ -333,6 +335,7 @@ def make_the_cut(input_files, sky_coord, cutout_size, single_outfile=True, outpu EXTENSION = [1, 2] # SCI and WHT OUTPUT_PREFIX = "hst_cutout_skycell-" MEMORY_ONLY = True # This code will modify the output before it is written. + SINGLE_OUTFILE = False # Making sure we have an array of images if type(input_files) == str: @@ -350,22 +353,24 @@ def make_the_cut(input_files, sky_coord, cutout_size, single_outfile=True, outpu # MULTIPLE FILES: For each file cutout, there is an HDUList comprised of a PHDU and one or more EHDUs. # The out_HDUList is then a list of HDULists. # SINGLE FILES: There is one bare minimum PHDU followed by all of the EHDUs. - out_HDUList = fits_cut(input_files, sky_coord, cutout_size, correct_wcs=CORRECT_WCS, - extension=EXTENSION, single_outfile=single_outfile, cutout_prefix=OUTPUT_PREFIX, - output_dir=".", memory_only=MEMORY_ONLY, verbose=True) + out_HDUList = [] + try: + out_HDUList = fits_cut(input_files, sky_coord, cutout_size, correct_wcs=CORRECT_WCS, + extension=EXTENSION, single_outfile=SINGLE_OUTFILE, cutout_prefix=OUTPUT_PREFIX, + output_dir=".", memory_only=MEMORY_ONLY, verbose=True) + except Exception as x_cept: + log.error("") + log.error("Exception encountered during the cutout process: {}".format(x_cept)) + log.error("No cutout files were created.") - # hst_cutout_skycell-p-ra<##>d<####>-dec<##>d<####>_detector[_filter].fits + # hst_cutout_skycell-p-ra<##>d<####>-dec<##>d<####>_detector_filter[_platescale][-ipppssoo].fits # Get the whole number and fractional components of the RA and Dec ra_whole = int(sky_coord.ra.value) - ra_frac = str(sky_coord.ra.value % 1).split(".")[1][0:6] + ra_frac = str(sky_coord.ra.value % 1).split(".")[1][0:4] dec_whole = abs(int(sky_coord.dec.value)) - dec_frac = str(sky_coord.dec.value % 1).split(".")[1][0:6] + dec_frac = str(sky_coord.dec.value % 1).split(".")[1][0:4] ns = "s" if sky_coord.dec.degree < 0.0 else "n" - if not single_outfile: - log.info("Returning cutouts as multiple FITS files.") - - # ORIG_FLE example is hst_skycell-p1253x05y09_acs_wfc_f658n_all_drc.fits filename_list = [] for HDU in out_HDUList: extlist = HDU[1:] @@ -373,56 +378,78 @@ def make_the_cut(input_files, sky_coord, cutout_size, single_outfile=True, outpu # Update the EXTNAME for all of the EHDUs for index in range(len(extlist)): input_filename = extlist[index].header["ORIG_FLE"] - tokens = extlist[index].header["ORIG_FLE"].split("_") - skycell = tokens[1].split("-")[1][1:5] + tokens = input_filename.split("_") + skycell = tokens[1].split("-")[1] detector = tokens[3] filter = tokens[4] - old_extname = extlist[index].header["O_EXT_NM"].strip().upper() + label_plus = tokens[5] + old_extname= extlist[index].header["O_EXT_NM"].strip().upper() extlist[index].header["EXTNAME"] = old_extname + "_CUTOUT_" + skycell + "_" + \ detector + "_" + filter - # If generating multiple output files, it is efficient to create them now - # NOTE: The multiple output cutout files are input to the CutoutsCombiner, so - # there is some additional keyword manipulation which is not done for the - # single output file. - if not single_outfile: - # SCI extensions are followed by WHT extensions - when the WHT extension - # has been updated, time to write out the file - if old_extname == "WHT": - - # Construct an MVM-style output filename with detector and filter - output_filename = OUTPUT_PREFIX + "p" + skycell + "-ra" + str(ra_whole) + \ - "d" + ra_frac + "-dec" + ns + str(dec_whole) + "d" + \ - dec_frac + "_" + detector + "_" + filter + ".fits" - - cutout_path = os.path.join(output_dir, output_filename) - - log.info("Cutout FITS filename: {}".format(cutout_path)) - - # Retain some keywords written in the PHDU of the cutout file - # by the astrocut software - ra_obj = HDU[0].header["RA_OBJ"] - dec_obj = HDU[0].header["DEC_OBJ"] - - # Replace the minimal primary header written by the astrocut - # software with the primary header from the corresponding input file, - # so we can retain a lot of information from the observation - HDU[0].header = fits.getheader(input_filename) - - # Put the new RA/DEC_OBJ keywords back - HDU[0].header["RA_OBJ"] = (ra_obj, "[deg] right ascension") - HDU[0].header["DEC_OBJ"] = (dec_obj, "[deg] declination") - - # Update PHDU FILENAME keyword with the new filename - HDU[0].header['FILENAME'] = output_filename - - # Populate a PHDU HISTORY keyword with the input filename - HDU[0].header["HISTORY"] = extlist[index].header["ORIG_FLE"] - - output_HDUs = fits.HDUList(HDU) - output_HDUs.writeto(cutout_path, overwrite=True) - - filename_list.append(output_filename) + # Determine if the file is WFC3/IR which has both a "fine" (default) and + # "coarse" platescale. + plate_scale = "_coarse" if label_plus.upper().find("COARSE") != -1 else "" + + # Since the multiple output cutout files can also be input to the CutoutsCombiner, + # there is some additional keyword manipulation done in the header. + # + # SCI extensions are followed by WHT extensions - when the WHT extension + # has been updated, it is time to write out the file. + if old_extname == "WHT": + + # Construct an MVM-style output filename with detector and filter + output_filename = OUTPUT_PREFIX + skycell + "-ra" + str(ra_whole) + \ + "d" + ra_frac + "-dec" + ns + str(dec_whole) + "d" + \ + dec_frac + "_" + detector + "_" + filter + plate_scale + ".fits" + + # Determine if the original file were a filter-level or exposure-level MVM product + # ORIG_FLE filter-level: hst_skycell-p1253x05y09_acs_wfc_f658n_all_drc.fits + # ORIG_FLE filter-level: hst_skycell-p0081x14y15_wfc3_ir_f128n_coarse-all_drz.fits + # ORIG_FLE filter-level: hst_skycell-p0081x14y15_wfc3_ir_f128n_all_drz.fits (fine scale) + # ORIG_FLE exposure-level: hst_skycell-p0081x14y15_wfc3_ir_f128n_coarse-all-ibp505mf_drz.fits + # NOTE: Be careful of the WFC3/IR filenames which can include "coarse". + ef_discriminant = label_plus.split("-")[-1] + if ef_discriminant.upper() != "ALL": + product_type="EXPOSURE" + output_filename = output_filename.replace(".fits", "-" + ef_discriminant + ".fits") + else: + product_type="FILTER" + + # Examples of output cutout filenames: + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_uvis_f275w.fits + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_wfc_f814w-jbp505jg.fits + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_ir_f128n_coarse.fits + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_ir_f128n_coarse-ibp505mf.fits + cutout_path = os.path.join(output_dir, output_filename) + + log.info("Cutout FITS filename: {}".format(cutout_path)) + + # Retain some keywords written in the PHDU of the cutout file + # by the astrocut software + ra_obj = HDU[0].header["RA_OBJ"] + dec_obj = HDU[0].header["DEC_OBJ"] + + # Replace the minimal primary header written by the astrocut + # software with the primary header from the corresponding input file, + # so we can retain a lot of information from the observation + HDU[0].header = fits.getheader(input_filename) + + # Put the new RA/DEC_OBJ keywords back + HDU[0].header["RA_OBJ"] = (ra_obj, "[deg] right ascension") + HDU[0].header["DEC_OBJ"] = (dec_obj, "[deg] declination") + + # Update PHDU FILENAME keyword with the new filename + HDU[0].header['FILENAME'] = output_filename + + # Insert the new keyword, ORIG_FLE, in the PHDU which is the + # *input* filename. This keyword is also in the EHDUs. + HDU[0].header["ORIG_FLE"] = input_filename + + output_HDUs = fits.HDUList(HDU) + output_HDUs.writeto(cutout_path, overwrite=True) + + filename_list.append(output_filename) # Clean up any files left by `˜astrocut.cutouts.fits_cut` try: @@ -437,50 +464,133 @@ def make_the_cut(input_files, sky_coord, cutout_size, single_outfile=True, outpu "Please delete these files to avoid confusion at your earliest convenience:") pprint(cruft_filenames) - # The single output file option was chosen - do not include a detector or - # filter in the output filename. - if single_outfile: - log.info("Returning cutout as single FITS file.") + return filename_list - # Construct an MVM-style output filename WITHOUT detector and filter - single_token = out_HDUList[0][1].header["ORIG_FLE"].split("_") - single_skycell = single_token[1].split("-")[1][1:5] - output_filename = OUTPUT_PREFIX + "p" + single_skycell + "-ra" + str(ra_whole) + \ - "d" + ra_frac + "-dec" + ns + str(dec_whole) + "d" + \ - dec_frac + ".fits" - cutout_path = os.path.join(output_dir, output_filename) +#def mvm_combine(cutout_files, img_combiner=None, output_dir=".", log_level=logutil.logging.INFO): +def mvm_combine(cutout_files, output_dir=".", log_level=logutil.logging.INFO): + """ + This function combines multiple MVM skycell cutout images from the same detector/filter combination + to create a single view of the requested data. All of the functions in this module are designed to + work in conjunction with one another, so the cutout images should be on the user's local disk. This + task is a high-level wrapper for the `˜astrocut.cutout_processing.combine` functionality. - log.info("Cutout FITS filename: {}".format(cutout_path)) - - # Get the primary header from the first input file, so we can retain - # a lot of information from at least the first observation for the PHDU - out_HDUList[0][0].header = fits.getheader(input_files[0]) + Specifically, this routine will combine filter-level cutouts from multiple skycells, all sharing + the same detector and filter. This routine will also combine exposure-level cutouts from + multiple skycells, all sharing the same detector, filter, and ipppssoo. Images which do not + share a detector and filter with any other image will be ignored. Individual exposures from + a single skycell will also be ignored. + + Parameters + ---------- + cutout_files : list + List of fits image cutout filenames where the cutouts are presumed to have been created + with `drizzlepac.haputils.hapcut_utils.make_the_cut`. - # Update PHDU FILENAME keyword with the new filename - out_HDUList[0][0].header["FILENAME"] = output_filename + # img_combiner : func + # The function to be used to combine the images - # Populate a PHDU HISTORY keyword with the names of all the input filenames - out_HDUList[0][0].header["HISTORY"] = input_files + output_dir : str + Default value '.'. The directory to save the cutout file(s) to. - out_HDUList[0].writeto(cutout_path, overwrite=True) - filename_list.append(output_filename) + log_level : int, optional + The desired level of verbosity in the log statements displayed on the screen and written to the + .log file. Default value is 20, or 'info'. - return filename_list + """ + img_combiner = None -def mvm_combine(cutout_files, img_combiner=None, output_dir=".", log_level=logutil.logging.INFO): + # set logging level to user-specified level + log.setLevel(log_level) + + # Make sure the cutout_files are really a list of MULTIPLE filenames + if type(cutout_files) == str or type(cutout_files) == list and len(cutout_files) < 2: + log.error("The 'mvm_combine' function requires a list of MULTIPLE cutout filenames where" \ + " the files were generated by 'make_the_cut'.") + + # Sort the cutout filenames by detector (primary) and filter (secondary) + cutout_files.sort(key = lambda x: (x.split("_")[3], x.split("_")[4])) + + # Report the cutout files submitted for the combination process + log.info("Input cutout files:") + for cf in cutout_files: + log.info("File: {}".format(cf)) + + # Examples of input cutout filenames + # Filter-level + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_uvis_f275w.fits + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_ir_f128n_coarse.fits + # Exposure-level + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_wfc_f814w-jbp505jg.fits + # hst_cutout_skycell-p0081x14y15-ra84d9207-decs69d8516_ir_f128n_coarse-ibp505mf.fits + # + # Combined filter-level files will be generated for each detector/filter combination + # Combined exposure-level files will be generated for each detector/filter combination + # where the ipppssoo is the same + # + # Walk the sorted input list and create filter-level and exposure-level dictionaries + filter_dict = {} + exposure_dict = {} + for cfile in cutout_files: + + # Since the filename could be modified, open the file and read the FILENAME keyword + hdu0 = fits.getheader(cfile, ext=0) + cf = hdu0["FILENAME"].replace(".fits", "") + + # Parse to get the important information + tokens = cf.split("_") + detector = tokens[3] + filter = tokens[4].split("-")[0] + str_tmp = tokens[-1].split("-") + ipppssoo = "" + if len(str_tmp) > 1: + ipppssoo = str_tmp[1] + + # Based upon type of input file, filter-level or exposure-level, populate + # the appropriate dictionary + det_filt_ippp = "" + det_filt = "" + if ipppssoo: + det_filt_ippp = detector + "_" + filter + "_" + ipppssoo + if det_filt_ippp not in exposure_dict: + exposure_dict[det_filt_ippp] = [cfile] + else: + exposure_dict[det_filt_ippp].append(cfile) + else: + det_filt = detector + "_" + filter + if det_filt not in filter_dict: + filter_dict[det_filt] = [cfile] + else: + filter_dict[det_filt].append(cfile) + + # FILTER-LEVEL COMBINATION + # For each detector/filter, generate the output filename and perform the combine + log.info("") + log.info("=== Combining filter-level files ===") + __combine_cutouts(filter_dict, type="FILTER", img_combiner=img_combiner, output_dir=output_dir, log_level=log_level) + + # EXPOSURE-LEVEL COMBINATION + log.info("") + log.info("=== Combining exposure-level files ===") + __combine_cutouts(exposure_dict, type="EXPOSURE", img_combiner=img_combiner, output_dir=output_dir, log_level=log_level) + + log.info("Cutout combination is done.") + + +def __combine_cutouts(input_dict, type="FILTER", img_combiner=None, output_dir=".", log_level=logutil.logging.INFO): """ - This function combines multiple MVM skycell cutout images to create a single view of the - requested data. All of the functions in this module are designed to work in conjunction - with one another, so the cutout images should be on the user's local disk. This task is - a high-level wrapper for the `˜astrocut.cutout_processing.combine` functionality. + This private function performs the actual combine of the multiple MVM skycell cutout images. Parameters ---------- - cutout_files : list - List of fits image cutout filenames where the cutouts are presumed to have been created - with `drizzlepac.haputils.hapcut_utils.make_the_cut`. + input_dict : dictionary + A dictionary where the key is the detector_filter or detector_filter_ipppssoo string and + the corresponding value is a list of filenames corresponding to the key. + + type : string + A string to indicate whether the input_dict variable is for a filter-level or exposure-level + dictionary img_combiner : func The function to be used to combine the images @@ -492,32 +602,50 @@ def mvm_combine(cutout_files, img_combiner=None, output_dir=".", log_level=logut The desired level of verbosity in the log statements displayed on the screen and written to the .log file. Default value is 20, or 'info'. - Note: The combiner only works with cutouts which have been generated as multiple output files, - rather than one large file with many cutout extensions. - """ - # Make sure the cutout_files are really a list of MULTIPLE filenames - if type(cutout_files) == str or type(cutout_files) == list and len(cutout_files) < 2: - log.error("The 'mvm_combine' function requires a list of MULTIPLE cutout filenames where" \ - " the files were generated by 'make_the_cut' with the parameter single_output=False.") + # set logging level to user-specified level + log.setLevel(log_level) - # Report the cutout files to be used for the combination process - log.info("Input cutout files:") - for cf in cutout_files: - log.info("{}".format(cf)) - - # Cutout filenames: hst_cutout_skycell-p0081-ra84d920799-decs69d851699[_uvis_f275w].fits - # Construct the combined filename based on the first file in the list as only RA and Dec are used - output_prefix = "hst_combined_skycells-" - filename = fits.getheader(cutout_files[0], ext=0)['FILENAME'] - tokens = filename.split("_")[2].split("-") - skycell = tokens[1][1:5] - ra = tokens[2] - dec = tokens[3] - output_filename = os.path.join(output_dir, output_prefix + ra + "-" + dec + ".fits") - - # Combine the SCI and then the WHT extensions in the specified files - log.info("Combining the SCI and then the WHT extensions of the input cutout files.") - combined_cutout = astrocut.CutoutsCombiner(cutout_files, img_combiner=img_combiner).combine(output_file=output_filename) - log.info("The combined output filename is {}.".format(output_filename)) + # Output prefix + OUTPUT_PREFIX = "hst_combined_skycells-" + + for key, file_list in input_dict.items(): + + # If there are multiple files to combine, then do it + if len(file_list) > 1: + + # Construct the combined filename based on the first file in the list + # Example: hst_combined_skycells-ra84d9207-decs69d8516_uvis_f275w.fits + filename = fits.getheader(file_list[0], ext=0)['FILENAME'] + fname = filename.replace(".fits", "") + sky_tokens = fname.split("_")[2].split("-") + skycell = sky_tokens[1][1:5] + ra = sky_tokens[2] + dec = sky_tokens[3] + + detector = key.split("_")[0] + filter = key.split("_")[1] + if type.upper() == "EXPOSURE": + exposure = key.split("_")[2] + output_filename = os.path.join(output_dir, OUTPUT_PREFIX + ra + "-" + dec + "_" + \ + detector + "_" + filter + "_" + exposure + ".fits") + else: + output_filename = os.path.join(output_dir, OUTPUT_PREFIX + ra + "-" + dec + "_" + \ + detector + "_" + filter + ".fits") + + # Combine the SCI and then the WHT extensions in the specified files + log.info("Combining the SCI and then the WHT extensions of the input cutout files.") + try: + combined_cutout = astrocut.CutoutsCombiner(file_list, img_combiner=img_combiner).combine(output_file=output_filename) + except Exception as x_cept: + log.warning("The cutout combine was not successful for files, {}, due to {}.".format(file_list, x_cept)) + log.warning("Processing continuuing on next possible set of data.") + continue + log.info("The combined output filename is {}.".format(output_filename)) + + # Only a single file + else: + log.warning("There is only one file for this detector/filter[/ipppssoo] combination, so there" \ + " is nothing to combine.") + log.warning("File {} will be ignored for combination purposes.".format(file_list)) diff --git a/drizzlepac/haputils/make_poller_files.py b/drizzlepac/haputils/make_poller_files.py index d9f675cf4..3bb77d261 100644 --- a/drizzlepac/haputils/make_poller_files.py +++ b/drizzlepac/haputils/make_poller_files.py @@ -1,7 +1,7 @@ #!/usr/bin/env python -"""Makes .out files used as input to runsinglehap.py, runmultihap.py based on the files found in the current -working dir""" +"""Makes .out files used as input to runsinglehap.py, runmultihap.py based on the files or rootnames listed +user-specified list file.""" import argparse import os @@ -13,6 +13,7 @@ __taskname__ = 'make_poller_files' + def generate_poller_file(input_list, poller_file_type='svm', output_poller_filename="poller_file.out", skycell_name=None): """Creates a properly formatted SVM or MVM poller file. @@ -45,7 +46,7 @@ def generate_poller_file(input_list, poller_file_type='svm', output_poller_filen output_list = [] for rootname in rootname_list: rootname = rootname.strip() - fullfilepath = locate_fitspath_from_rootname(rootname) + fullfilepath = locate_fitsfile(rootname) if len(fullfilepath) > 0: print("Rootname {}: Found fits file {}".format(rootname, fullfilepath)) imgname = fullfilepath.split("/")[-1] @@ -60,7 +61,7 @@ def generate_poller_file(input_list, poller_file_type='svm', output_poller_filen imghdu = fits.open(fullfilepath) imghdr = imghdu[0].header linelist.append("{}".format(imghdr['proposid'])) - linelist.append(imgname[1:4].upper()) + linelist.append(imgname.split("_")[-2][1:4].upper()) linelist.append(imghdr['linenum'].split(".")[0]) linelist.append("{}".format(imghdr['exptime'])) if imghdr['INSTRUME'].lower() == "acs": @@ -70,7 +71,10 @@ def generate_poller_file(input_list, poller_file_type='svm', output_poller_filen linelist.append(filter.upper()) linelist.append(imghdr['detector'].upper()) if poller_file_type == 'mvm': # Additional stuff to add to MVM poller files - linelist.append("skycell-{}".format(skycell_name)) + if skycell_name.startswith("skycell-"): + linelist.append("{}".format(skycell_name)) + if skycell_name.startswith("p"): + linelist.append("skycell-{}".format(skycell_name)) linelist.append("NEW") linelist.append(fullfilepath) imghdu.close() @@ -93,33 +97,51 @@ def generate_poller_file(input_list, poller_file_type='svm', output_poller_filen # ============================================================================================================ -def locate_fitspath_from_rootname(rootname): - """returns full file name (fullpath + filename) for a specified rootname. +def locate_fitsfile(search_string): + """returns full file name (fullpath + filename) for a specified rootname or filename. The search + algorithm looks for the file in the following order: + + - Search for a _flc.fits file in the current working directory + - Search for a _flt.fits file in the current working directory + - Search for a _flc.fits file in subdirectory in the path specified in $DATA_PATH + - Search for a _flt.fits file in subdirectory in the path specified in $DATA_PATH Parameters ---------- - rootname : str - rootname to process + search_string : str + rootname or filename to locate Returns ------- fullfilepath : str - full path + image name of specified rootname. + full file path + image name of specified search_string. """ - # Look for files in CWD first - if os.path.exists("{}_flc.fits".format(rootname)): - return "{}_flc.fits".format(rootname) - if os.path.exists("{}_flt.fits".format(rootname)): - return "{}_flt.fits".format(rootname) - # If not found in CWD, look elsewhere... - if not os.getenv("DATA_PATH"): - sys.exit("ERROR: Undefined online cache data root path. Please set environment variable 'DATA_PATH'") - filenamestub = "{}/{}/{}/{}".format(os.getenv("DATA_PATH"), rootname[:4], rootname, rootname) - if os.path.exists("{}_flc.fits".format(filenamestub)): - fullfilepath = "{}_flc.fits".format(filenamestub) - else: - fullfilepath = "{}_flt.fits".format(filenamestub) - return fullfilepath + if search_string.endswith("_flt.fits") or search_string.endswith("_flc.fits"): # Process search_string as a full filename + # Look for files in CWD first + if os.path.exists(search_string): + return os.getcwd()+"/"+search_string + # If not found in CWD, look elsewhere... + if not os.getenv("DATA_PATH"): + sys.exit("ERROR: Undefined online cache data root path. Please set environment variable 'DATA_PATH'") + fullfilepath = "{}/{}/{}/{}".format(os.getenv("DATA_PATH"), search_string[:4], search_string, search_string) + if os.path.exists(search_string): + return fullfilepath + + else: # Process search_string as a rootname + # Look for files in CWD first + if os.path.exists("{}_flc.fits".format(search_string)): + return "{}/{}_flc.fits".format(os.getcwd(), search_string) + if os.path.exists("{}_flt.fits".format(search_string)): + return "{}/{}_flt.fits".format(os.getcwd(), search_string) + # If not found in CWD, look elsewhere... + if not os.getenv("DATA_PATH"): + sys.exit("ERROR: Undefined online cache data root path. Please set environment variable 'DATA_PATH'") + filenamestub = "{}/{}/{}/{}".format(os.getenv("DATA_PATH"), search_string[:4], search_string, search_string) + if os.path.exists("{}_flc.fits".format(filenamestub)): + fullfilepath = "{}_flc.fits".format(filenamestub) + else: + fullfilepath = "{}_flt.fits".format(filenamestub) + return fullfilepath # ============================================================================================================ @@ -129,7 +151,8 @@ def locate_fitspath_from_rootname(rootname): parser = argparse.ArgumentParser(description='Create a HAP SVM or MVM poller file') parser.add_argument('input_list', - help='Name of a file containing a list of rootnames (9 characters, usually ending ' + help='Name of a file containing a list of calibrated fits files (ending with ' + '"_flt.fits" or "_flc.fits") or rootnames (9 characters, usually ending ' 'with a "q" to process. The corresponding flc.fits or flt.fits files must ' 'exist in the online cache') parser.add_argument('-o', '--output_poller_filename', required=False, default="poller_file.out", diff --git a/drizzlepac/haputils/poller_utils.py b/drizzlepac/haputils/poller_utils.py index 1aa91d162..7236d99a8 100755 --- a/drizzlepac/haputils/poller_utils.py +++ b/drizzlepac/haputils/poller_utils.py @@ -283,7 +283,7 @@ def interpret_mvm_input(results, log_level, layer_method='all', exp_limit=2.0, u obset_table = copy.deepcopy(user_table) # Add INSTRUMENT column - instr = [INSTRUMENT_DICT[fname[0]] for fname in obset_table['filename']] + instr = [INSTRUMENT_DICT[fname.split("_")[-2][0]] for fname in obset_table['filename']] # convert input to an Astropy Table for parsing obset_table.add_column(Column(instr, name='instrument')) @@ -457,7 +457,7 @@ def parse_mvm_tree(det_tree, all_mvm_exposures, log_level): # Parse the first filename[1] to determine if the products are flt or flc if det_indx != prev_det_indx: filetype = "drc" - if filename[1][10:13].lower().endswith("flt"): + if filename[1][-8:-5].lower().endswith("flt"): filetype = "drz" prev_det_indx = det_indx @@ -651,7 +651,7 @@ def parse_obset_tree(det_tree, log_level): # easier for the user by having the same WCS in both the direct and # Grism/Prism products. # - # The GrismExposureProduct is only an attibutes of the TotalProduct. + # The GrismExposureProduct is only an attribute of the TotalProduct. prod_list = prod_info.split(" ") # Determine if this image is a Grism/Prism or a nominal direct exposure @@ -659,8 +659,9 @@ def parse_obset_tree(det_tree, log_level): if prod_list[5].lower().find('g') != -1 or prod_list[5].lower().find('pr') != -1: is_grism = True filt_indx -= 1 - grism_sep_obj = GrismExposureProduct(prod_list[0], prod_list[1], prod_list[2], prod_list[3], - filename[1], prod_list[5], prod_list[6], log_level) + grism_sep_obj = GrismExposureProduct(prod_list[0], prod_list[1], prod_list[2], + prod_list[3], filename[1], prod_list[5], + prod_list[6], log_level) else: sep_obj = ExposureProduct(prod_list[0], prod_list[1], prod_list[2], prod_list[3], filename[1], prod_list[5], prod_list[6], log_level) @@ -1199,7 +1200,8 @@ def add_primary_fits_header_as_attr(hap_obj, log_level=logutil.logging.NOTSET): object to update log_level : int, optional. - The desired level of verboseness in the log statements displayed on the screen and written to the .log file. + The desired level of verboseness in the log statements displayed on the screen and written to the + .log file. Returns ------- diff --git a/drizzlepac/haputils/product.py b/drizzlepac/haputils/product.py index 6a3057dc2..335f41e2f 100755 --- a/drizzlepac/haputils/product.py +++ b/drizzlepac/haputils/product.py @@ -184,6 +184,14 @@ def align_to_gaia(self, catalog_list=[], output=True, align_table.find_alignment_sources(output=output, crclean=crclean) is_good_fit = False + # determine minimum number of sources available for alignment + source_nums = [] + for img_cats in align_table.extracted_sources.values(): + sources = 0 + for chip in img_cats: + sources += len(img_cats[chip]) + source_nums.append(sources) + min_sources = min(source_nums) # Loop for available catlogs, the first successful fit for a # (catalog, fitting method, and fit geometry) is satisfactory to break out of the looping. @@ -206,7 +214,8 @@ def align_to_gaia(self, catalog_list=[], output=True, full_catalog=True) # Add a weight column which is based upon proper motion measurements - if ref_catalog.meta['converted'] and ref_catalog['RA_error'][0] != np.nan: + if 'converted' in ref_catalog.meta and ref_catalog.meta['converted'] \ + and ref_catalog['RA_error'][0] != np.nan: ref_weight = np.sqrt(ref_catalog['RA_error'] ** 2 + ref_catalog['DEC_error'] ** 2) ref_weight = np.nan_to_num(ref_weight, copy=False, nan=0.0, posinf=0.0, neginf=0.0) else: @@ -242,7 +251,7 @@ def align_to_gaia(self, catalog_list=[], output=True, # If there are not enough references sources for the specified fitgeom, # downgrade the fitgeom until a valid one is found. Also, if the fit done - # with the fitgeom was unsatisfactory, downgrade if possible and try again. + # with the fitgeom was unsatisfactory, downgrade if possible and try again while num_ref_sources < mosaic_fitgeom_list[mosaic_fitgeom_index][1]: log.warning("Not enough reference sources for alignment using catalog '{}' with fit method '{}' and fit geometry '{}'.".format(catalog_item, method_name, mosaic_fitgeom_list[mosaic_fitgeom_index][0])) mosaic_fitgeom_index -= 1 @@ -251,7 +260,6 @@ def align_to_gaia(self, catalog_list=[], output=True, break while mosaic_fitgeom_index > -1: - # In case of a downgrade, make sure the fit geometry "name" is based upon the index mosaic_fitgeom = mosaic_fitgeom_list[mosaic_fitgeom_index][0] log.info("Proceeding to use the '{}' fit geometry.".format(mosaic_fitgeom)) @@ -305,8 +313,19 @@ def align_to_gaia(self, catalog_list=[], output=True, format(catalog_item, method_name, mosaic_fitgeom)) traceback.print_exc() + # Try again with a different fit geometry algorithm mosaic_fitgeom_index -= 1 + log.info("Resetting mosaic index for next fit method") + # Only if there are too few sources to perform an alignment with the current + # method should the next method be attempted. Otherwise, it is assumed that the + # image sources are either cosmic-rays or features that do not match across + # the filters (e.g., F300X vs F475W from 'id7r03') and can't be aligned. Trying with a + # less strict method (rshift vs rscale) will only increase the probability of + # matching incorrectly and introducing an error to the relative alignment. + if min_sources > mosaic_fitgeom_list[mosaic_fitgeom_index][1]: + mosaic_fitgeom_index = -1 + log.warning("Too many (bad?) sources to try any more fitting with this catalog.") # Break out of the "fit method" for loop if is_good_fit: @@ -794,10 +813,14 @@ def __init__(self, prop_id, obset_id, instrument, detector, filename, layer, fil # layer: [filter_str, pscale_str, exptime_str, epoch_str] # e.g.: [f160w, coarse, all, all] # + if filename.startswith("hst"): + exposure_name = filename.split("_")[-2] + else: + exposure_name = self.exposure_name if layer[1] == 'coarse': - layer_vals = [layer[1], layer[2], self.exposure_name] + layer_vals = [layer[1], layer[2], exposure_name] else: - layer_vals = ['all', self.exposure_name] + layer_vals = ['all', exposure_name] layer_str = '-'.join(layer_vals) cell_id = "p{}{}".format(prop_id, obset_id) @@ -906,8 +929,11 @@ def copy_exposure(self, filename): New MVM-compatible HAP filename for input exposure """ - suffix = filename.split("_")[1] - sce_filename = '_'.join([self.product_basename, suffix]) + if filename.startswith("hst"): + sce_filename = "{}_{}".format(self.product_basename, filename.split("_")[-1]) + else: + suffix = filename.split("_")[1] + sce_filename = '_'.join([self.product_basename, suffix]) log.info("Copying {} to MVM input: \n {}".format(filename, sce_filename)) try: shutil.copy(filename, sce_filename) @@ -1079,6 +1105,7 @@ def wcs_drizzle_product(self, meta_wcs): .format(meta_wcs, self.trl_logname)) edp_filenames = [element.full_filename for element in self.edp_list] + astrodrizzle.AstroDrizzle(input=edp_filenames, output=self.drizzle_filename, **drizzle_pars) diff --git a/drizzlepac/haputils/which_skycell.py b/drizzlepac/haputils/which_skycell.py new file mode 100644 index 000000000..12e24bed9 --- /dev/null +++ b/drizzlepac/haputils/which_skycell.py @@ -0,0 +1,61 @@ +#!/usr/bin/env python +"""Reports the name(s) of the projection cell(s)/skycell(s) that the observations in the user-specified input + file occupy.""" + +import argparse +import glob +import os +import sys + +from drizzlepac.haputils import cell_utils + + +__version__ = 0.1 +__version_date__ = '08-Sept-2021' +# ------------------------------------------------------------------------------------------------------------ + + +def report_skycells(img_list): + """reports the names of the projection cell(s)/skycell(s) that the observations in the user-specified + input list occupy. + + Parameters + ---------- + img_list : list + list of the images to process + + Returns + ------- + Nothing! + """ + skycell_dict = cell_utils.get_sky_cells(img_list) + print("\n") + for skycell_name in skycell_dict.keys(): + print("Skycell {} contains {} image(s)".format(skycell_name, len(skycell_dict[skycell_name].members))) + for imgname in skycell_dict[skycell_name].members: + print(" {}".format(imgname)) + print("\n") + +# ------------------------------------------------------------------------------------------------------------ + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='identifies the projection cell(s) and skycell(s) that the ' + 'user-specified observations occupy') + parser.add_argument('-i', '--input_list', required=False, default="NONE", + help='Name of a file containing a list of calibrated fits files (ending with ' + '"_flt.fits" or "_flc.fits") to process. If not explicitly specified, all ' + 'flc/flt fits files in the current working directory will be processed.') + in_args = parser.parse_args() + if in_args.input_list is "NONE": + img_list = glob.glob("*_fl?.fits") + if len(img_list) > 0: + report_skycells(img_list) + else: + sys.exit("ERROR: No flc/flt fits files found in current path {}/".format(os.getcwd())) + elif os.path.exists(in_args.input_list): + with open(in_args.input_list) as fin: + img_list = fin.read().splitlines() + report_skycells(img_list) + else: + sys.exit("ERROR: Input file {} does not exist!".format(in_args.input_list)) diff --git a/drizzlepac/make_custom_mosaic.py b/drizzlepac/make_custom_mosaic.py index 345da27fc..ac0d5d1ac 100644 --- a/drizzlepac/make_custom_mosaic.py +++ b/drizzlepac/make_custom_mosaic.py @@ -177,7 +177,7 @@ def create_poller_file(img_list, proj_cell_dict): poller_filename = "temp_{}_mvm.out".format(skycell_name) rootname_list = [] for imgname in img_list: - rootname_list.append(imgname.split("_")[0]+"\n") + rootname_list.append(imgname+"\n") tf = tempfile.NamedTemporaryFile(mode='w+t', dir=os.getcwd()) with open(tf.name, 'w') as f: f.writelines(rootname_list) @@ -313,7 +313,7 @@ def determine_projection_cell(img_list): # ------------------------------------------------------------------------------------------------------------ -def perform(input_image_source, log_level='info', output_file_prefix=None): +def perform(input_image_source, log_level='info', output_file_prefix=None, skip_gaia_alignment=False): """Main calling subroutine Parameters @@ -340,6 +340,10 @@ def perform(input_image_source, log_level='info', output_file_prefix=None): skycell = 1974, ra = 201.9512, and dec = +26.0012, The filename prefix would be "skycell-p1974-ra201d9512-decn26d0012". + skip_gaia_alignment : bool, optional + Skip alignment of all input images to known Gaia/HSC sources in the input image footprint? If set to + 'True', The existing input image alignment solution will be used instead. The default is False. + Returns ------- return_value : int @@ -384,7 +388,8 @@ def perform(input_image_source, log_level='info', output_file_prefix=None): return_value = hapmultisequencer.run_mvm_processing(poller_filename, custom_limits=custom_limits, log_level=logging_level, - output_file_prefix=output_file_prefix) + output_file_prefix=output_file_prefix, + skip_gaia_alignment=skip_gaia_alignment) except Exception: if return_value == 0: @@ -443,8 +448,15 @@ def main(): 'Note that the "" denotes if the declination is north (positive) or south ' '(negative). Example: For skycell = 1974, ra = 201.9512, and dec = +26.0012, ' 'The filename prefix would be "skycell-p1974-ra201d9512-decn26d0012".') + parser.add_argument('-s', '--skip_gaia_alignment', required=False, action='store_true', + help='Skip alignment of all input images to known Gaia/HSC sources in the input ' + 'image footprint? If this option is turned on, the existing input image ' + 'alignment solution will be used instead. The default is False.') user_args = parser.parse_args() - return_value = perform(user_args.input_image_source, user_args.log_level, user_args.output_file_prefix) + return_value = perform(user_args.input_image_source, + log_level=user_args.log_level, + output_file_prefix=user_args.output_file_prefix, + skip_gaia_alignment=user_args.skip_gaia_alignment) # ------------------------------------------------------------------------------------------------------------ diff --git a/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_alignment_all.json b/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_alignment_all.json index ffa4ff39b..cf6ed6311 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_alignment_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_alignment_all.json @@ -22,7 +22,7 @@ { "box_size":27, "win_size":3, - "nsigma":5.0, + "nsigma":3.0, "centering_mode": "starfind", "bkg_estimator": "MedianBackground", "rms_estimator": "StdBackgroundRMS", diff --git a/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_catalog_generation_all.json b/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_catalog_generation_all.json index 31490a003..ba066940a 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_catalog_generation_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/acs/hrc/acs_hrc_catalog_generation_all.json @@ -32,8 +32,8 @@ "fwhm": 0.073, "source_box": 5, "segm_nsigma": 3.0, - "nlevels": 32, - "contrast": 0.005, + "nlevels": 64, + "contrast": 0.001, "border": 10, "rw2d_size": 15, "rw2d_nsigma": 3.0, diff --git a/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_alignment_all.json b/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_alignment_all.json index ffa4ff39b..492d8be07 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_alignment_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_alignment_all.json @@ -22,14 +22,14 @@ { "box_size":27, "win_size":3, - "nsigma":5.0, + "nsigma":3.0, "centering_mode": "starfind", "bkg_estimator": "MedianBackground", "rms_estimator": "StdBackgroundRMS", "num_sources": 250, "deblend": false, "fwhmpsf": 0.13, - "classify": true, + "classify": false, "threshold": -1.1 }, "generate_astrometric_catalog": diff --git a/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_catalog_generation_all.json b/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_catalog_generation_all.json index 8910f1f5e..bcfe3f2f3 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_catalog_generation_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/acs/sbc/acs_sbc_catalog_generation_all.json @@ -32,8 +32,8 @@ "fwhm": 0.065, "source_box": 5, "segm_nsigma": 3.0, - "nlevels": 32, - "contrast": 0.005, + "nlevels": 64, + "contrast": 0.001, "border": 10, "rw2d_size": 15, "rw2d_nsigma": 3.0, diff --git a/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_alignment_all.json b/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_alignment_all.json index 0f0d20f2c..86653a427 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_alignment_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_alignment_all.json @@ -22,7 +22,7 @@ { "box_size":27, "win_size":3, - "nsigma":5.0, + "nsigma":3.0, "centering_mode": "starfind", "bkg_estimator": "SExtractorBackground", "rms_estimator": "StdBackgroundRMS", diff --git a/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_catalog_generation_all.json b/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_catalog_generation_all.json index 90a5578b2..62576b8e7 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_catalog_generation_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/acs/wfc/acs_wfc_catalog_generation_all.json @@ -32,8 +32,8 @@ "fwhm": 0.13, "source_box": 5, "segm_nsigma": 3.0, - "nlevels": 32, - "contrast": 0.005, + "nlevels": 64, + "contrast": 0.001, "border": 10, "rw2d_size": 15, "rw2d_nsigma": 3.0, diff --git a/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_alignment_all.json b/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_alignment_all.json index dd4bcbbde..a16c73eec 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_alignment_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_alignment_all.json @@ -22,7 +22,7 @@ { "box_size":27, "win_size":3, - "nsigma":5.0, + "nsigma":3.0, "centering_mode": "starfind", "bkg_estimator": "MedianBackground", "rms_estimator": "StdBackgroundRMS", diff --git a/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_catalog_generation_all.json b/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_catalog_generation_all.json index 64345aac0..e4d2d90fb 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_catalog_generation_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/wfc3/ir/wfc3_ir_catalog_generation_all.json @@ -32,7 +32,7 @@ "fwhm": 0.14, "source_box": 5, "segm_nsigma": 3.0, - "nlevels": 32, + "nlevels": 64, "contrast": 0.001, "border": 10, "rw2d_size": 11, diff --git a/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_alignment_all.json b/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_alignment_all.json index bbb5d9495..1aeb8897a 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_alignment_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_alignment_all.json @@ -22,7 +22,7 @@ { "box_size":27, "win_size":3, - "nsigma":5.0, + "nsigma":3.0, "centering_mode": "starfind", "bkg_estimator": "SExtractorBackground", "rms_estimator": "StdBackgroundRMS", diff --git a/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_catalog_generation_all.json b/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_catalog_generation_all.json index 0123a3e01..1dbc852dc 100644 --- a/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_catalog_generation_all.json +++ b/drizzlepac/pars/hap_pars/default_parameters/wfc3/uvis/wfc3_uvis_catalog_generation_all.json @@ -32,8 +32,8 @@ "fwhm": 0.076, "source_box": 5, "segm_nsigma": 3.0, - "nlevels": 32, - "contrast": 0.005, + "nlevels": 64, + "contrast": 0.001, "border": 10, "rw2d_size": 15, "rw2d_nsigma": 3.0, diff --git a/drizzlepac/runmultihap.py b/drizzlepac/runmultihap.py index 9a1e31d73..b5d7b5ad0 100755 --- a/drizzlepac/runmultihap.py +++ b/drizzlepac/runmultihap.py @@ -104,12 +104,17 @@ def main(): 'Specifying "critical" will only record/display "critical" log statements, and ' 'specifying "error" will record/display both "error" and "critical" log ' 'statements, and so on.') + parser.add_argument('-s', '--skip_gaia_alignment', required=False, action='store_true', + help='Skip alignment of all input images to known Gaia/HSC sources in the input ' + 'image footprint? If this option is turned on, the existing input image ' + 'alignment solution will be used instead. The default is False.') user_args = parser.parse_args() print("Multi-visit processing started for: {}".format(user_args.input_filename)) rv = perform(user_args.input_filename, diagnostic_mode=user_args.diagnostic_mode, - log_level=user_args.log_level) + log_level=user_args.log_level, + skip_gaia_alignment=user_args.skip_gaia_alignment) print("Return Value: ", rv) return rv