diff --git a/astrophot/cli/cli_segmap.py b/astrophot/cli/cli_segmap.py new file mode 100644 index 00000000..019ac805 --- /dev/null +++ b/astrophot/cli/cli_segmap.py @@ -0,0 +1,338 @@ +# ============================================================================= +# Fit all objects identified in a segmentation map +# +# This is a quick script to fit all the objects identified in a segmentation map +# using a single model type. The script will load the target image, mask, +# psf, and variance image (if available) and fit the models to the target image. +# +# Run this script with: +# ~$ python segmap_astrophot_model.py target.fits segmap.fits [OPTIONS] +# ============================================================================= + +import astrophot as ap +import numpy as np +from astropy.io import fits +from astropy.wcs import WCS +import argparse + +try: + import yaml +except ImportError: + from astropy.io.misc import yaml + + +def to_serializable(value): + if value is None: + return None + if hasattr(value, "detach"): + value = value.detach().cpu().numpy() + value = np.array(value) + return value.item() if value.shape == () else value.tolist() + + +def collect_model_parameters(model): + params = model.dynamic_params + names = tuple(param.name for param in params) + values = tuple(to_serializable(param.npvalue) for param in params) + uncertainties = tuple(to_serializable(param.uncertainty) for param in params) + return { + name: {"value": val, "uncertainty": unc} + for name, val, unc in zip(names, values, uncertainties) + } + + +def main(): + parser = argparse.ArgumentParser( + description="Fit a model to a series of targets in an image using AstroPhot." + ) + + # Core Parameters + parser.add_argument("target_file", type=str, help="Path to the target FITS file") + parser.add_argument("segmap_file", type=str, help="Path to the segmentation map FITS file") + parser.add_argument( + "--cat", + type=str, + default=None, + help="Path to a catalogue yaml file with initial parameters for each segmentation map id", + ) + parser.add_argument( + "--name", type=str, default="astrophot_model", help="Prefix name used for models" + ) + parser.add_argument("--psf", type=str, default=None, help="Path to the PSF FITS file") + parser.add_argument( + "--psf_upsample", type=int, default=1, help="PSF upsampling factor for convolution (int)" + ) + parser.add_argument("--zeropoint", type=float, default=None, help="Magnitude zeropoint") + parser.add_argument( + "--initial_sky", + type=float, + default=None, + help="Initial sky value for the I0 parameter", + ) + parser.add_argument( + "--sky_locked", action="store_true", help="Lock the sky model during fitting" + ) + parser.add_argument( + "--model_type", + type=str, + default="sersic_galaxy_model", + help="Type of AstroPhot model to fit. Replace spaces with underscores, e.g. 'sersic_galaxy_model' or 'exponential_disk_model'", + ) + parser.add_argument("--verbose", type=int, default=1, help="Verbosity level for fitting output") + parser.add_argument( + "--dump", + action="store_true", + help="Dump this file to the current directory as 'single_astrophot_model.py' for editing and running", + ) + + # Window Parameters + parser.add_argument( + "--window_expand_scale", + type=float, + default=1.0, + help="Scale factor to expand windows from segmentation map for final fit", + ) + parser.add_argument( + "--window_expand_border", + type=int, + default=0, + help="Number of pixels to expand windows from segmentation map for final fit", + ) + parser.add_argument( + "--window_min_size", + type=int, + default=1, + help="Minimum size of windows to include in the fit. Before expanding. Number of pixels.", + ) + parser.add_argument( + "--filter_ids", + type=int, + nargs="+", + default=[0], + help="List of segmentation map ids to exclude from fitting. Default is [0] to exclude the background", + ) + + # Variance and Mask Parameters + parser.add_argument("--variance", type=str, default=None, help="Path to the variance FITS file") + parser.add_argument( + "--variance_hdu", type=int, default=0, help="FITS file index for variance data" + ) + parser.add_argument("--mask", type=str, default=None, help="Path to the mask FITS file") + parser.add_argument("--mask_hdu", type=int, default=0, help="FITS file index for mask data") + + # Extra Parameters + parser.add_argument( + "--no_save_images", + action="store_false", + dest="save_images", + help="Disable saving the model and residual images", + ) + parser.add_argument( + "--no_save_cov", + action="store_false", + dest="save_covariance_matrix", + help="Disable saving the covariance matrix", + ) + parser.add_argument( + "--target_hdu", type=int, default=0, help="FITS file index for target image data" + ) + parser.add_argument( + "--segmap_hdu", type=int, default=0, help="FITS file index for segmentation map data" + ) + parser.add_argument("--psf_hdu", type=int, default=0, help="FITS file index for PSF data") + parser.add_argument( + "--sky_model_type", + type=str, + default="flat", + help="Type of sky model to fit, options include 'flat' or 'plane'", + ) + + # Parse known arguments, leave the rest for the dynamic dictionary + args = parser.parse_args() + + if args.dump: + with open("segmap_astrophot_model.py", "w", encoding="utf-8") as f: + f.write( + "# This file was generated by dumping the segmap_model_cli.py script. You can edit and run this file directly.\n\n" + ) + with open(__file__, "r", encoding="utf-8") as original: + f.write(original.read()) + print("Dumped to segmap_astrophot_model.py") + + # Load Target Data + # --------------------------------------------------------------------- + if args.verbose > 0: + print("Loading target image...") + with fits.open(args.target_file) as hdu: + target_data = np.array(hdu[args.target_hdu].data, dtype=np.float64) + + # Load Variance Data + # --------------------------------------------------------------------- + variance_data = None + if args.variance is not None: + if args.verbose > 0: + print("Loading variance image...") + with fits.open(args.variance) as hdu: + variance_data = np.array(hdu[args.variance_hdu].data, dtype=np.float64) + + # Load Mask Data + # --------------------------------------------------------------------- + mask_data = None + if args.mask is not None: + if args.verbose > 0: + print("Loading mask image...") + with fits.open(args.mask) as hdu: + mask_data = np.array(hdu[args.mask_hdu].data) + + # Load PSF + # --------------------------------------------------------------------- + psf_data = None + if args.psf is not None: + if args.verbose > 0: + print("Loading PSF...") + with fits.open(args.psf) as hdu: + psf_data = np.array(hdu[args.psf_hdu].data, dtype=np.float64) + psf_data = ap.PSFImage(data=psf_data, upsample=args.psf_upsample) + + # Make Target + # --------------------------------------------------------------------- + target_wcs = WCS(fits.getheader(args.target_file, args.target_hdu)) + target = ap.TargetImage( + data=target_data, + wcs=target_wcs, + zeropoint=args.zeropoint, + variance=variance_data, + mask=mask_data, + psf=psf_data, + ) + + # Load Segmentation Map + # --------------------------------------------------------------------- + if args.verbose > 0: + print("Loading segmentation map...") + with fits.open(args.segmap_file) as hdu: + segmap_data = np.array(hdu[args.segmap_hdu].data, dtype=np.int32) + + # Load Catalogue of Initial Parameters + # --------------------------------------------------------------------- + override_init_params = {} + if args.cat is not None: + if args.verbose > 0: + print("Loading catalogue of initial parameters...") + with open(args.cat, "r", encoding="utf-8") as cat_file: + override_init_params = yaml.safe_load(cat_file) + + # Initialization from segmap + # --------------------------------------------------------------------- + if args.verbose > 0: + print("Parsing segmentation map") + windows = ap.utils.initialize.windows_from_segmentation_map( + segmap_data, skip_index=args.filter_ids + ) + + windows = ap.utils.initialize.filter_windows( + windows, + min_area=args.window_min_size, + image=target, + ) + windows = ap.utils.initialize.scale_windows( + windows, + image=target, + expand_scale=args.window_expand_scale, + expand_border=args.window_expand_border, + ) + + centers = ap.utils.initialize.centroids_from_segmentation_map(segmap_data, target) + if "galaxy" in args.model_type: + PAs = ap.utils.initialize.PA_from_segmentation_map(segmap_data, target, centers) + qs = ap.utils.initialize.q_from_segmentation_map(segmap_data, target, centers) + else: + PAs = None + qs = None + init_params = {} + for window in windows: + init_params[window] = {"center": centers[window]} + if "galaxy" in args.model_type: + init_params[window]["PA"] = PAs[window] + init_params[window]["q"] = qs[window] + if window in override_init_params: + init_params[window].update(override_init_params[window]) + + # Create and fit Models + # --------------------------------------------------------------------- + results = {} + model_images = {} + residual_images = {} + cov_matrices = {} + for window in windows: + if args.verbose > 0: + print(f"Fitting model for window {window}...") + W = ap.Window(windows[window], target) + subtarget = target[W] + sub_segmap = segmap_data[target.get_indices(W)] + # Mask interlopers + subtarget.mask = np.logical_or( + subtarget.mask, ~(np.logical_or(sub_segmap == 0, sub_segmap == window)).T + ) + # Sky model + sky = ap.Model( + name="sky", + model_type=args.sky_model_type + " sky model", + target=subtarget, + I0=args.initial_sky, + ) + if args.sky_locked: + sky.to_static() + # Primary object model + object_model = ap.Model( + name=f"{args.model_type}_{window}", + model_type=args.model_type.replace("_", " "), + target=subtarget, + **init_params[window], + ) + # Combine into group model for fitting + model = ap.Model( + name=f"{args.name}_{window}", + model_type="group model", + target=subtarget, + models=[sky, object_model], + ) + # Fitting + model.initialize() + result = ap.fit.LM(model, verbose=args.verbose - 1).fit() + # Collect results + results[model.name] = { + "parameters": collect_model_parameters(object_model), + "total_flux": to_serializable(object_model.total_flux()), + "total_flux_uncertainty": to_serializable(object_model.total_flux_uncertainty()), + "sky_parameters": collect_model_parameters(sky), + "message": result.message, + } + if args.save_covariance_matrix: + cov_matrices[model.name] = to_serializable(result.covariance_matrix) + if args.save_images: + model_images[model.name] = to_serializable(model().data) + residual_images[model.name] = to_serializable((subtarget - model()).data) + if args.zeropoint is not None: + results[model.name]["total_magnitude"] = to_serializable(object_model.total_magnitude()) + results[model.name]["total_magnitude_uncertainty"] = to_serializable( + object_model.total_magnitude_uncertainty() + ) + + # Save outputs + # --------------------------------------------------------------------- + if args.save_covariance_matrix: + np.savez(f"{args.name}_covariance_matrix.npz", **cov_matrices) + + if args.save_images: + hdul = fits.HDUList() + for name, image in model_images.items(): + hdul.append(fits.ImageHDU(data=image, name=f"{name}_model")) + hdul.writeto(f"{args.name}_model_images.fits", overwrite=True) + hdul = fits.HDUList() + for name, image in residual_images.items(): + hdul.append(fits.ImageHDU(data=image, name=f"{name}_residual")) + hdul.writeto(f"{args.name}_residual_images.fits", overwrite=True) + + with open(f"{args.name}_parameters.yaml", "w", encoding="utf-8") as output_file: + yaml.dump(results, output_file, default_flow_style=False) diff --git a/astrophot/cli/cli_single.py b/astrophot/cli/cli_single.py new file mode 100644 index 00000000..e84dae84 --- /dev/null +++ b/astrophot/cli/cli_single.py @@ -0,0 +1,310 @@ +# ============================================================================= +# Fit a single model to a target image +# +# This is a quick script to fit a single AstroPhot model to a target image. The +# script will load the target image, mask, psf, and variance image (if +# available) and fit the model to the target image. The script will save the +# model image, residual image, and covariance matrix to the current directory. +# This script is intended for quick easy fits, and as a starting point to build +# a more complex analysis. +# +# Run this script with: +# >>> python single_astrophot_model.py target_image.fits [OPTIONS] +# ============================================================================= + +import astrophot as ap +import numpy as np +from astropy.io import fits +from astropy.wcs import WCS +import argparse +import ast + +try: + import yaml +except ImportError: + from astropy.io.misc import yaml + + +def parse_arbitrary_args(unknown_args): + """ + Parses unrecognized command-line arguments into a dictionary. + Attempts to cast values to native Python types (int, float, list, etc.). + """ + params = {} + i = 0 + while i < len(unknown_args): + arg = unknown_args[i] + if arg.startswith("--"): + key = arg[2:] + # Handle --key=value format + if "=" in key: + key, val_str = key.split("=", 1) + # Handle --key value format + elif i + 1 < len(unknown_args) and not unknown_args[i + 1].startswith("--"): + val_str = unknown_args[i + 1] + i += 1 + else: + val_str = "True" # Assume a boolean flag if no value is provided + + # Safely evaluate strings into numbers, lists, or booleans + try: + val = ast.literal_eval(val_str) + except (ValueError, SyntaxError): + val = val_str # Fallback to keeping it as a string + + params[key] = val + i += 1 + return params + + +def to_serializable(value): + if value is None: + return None + if hasattr(value, "detach"): + value = value.detach().cpu().numpy() + value = np.array(value) + return value.item() if value.shape == () else value.tolist() + + +def collect_model_parameters(model): + params = model.dynamic_params + names = tuple(param.name for param in params) + values = tuple(to_serializable(param.npvalue) for param in params) + uncertainties = tuple(to_serializable(param.uncertainty) for param in params) + return { + name: {"value": val, "uncertainty": unc} + for name, val, unc in zip(names, values, uncertainties) + } + + +def main(): + parser = argparse.ArgumentParser(description="Fit a model to a target image using AstroPhot.") + + # Core Parameters + parser.add_argument("target_file", type=str, help="Path to the target FITS file") + parser.add_argument( + "--window", + type=int, + default=None, + nargs=4, + metavar=("IMIN", "IMAX", "JMIN", "JMAX"), + help="Window border for the fit in pixel coordinates, given as: --window imin imax jmin jmax", + ) + parser.add_argument( + "--name", type=str, default="astrophot_model", help="Name used for saving files" + ) + parser.add_argument("--psf", type=str, default=None, help="Path to the PSF FITS file") + parser.add_argument( + "--psf_upsample", type=int, default=1, help="PSF upsampling factor for convolution (int)" + ) + parser.add_argument("--zeropoint", type=float, default=None, help="Magnitude zeropoint") + parser.add_argument( + "--initial_sky", + type=float, + default=None, + help="Initial sky value for the I0 parameter", + ) + parser.add_argument( + "--sky_locked", action="store_true", help="Lock the sky model during fitting" + ) + parser.add_argument( + "--model_type", + type=str, + default="sersic_galaxy_model", + help="Type of AstroPhot model to fit. Replace spaces with underscores, e.g. 'sersic_galaxy_model' or 'exponential_disk_model'", + ) + parser.add_argument("--verbose", type=int, default=1, help="Verbosity level for fitting output") + parser.add_argument( + "--dump", + action="store_true", + help="Dump this file to the current directory as 'single_astrophot_model.py' for editing and running", + ) + + # Variance and Mask Parameters + parser.add_argument("--variance", type=str, default=None, help="Path to the variance FITS file") + parser.add_argument( + "--variance_hdu", type=int, default=0, help="FITS file index for variance data" + ) + parser.add_argument("--mask", type=str, default=None, help="Path to the mask FITS file") + parser.add_argument("--mask_hdu", type=int, default=0, help="FITS file index for mask data") + + # Extra Parameters + parser.add_argument( + "--no_save_images", + action="store_false", + dest="save_images", + help="Disable saving the model and residual images", + ) + parser.add_argument( + "--no_save_cov", + action="store_false", + dest="save_covariance_matrix", + help="Disable saving the covariance matrix", + ) + parser.add_argument( + "--target_hdu", type=int, default=0, help="FITS file index for target image data" + ) + parser.add_argument("--psf_hdu", type=int, default=0, help="FITS file index for PSF data") + parser.add_argument( + "--sky_model_type", + type=str, + default="flat", + help="Type of sky model to fit, options include 'flat' or 'plane'", + ) + + # Parse known arguments, leave the rest for the dynamic dictionary + args, unknown = parser.parse_known_args() + + # Parse the arbitrary parameters for the model + initial_params = parse_arbitrary_args(unknown) + + if args.verbose == 0: + ap.config.set_logging_output(stdout=False, filename=None) + + if args.dump: + with open("single_astrophot_model.py", "w", encoding="utf-8") as f: + f.write( + "# This file was generated by dumping the single_model_cli.py script. You can edit and run this file directly.\n\n" + ) + with open(__file__, "r", encoding="utf-8") as original: + f.write(original.read()) + print("Dumped to single_astrophot_model.py") + + # Load Target Data + # --------------------------------------------------------------------- + if args.verbose > 0: + print("Loading target image...") + with fits.open(args.target_file) as hdu: + target_data = np.array(hdu[args.target_hdu].data, dtype=np.float64) + + # Load Variance Data + # --------------------------------------------------------------------- + variance_data = None + if args.variance is not None: + if args.verbose > 0: + print("Loading variance image...") + with fits.open(args.variance) as hdu: + variance_data = np.array(hdu[args.variance_hdu].data, dtype=np.float64) + + # Load Mask Data + # --------------------------------------------------------------------- + mask_data = None + if args.mask is not None: + if args.verbose > 0: + print("Loading mask image...") + with fits.open(args.mask) as hdu: + mask_data = np.array(hdu[args.mask_hdu].data) + + # Load PSF + # --------------------------------------------------------------------- + psf_data = None + if args.psf is not None: + if args.verbose > 0: + print("Loading PSF...") + with fits.open(args.psf) as hdu: + psf_data = np.array(hdu[args.psf_hdu].data, dtype=np.float64) + psf_data = ap.PSFImage(data=psf_data, upsample=args.psf_upsample) + + # Make Target + # --------------------------------------------------------------------- + target_wcs = WCS(fits.getheader(args.target_file, args.target_hdu)) + target = ap.TargetImage( + data=target_data, + wcs=target_wcs, + zeropoint=args.zeropoint, + variance=variance_data, + mask=mask_data, + psf=psf_data, + ) + + # Create Model + # --------------------------------------------------------------------- + model_object = ap.Model( + name=args.name, + model_type=args.model_type.replace("_", " "), + target=target, + **initial_params, + window=args.window, + ) + + model_sky = ap.Model( + name="sky", + model_type=args.sky_model_type + " sky model", + target=target, + I0=args.initial_sky, + window=args.window, + ) + + if args.sky_locked: + model_sky.to_static() + + model = ap.Model( + name="astrophot_model", + model_type="group model", + target=target, + models=[model_sky, model_object], + ) + + # Fit the model + # --------------------------------------------------------------------- + if args.verbose > 0: + print("Initializing model...") + model.initialize() + if args.verbose > 0: + print("Fitting model...") + result = ap.fit.LM(model, verbose=args.verbose).fit() + + # Report Total Magnitude or Flux + # ---------------------------------------------------------------------- + total_flux = to_serializable(model_object.total_flux()) + total_flux_uncertainty = to_serializable(model_object.total_flux_uncertainty()) + total_magnitude = None + total_magnitude_uncertainty = None + if args.zeropoint is not None: + total_magnitude = to_serializable(model_object.total_magnitude()) + total_magnitude_uncertainty = to_serializable(model_object.total_magnitude_uncertainty()) + if args.verbose > 0: + print(model) + if args.zeropoint is not None: + print(f"Total Magnitude: {total_magnitude} +- {total_magnitude_uncertainty}") + else: + print(f"Total Flux: {total_flux} +- {total_flux_uncertainty}") + + # Save Results + # ---------------------------------------------------------------------- + output_summary = { + model_object.name: { + "model_type": model_object.model_type, + "parameters": collect_model_parameters(model_object), + "total_flux": total_flux, + "total_flux_uncertainty": total_flux_uncertainty, + "total_magnitude": total_magnitude, + "total_magnitude_uncertainty": total_magnitude_uncertainty, + "note": "Total flux/mag is within the fitting window, not extended to infinity", + }, + "sky_model": { + "model_type": model_sky.model_type, + "parameters": collect_model_parameters(model_sky), + }, + } + if len(model_object.static_params) > 0: + output_summary[model_object.name]["static_parameters"] = [ + param.name for param in model_object.static_params + ] + + with open(f"{args.name}_parameters.yaml", "w", encoding="utf-8") as output_file: + yaml.dump(output_summary, output_file, default_flow_style=False) + + if args.save_images: + model().save(f"{args.name}_model_image.fits") + (target[model.window] - model()).save(f"{args.name}_residual_image.fits") + + if args.save_covariance_matrix: + cov = result.covariance_matrix + if hasattr(cov, "detach"): + cov = cov.detach().cpu().numpy() + np.save(f"{args.name}_covariance_matrix.npy", np.array(cov)) + + +if __name__ == "__main__": + main() diff --git a/docs/source/fastfit.rst b/docs/source/fastfit.rst index 08822f09..06329058 100644 --- a/docs/source/fastfit.rst +++ b/docs/source/fastfit.rst @@ -1,86 +1,91 @@ -=================== -Fit Something Fast! -=================== +============================ +Fit Fast! with AstroPhot CLI +============================ -Here are some scripts which can be used like configuration files. Simply fill in -the blanks at the top of the file and a basic fitting routine will run on your -data. This is useful to get results fast on a fairly standard dataset, or to use -as a starting point when building a more sophisticated analysis routine. +Sometimes you just have a FITS file and want to quickly fit a simple model to +it. The AstroPhot Command-Line-Interface (CLI) is designed for this purpose. It +provides a simple interface to fit an isolated model, or a catalogue of objects +all in a single line command. Fit a single isolated object ---------------------------- -Get the script here: :download:`single_model_fit.py ` +This script will fit a single model (plus a sky level) to an image. The simplest +usage is just: -This script will fit a single object with a single model (plus a sky model if -requested). +.. code:: bash -basic usage is to fill in these blanks at the top of the file. Even just filling -the ``target_file`` is enough to get started: + ~$ astrophot_single target_file.fits -.. code:: python +Which will fit a sersic model to the whole image, assuming a variance of 1 for +every pixel. This is almost certainly not what you want, but as a starting point +it is a good way check that everything is working. Here is a more realistic example: - name = "object_name" # used for saving files - target_file = ".fits" # can be a numpy array instead - mask_file = None # ".fits" # can be a numpy array instead - psf_file = None # ".fits" # can be a numpy array instead - variance_file = None # ".fits" # or numpy array or "auto" - pixelscale = 0.1 # arcsec/pixel - zeropoint = 22.5 # mag - initial_params = None # e.g. {"center": [3, 3], "q": {"value": 0.8, "locked": True}} - window = None # None to fit whole image, otherwise ((xmin,xmax),(ymin,ymax)) pixels - initial_sky = None # If None, sky will be estimated - sky_locked = False - model_type = "sersic galaxy model" +.. code:: bash -then run the script from the command line as a python file: + ~$ astrophot_single target_file.fits --psf psf.fits --mask mask.fits --variance variance.fits --window 10 100 10 100 --zeropoint 22.5 -.. code:: bash +In this version we pass FITS files for the PSF, mask, and variance, we only fit +a 90x90 pixel window of the image, and we set the zeropoint for magnitude +calculations. See the :doc:`tutorials/GettingStarted` tutorial or use ``--help`` +for more information on these inputs. - >>> python single_model_fit.py +The final output will be a yaml file with the fitted parameters. Depending on +the options you choose, it may also output fits files for the model and residual +images as well as a numpy file with a covariance matrix of uncertainties for the +parameters. -This will output the fitted parameters and save the model and residual images as -fits files. See the :doc:`tutorials/GettingStarted` tutorial for more -information. +**Note**: Unfortunately, this script is really slowed down by the fact that +PyTorch needs to be loaded for each run, which can take several seconds. For +faster bulk processing see the section below. Fit all objects in an image --------------------------- -Get the script here: :download:`segmap_models_fit.py ` - -This script will fit all objects in an image with a single model type. It will -also fit a sky model (if requested) and a single special model as the "primary -obejct" (if requested). - -basic usage is to fill in these blanks at the top of the file. Even just filling -the ``target_file`` and ``segmap_file`` is enough to get started: - -.. code:: python - - name = "field_name" # used for saving files - target_file = ".fits" # can be a numpy array instead - segmap_file = ".fits" # can be a numpy array instead - mask_file = None # ".fits" # can be a numpy array instead - psf_file = None # ".fits" # can be a numpy array instead - variance_file = None # ".fits" # or numpy array or "auto" - pixelscale = 0.1 # arcsec/pixel - zeropoint = 22.5 # mag - initial_sky = None # If None, sky will be estimated. Recommended to set manually - sky_locked = False - model_type = "sersic galaxy model" # model type for segmap entries - segmap_filter = {} # in pixels or ADU: min_size, min_area, min_flux - segmap_filter_ids = [] # list of segmap ids to remove from fit - segmap_override_init_params = {} # Override some initial parameters for segmap models - primary_key = None # segmentation map id, use None to have no primary object - primary_name = "primary object" # name for primary object - primary_model_type = "sersic galaxy model" - primary_initial_params = None # {"center": [3, 3], "q": {"value": 0.8, "locked": True}} - -then run the script from the command line as a python file: +This script will fit all objects in an image with a single model type. It +functions similarly to the single object fitting script, but instead of fitting +a single model it will loop through and fit everything identified in a +segmentation map. Note that this is not a group model fitting script. So each +segmentation will be fit separately rather than as a joint likelihood. If you +wish to account for blending and fit a group model, see the +:doc:`tutorials/GroupModels` tutorial. + +The simplest use case is just: .. code:: bash - >>> python segmap_models_fit.py + ~$ astrophot_segmap target_file.fits segmap_file.fits + +Just like in the single object case, this is not likely to be what you want, but +it is a good starting point to check that everything is working. A more +realistic example is: + +.. code:: bash -This will output the fitted parameters and save the model and residual images as -fits files. See the :doc:`tutorials/GroupModels` tutorial for more information. + ~$ astrophot_segmap target_file.fits segmap_file.fits --psf psf.fits --mask mask.fits --variance variance.fits --window_expand_border 5 --zeropoint 22.5 + +In this version we pass FITS files for the PSF, mask, and variance, we expand +the fitting windows by 5 pixels in each direction, and we set the zeropoint for +magnitude calculations. See the :doc:`tutorials/GettingStarted` tutorial or use +``--help`` for more information on these inputs. + +If you are looking for a bit more control, you can also pass a catalogue of +initial parameters for each object using ``--cat catalogue.yaml``. This is a +yaml file with keys corresponding to the segmentation ids and values that are +dictionaries of parameters to override the default initial parameters. For +example: + +.. code:: yaml + + 1: + center: [10, 10] + q: + value: 0.8 + locked: true + 2: + center: [20, 20] + q: 0.6 + +The final output will be a yaml file with the fitted parameters for each object. +Depending on your configuration, the covariance matrix of uncertainties and fits +files for the model and residuals may also be included. diff --git a/docs/source/prebuilt/segmap_models_fit.py b/docs/source/prebuilt/segmap_models_fit.py deleted file mode 100644 index 389fe8ae..00000000 --- a/docs/source/prebuilt/segmap_models_fit.py +++ /dev/null @@ -1,228 +0,0 @@ -# ============================================================================= -# Fit all objects identified in a segmentation map -# -# This is a quick script to fit all the objects identified in a segmentation map -# using a single model type. You should set the parameters under PARAMETERS to -# be appropriate for your data. The script will load the target image, mask, -# psf, and variance image (if available) and fit the models to the target image. -# -# First a fit will be run on tight windows exactly enclosing the segmentations -# for each object. Then the windows will be expanded by the set factors and the -# fit will be run again. This is more stable than fitting the expanded windows -# from the start since it reduces the effects of overlap -# -# Run this script with: -# >>> python segmap_models_fit.py -# ============================================================================= - -import astrophot as ap -import numpy as np -from astropy.io import fits -import matplotlib.pyplot as plt - -# PARAMETERS -###################################################################### -name = "field_name" # used for saving files -target_file = ".fits" # can be a numpy array instead -segmap_file = ".fits" # can be a numpy array instead -psf_file = None # ".fits" # can be a numpy array instead -zeropoint = 22.5 # mag -initial_sky = None # If None, sky will be estimated. Recommended to set manually -sky_locked = False -model_type = "sersic galaxy model" # model type for segmap entries -segmap_filter = {} # in pixels or ADU: min_size, min_area, min_flux -segmap_filter_ids = [] # list of segmap ids to remove from fit -segmap_override_init_params = {} # Override some initial parameters for segmap models -primary_key = None # segmentation map id, use None to have no primary object -primary_name = "primary object" # name for primary object -primary_model_type = "spline galaxy model" -primary_initial_params = {} # {"center": [3, 3], "q": 0.8} -# Extra parameters -###################################################################### -save_model_image = True -save_residual_image = True -target_hdu = 0 # FITS file index for image data -segmap_hdu = 0 -psf_hdu = 0 -window_expand_scale = 2 # Windows from segmap will be expanded by this factor -window_expand_border = 10 # Windows from segmap will be expanded by this number of pixels -sky_model_type = "flat sky model" -print_all_models = True -###################################################################### - -# load target and segmentation map -# --------------------------------------------------------------------- -print("loading target and segmentation map") -target = ap.TargetImage( - filename=target_file, - hduext=target_hdu, - zeropoint=zeropoint, -) - -if isinstance(segmap_file, str): - hdu = fits.open(segmap_file) - segmap_data = np.array(hdu[segmap_hdu].data, dtype=np.int32) -else: - segmap_data = segmap_file - -# load psf -# --------------------------------------------------------------------- -# PSF -if isinstance(psf_file, str): - print("loading psf") - hdu = fits.open(psf_file) - psf_data = np.array(hdu[psf_hdu].data, dtype=np.float64) - target.psf = target.psf_image(data=psf_data) -elif psf_file is None: - psf = None -else: - target.psf = target.psf_image(data=psf_file) - -# Initialization from segmap -# --------------------------------------------------------------------- -print("Parsing segmentaiton map") -windows = ap.utils.initialize.windows_from_segmentation_map(segmap_data) -if len(segmap_filter) > 0: - windows = ap.utils.initialize.filter_windows( - windows, - **segmap_filter, - image=target, - ) - -for ids in segmap_filter_ids: - del windows[ids] -centers = ap.utils.initialize.centroids_from_segmentation_map(segmap_data, target) -if "galaxy" in model_type: - PAs = ap.utils.initialize.PA_from_segmentation_map(segmap_data, target, centers) - qs = ap.utils.initialize.q_from_segmentation_map(segmap_data, target, centers) -else: - PAs = None - qs = None -init_params = {} -for window in windows: - init_params[window] = {"center": centers[window]} - if "galaxy" in model_type: - init_params[window]["PA"] = PAs[window] - init_params[window]["q"] = qs[window] - init_params[window].update(segmap_override_init_params) - -# Create Models -# --------------------------------------------------------------------- -print("Creating models") -models = [] -models.append( - ap.Model( - name="sky", - model_type=sky_model_type, - target=target, - I0=initial_sky if initial_sky is not None else {}, - ) -) -if sky_locked: - models[0].to_static() -primary_model = None -for window in windows: - if primary_key is not None and window == primary_key: - print(primary_name, window) - if "center" not in primary_initial_params: - primary_initial_params["center"] = init_params[window]["center"] - if ( - "PA" not in primary_initial_params - and PAs is not None - and "galaxy" in primary_model_type - ): - primary_initial_params["PA"] = PAs[window] - if "q" not in primary_initial_params and qs is not None and "galaxy" in primary_model_type: - primary_initial_params["q"] = qs[window] - model = ap.Model( - name=primary_name, - model_type=primary_model_type, - target=target, - **primary_initial_params, - window=windows[window], - ) - primary_model = model - else: - print(window) - model = ap.Model( - name=f"{model_type}_{window}", - model_type=model_type, - target=target, - window=windows[window], - **init_params[window], - ) - models.append(model) -model = ap.Model( - name=f"{name}_model", - model_type="group model", - target=target, - models=models, -) - -# Fit the model -# --------------------------------------------------------------------- -print("Initializing model") -model.initialize() -print("Fitting model round 1") -result = ap.fit.Iter(model, verbose=1).fit() -print("expanding windows") -windows = ap.utils.initialize.scale_windows( - windows, - image=target, - expand_scale=window_expand_scale, - expand_border=window_expand_border, -) -for i, window in enumerate(windows): - models[i + 1].window = windows[window] -print("Fitting round 2") -result = ap.fit.Iter(model, verbose=1).fit() - -# Report Results -# ---------------------------------------------------------------------- -if not sky_locked: - print(models[0]) - -if not primary_model is None: - print(primary_model) - totmag = primary_model.total_magnitude().detach().cpu().numpy() - print(f"Total Magnitude: {totmag}") - if hasattr(primary_model, "radial_model"): - fig, ax = plt.subplots(figsize=(8, 8)) - ap.plots.radial_light_profile(fig, ax, primary_model) - plt.savefig(f"{name}_radial_light_profile.jpg") - plt.close() - with open(f"{name}_primary_params.csv", "w") as f: - f.write("Name,Total Magnitude," + ",".join(primary_model.build_params_array_names()) + "\n") - f.write("string,mag," + ",".join(primary_model.build_params_array_units()) + "\n") - params = primary_model.get_values().detach().cpu().numpy() - f.write(",".join([str(x) for x in params]) + "\n") - -if print_all_models: - print(model) - segmap_params = [] - for segmodel in models[1:]: - if segmodel.name == primary_name: - continue - totmag = segmodel.total_magnitude().detach().cpu().numpy() - segmap_params.append( - [segmodel.name, totmag] + list(segmodel.get_values().detach().cpu().numpy()) - ) - with open(f"{name}_segmap_params.csv", "w") as f: - f.write("Name,Total Magnitude," + ",".join(segmodel.build_params_array_names()) + "\n") - f.write("string,mag," + ",".join(segmodel.build_params_array_units()) + "\n") - for row in segmap_params: - f.write(",".join([str(x) for x in row]) + "\n") - -model.save_state(f"{name}_parameters.hdf5") -if save_model_image: - model().save(f"{name}_model_image.fits") - fig, ax = plt.subplots() - ap.plots.model_image(fig, ax, model) - plt.savefig(f"{name}_model_image.jpg") - plt.close() -if save_residual_image: - (target - model()).save(f"{name}_residual_image.fits") - fig, ax = plt.subplots() - ap.plots.residual_image(fig, ax, model, normalize_residuals=True) - plt.savefig(f"{name}_residual_image.jpg") - plt.close() diff --git a/docs/source/prebuilt/single_model_fit.py b/docs/source/prebuilt/single_model_fit.py deleted file mode 100644 index ffeea84d..00000000 --- a/docs/source/prebuilt/single_model_fit.py +++ /dev/null @@ -1,130 +0,0 @@ -# ============================================================================= -# Fit a single model to a target image -# -# This is a quick script to fit a single model to a target image. You should set -# the parameters under PARAMETERS to be appropriate for your data. The script -# will load the target image, mask, psf, and variance image (if available) and -# fit the model to the target image. The script will save the model image, -# residual image, and covariance matrix to the current directory. This script is -# intended for quick easy fits, users more comfortable with configuration file -# style behavior, and as a starting point to build a more complex analysis. -# -# Run this script with: -# >>> python single_model_fit.py -# ============================================================================= - -import astrophot as ap -import numpy as np -from astropy.io import fits -import matplotlib.pyplot as plt - -# PARAMETERS -###################################################################### -name = "object_name" # used for saving files -target_file = ".fits" # can be a numpy array instead -psf_file = None # ".fits" # can be a numpy array instead -zeropoint = 22.5 # mag -initial_params = {} # e.g. {"center": [3, 3], "q": 0.8} -window = None # None to fit whole image, otherwise (xmin,xmax,ymin,ymax) pixels -initial_sky = None # If None, sky will be estimated -sky_locked = False -model_type = "sersic galaxy model" -# Extra parameters -###################################################################### -save_model_image = True -save_residual_image = True -save_covariance_matrix = True -target_hdu = 0 # FITS file index for image data -psf_hdu = 0 -sky_model_type = "flat sky model" -###################################################################### - -# load target -# --------------------------------------------------------------------- -print("loading target") -target = ap.TargetImage( - filename=target_file, - hduext=target_hdu, - zeropoint=zeropoint, -) - -# PSF -if isinstance(psf_file, str): - print("loading psf") - hdu = fits.open(psf_file) - psf_data = np.array(hdu[psf_hdu].data, dtype=np.float64) - target.psf = target.psf_image(data=psf_data) -elif psf_file is None: - psf = None -else: - target.psf = target.psf_image(data=psf_file) - -# Create Model -# --------------------------------------------------------------------- -model_object = ap.Model( - name=name, - model_type=model_type, - target=target, - psf_convolve=True if psf_file is not None else False, - **initial_params, - window=window, -) -model_sky = ap.Model( - name="sky", - model_type=sky_model_type, - target=target, - I0=initial_sky if initial_sky is not None else {}, - window=window, -) -if sky_locked: - model_sky.to_static() -model = ap.Model( - name="astrophot_model", - model_type="group model", - target=target, - models=[model_sky, model_object], -) - -# Fit the model -# --------------------------------------------------------------------- -print("Initializing model") -model.initialize() -print("Fitting model") -result = ap.fit.LM(model, verbose=1).fit() - -# Report Results -# ---------------------------------------------------------------------- -print(model) -totmag = model_object.total_magnitude().detach().cpu().numpy() -totmag_err = model_object.total_magnitude_uncertainty().detach().cpu().numpy() -print(f"Total Magnitude: {totmag} +- {totmag_err}") - -model.save_state(f"{name}_parameters.hdf5") -if save_model_image: - model().save(f"{name}_model_image.fits") - fig, ax = plt.subplots() - ap.plots.model_image(fig, ax, model) - plt.savefig(f"{name}_model_image.jpg") - plt.close() - if hasattr(model_object, "radial_model"): - fig, ax = plt.subplots(figsize=(8, 8)) - ap.plots.radial_light_profile(fig, ax, model_object) - plt.savefig(f"{name}_radial_light_profile.jpg") - plt.close() -if save_residual_image: - (target - model()).save(f"{name}_residual_image.fits") - fig, ax = plt.subplots() - ap.plots.residual_image(fig, ax, model, normalize_residuals=True) - plt.savefig(f"{name}_residual_image.jpg") - plt.close() - -if save_covariance_matrix: - np.save(f"{name}_covariance_matrix.npy", result.covariance_matrix.detach().cpu().numpy()) - fig, ax = ap.plots.covariance_matrix( - result.covariance_matrix.detach().cpu().numpy(), - model.get_values().detach().cpu().numpy(), - model.build_params_array_names(), - ) - fig.suptitle("Parameter Covariance") - plt.savefig(f"{name}_covariance_matrix.pdf") - plt.close() diff --git a/pyproject.toml b/pyproject.toml index 58db01a5..b5a9aa4f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -66,6 +66,10 @@ dev = [ "pyvo" ] +[project.scripts] +astrophot_single = "astrophot.cli.cli_single:main" +astrophot_segmap = "astrophot.cli.cli_segmap:main" + [tool.hatch.version] source = "vcs" diff --git a/tests/test_cli.py b/tests/test_cli.py new file mode 100644 index 00000000..562e017e --- /dev/null +++ b/tests/test_cli.py @@ -0,0 +1,131 @@ +import os +import shutil +import subprocess +import numpy as np +from astropy.io import fits +from photutils.segmentation import detect_sources, deblend_sources + + +def test_single_cli(): + + # Target image for testing + target_image_src = os.path.join( + os.path.split(os.path.dirname(__file__))[0], + "docs", + "source", + "tutorials", + "target_image.fits", + ) + + # Copy them to current working directory (which is where pytest runs, usually root or tests/) + target_image = "target_image.fits" + + shutil.copy(target_image_src, target_image) + + try: + # Run the CLI script + subprocess.run( + [ + "astrophot_single", + target_image, + "--name", + "test_cli_output", + "--verbose", + "1", + "--model_type", + "sersic_galaxy_model", + "--zeropoint", + "25.0", + "--dump", + ], + check=True, + ) + + # Verify outputs are created + assert os.path.exists("test_cli_output_parameters.yaml") + assert os.path.exists("test_cli_output_model_image.fits") + assert os.path.exists("test_cli_output_residual_image.fits") + assert os.path.exists("test_cli_output_covariance_matrix.npy") + assert os.path.exists("single_astrophot_model.py") + + finally: + # Cleanup + for file in [ + target_image, + "test_cli_output_parameters.yaml", + "test_cli_output_model_image.fits", + "test_cli_output_residual_image.fits", + "test_cli_output_covariance_matrix.npy", + "single_astrophot_model.py", + ]: + if os.path.exists(file): + os.remove(file) + + +def test_segmap_cli(): + + # Target image for testing + target_image_src = os.path.join( + os.path.split(os.path.dirname(__file__))[0], + "docs", + "source", + "tutorials", + "group_target_image.fits", + ) + + # Copy them to current working directory (which is where pytest runs, usually root or tests/) + target_image = "group_target_image.fits" + + shutil.copy(target_image_src, target_image) + + with fits.open(target_image) as hdu: + target_data = np.array(hdu[0].data, dtype=np.float64) + initsegmap = detect_sources(target_data, threshold=0.02, npixels=6) + segmap = deblend_sources(target_data, initsegmap, npixels=5).data + fits.writeto("group_segmap.fits", segmap, overwrite=True) + + try: + # Run the CLI script + subprocess.run( + [ + "astrophot_segmap", + target_image, + "group_segmap.fits", + "--name", + "test_segmap_cli", + "--verbose", + "1", + "--model_type", + "sersic_galaxy_model", + "--window_expand_border", + "5", + "--filter_ids", + "0", + "8", + "--zeropoint", + "25.0", + "--dump", + ], + check=True, + ) + + # Verify outputs are created + assert os.path.exists("test_segmap_cli_parameters.yaml") + assert os.path.exists("test_segmap_cli_model_images.fits") + assert os.path.exists("test_segmap_cli_residual_images.fits") + assert os.path.exists("test_segmap_cli_covariance_matrix.npz") + assert os.path.exists("segmap_astrophot_model.py") + + finally: + # Cleanup + for file in [ + target_image, + "test_segmap_cli_parameters.yaml", + "test_segmap_cli_model_images.fits", + "test_segmap_cli_residual_images.fits", + "test_segmap_cli_covariance_matrix.npz", + "segmap_astrophot_model.py", + "group_segmap.fits", + ]: + if os.path.exists(file): + os.remove(file) diff --git a/tests/test_segmap_cli_output_covariance_matrix.npz b/tests/test_segmap_cli_output_covariance_matrix.npz new file mode 100644 index 00000000..1895b3fb Binary files /dev/null and b/tests/test_segmap_cli_output_covariance_matrix.npz differ