level-2b_tweakreg_align-to-gaia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the level-2b pipeline (```calwebb_image2```)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Documentation: https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_image2.html#calwebb-image2-stage-2-imaging-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-04-29 19:10:26,140 - stpipe.Image2Pipeline - INFO - Image2Pipeline instance created.\n",
      "2022-04-29 19:10:26,142 - stpipe.Image2Pipeline.bkg_subtract - INFO - BackgroundStep instance created.\n",
      "2022-04-29 19:10:26,144 - stpipe.Image2Pipeline.assign_wcs - INFO - AssignWcsStep instance created.\n",
      "2022-04-29 19:10:26,145 - stpipe.Image2Pipeline.flat_field - INFO - FlatFieldStep instance created.\n",
      "2022-04-29 19:10:26,147 - stpipe.Image2Pipeline.photom - INFO - PhotomStep instance created.\n",
      "2022-04-29 19:10:26,149 - stpipe.Image2Pipeline.resample - INFO - ResampleStep instance created.\n",
      "2022-04-29 19:10:26,201 - stpipe.Image2Pipeline - INFO - Step Image2Pipeline running with args ('v1/jw01144030001_02101_00004_nrcblong_cal.fits',).\n",
      "2022-04-29 19:10:26,206 - stpipe.Image2Pipeline - INFO - Step Image2Pipeline parameters are: {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': True, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': '', 'save_bsub': False, 'steps': {'bkg_subtract': {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': '', 'save_combined_background': False, 'sigma': 3.0, 'maxiters': None}, 'assign_wcs': {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': '', 'sip_approx': True, 'sip_max_pix_error': 0.25, 'sip_degree': None, 'sip_max_inv_pix_error': 0.25, 'sip_inv_degree': None, 'sip_npoints': 32, 'slit_y_low': -0.55, 'slit_y_high': 0.55}, 'flat_field': {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': '', 'save_interpolated_flat': False, 'user_supplied_flat': None, 'inverse': False}, 'photom': {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': '', 'inverse': False, 'source_type': None}, 'resample': {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': '', 'pixfrac': 1.0, 'kernel': 'square', 'fillval': 'INDEF', 'weight_type': 'ivm', 'output_shape': None, 'crpix': None, 'crval': None, 'rotation': None, 'pixel_scale_ratio': 1.0, 'pixel_scale': None, 'single': False, 'blendheaders': True, 'allowed_memory': None}}}\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reducing file v1/jw01144030001_02101_00004_nrcblong_cal.fits with distortion file /grp/crds/cache/references/jwst/jwst_nircam_distortion_0133.asdf\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-04-29 19:10:27,494 - stpipe.Image2Pipeline - INFO - Prefetching reference files for dataset: 'jw01144030001_02101_00004_nrcblong_cal.fits' reftypes = ['area', 'camera', 'collimator', 'dflat', 'disperser', 'distortion', 'drizpars', 'fflat', 'filteroffset', 'flat', 'fore', 'fpa', 'ifufore', 'ifupost', 'ifuslicer', 'msa', 'ote', 'photom', 'regions', 'sflat', 'specwcs', 'wavelengthrange', 'wfssbkg']\n",
      "2022-04-29 19:10:31,131 - stpipe.Image2Pipeline - INFO - Prefetch for AREA reference file is '/grp/crds/cache/references/jwst/jwst_nircam_area_0026.fits'.\n",
      "2022-04-29 19:10:31,156 - stpipe.Image2Pipeline - INFO - Prefetch for CAMERA reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,157 - stpipe.Image2Pipeline - INFO - Prefetch for COLLIMATOR reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,157 - stpipe.Image2Pipeline - INFO - Prefetch for DFLAT reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,158 - stpipe.Image2Pipeline - INFO - Prefetch for DISPERSER reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,159 - stpipe.Image2Pipeline - INFO - Prefetch for DISTORTION reference file is '/grp/crds/cache/references/jwst/jwst_nircam_distortion_0137.asdf'.\n",
      "2022-04-29 19:10:31,187 - stpipe.Image2Pipeline - INFO - Prefetch for DRIZPARS reference file is '/grp/crds/cache/references/jwst/jwst_nircam_drizpars_0001.fits'.\n",
      "2022-04-29 19:10:31,216 - stpipe.Image2Pipeline - INFO - Prefetch for FFLAT reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,218 - stpipe.Image2Pipeline - INFO - Prefetch for FILTEROFFSET reference file is '/grp/crds/cache/references/jwst/jwst_nircam_filteroffset_0001.asdf'.\n",
      "2022-04-29 19:10:31,247 - stpipe.Image2Pipeline - INFO - Prefetch for FLAT reference file is '/grp/crds/cache/references/jwst/jwst_nircam_flat_0267.fits'.\n",
      "2022-04-29 19:10:31,277 - stpipe.Image2Pipeline - INFO - Prefetch for FORE reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,278 - stpipe.Image2Pipeline - INFO - Prefetch for FPA reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,279 - stpipe.Image2Pipeline - INFO - Prefetch for IFUFORE reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,279 - stpipe.Image2Pipeline - INFO - Prefetch for IFUPOST reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,280 - stpipe.Image2Pipeline - INFO - Prefetch for IFUSLICER reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,280 - stpipe.Image2Pipeline - INFO - Prefetch for MSA reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,281 - stpipe.Image2Pipeline - INFO - Prefetch for OTE reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,281 - stpipe.Image2Pipeline - INFO - Prefetch for PHOTOM reference file is '/grp/crds/cache/references/jwst/jwst_nircam_photom_0095.fits'.\n",
      "2022-04-29 19:10:31,307 - stpipe.Image2Pipeline - INFO - Prefetch for REGIONS reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,307 - stpipe.Image2Pipeline - INFO - Prefetch for SFLAT reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,308 - stpipe.Image2Pipeline - INFO - Prefetch for SPECWCS reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,309 - stpipe.Image2Pipeline - INFO - Prefetch for WAVELENGTHRANGE reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,309 - stpipe.Image2Pipeline - INFO - Prefetch for WFSSBKG reference file is 'N/A'.\n",
      "2022-04-29 19:10:31,310 - stpipe.Image2Pipeline - INFO - Starting calwebb_image2 ...\n",
      "2022-04-29 19:10:31,311 - stpipe.Image2Pipeline - INFO - Processing product v1/jw01144030001_02101_00004_nrcblong\n",
      "2022-04-29 19:10:31,312 - stpipe.Image2Pipeline - INFO - Working on input v1/jw01144030001_02101_00004_nrcblong_cal.fits ...\n",
      "2022-04-29 19:10:31,505 - stpipe.Image2Pipeline.assign_wcs - INFO - Step assign_wcs running with args (<ImageModel(2048, 2048) from jw01144030001_02101_00004_nrcblong_cal.fits>,).\n",
      "2022-04-29 19:10:31,507 - stpipe.Image2Pipeline.assign_wcs - INFO - Step assign_wcs parameters are: {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': 'v1', 'sip_approx': True, 'sip_max_pix_error': 0.25, 'sip_degree': None, 'sip_max_inv_pix_error': 0.25, 'sip_inv_degree': None, 'sip_npoints': 32, 'slit_y_low': -0.55, 'slit_y_high': 0.55}\n",
      "2022-04-29 19:10:32,500 - stpipe.Image2Pipeline.assign_wcs - INFO - Update S_REGION to POLYGON ICRS  80.481693684 -69.491095700 80.584994473 -69.491741222 80.587675519 -69.456282117 80.485369803 -69.454799297\n",
      "2022-04-29 19:10:32,501 - stpipe.Image2Pipeline.assign_wcs - INFO - assign_wcs updated S_REGION to POLYGON ICRS  80.481693684 -69.491095700 80.584994473 -69.491741222 80.587675519 -69.456282117 80.485369803 -69.454799297\n",
      "2022-04-29 19:10:32,501 - stpipe.Image2Pipeline.assign_wcs - INFO - COMPLETED assign_wcs\n",
      "2022-04-29 19:10:32,582 - stpipe.Image2Pipeline.assign_wcs - INFO - Step assign_wcs done\n",
      "2022-04-29 19:10:32,658 - stpipe.Image2Pipeline.flat_field - INFO - Step flat_field running with args (<ImageModel(2048, 2048) from jw01144030001_02101_00004_nrcblong_cal.fits>,).\n",
      "2022-04-29 19:10:32,660 - stpipe.Image2Pipeline.flat_field - INFO - Step flat_field parameters are: {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': 'v1', 'save_interpolated_flat': False, 'user_supplied_flat': None, 'inverse': False}\n",
      "2022-04-29 19:10:39,000 - stpipe.Image2Pipeline.flat_field - INFO - Step flat_field done\n",
      "2022-04-29 19:10:39,092 - stpipe.Image2Pipeline.photom - INFO - Step photom running with args (<ImageModel(2048, 2048) from jw01144030001_02101_00004_nrcblong_cal.fits>,).\n",
      "2022-04-29 19:10:39,094 - stpipe.Image2Pipeline.photom - INFO - Step photom parameters are: {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': False, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': 'v1', 'inverse': False, 'source_type': None}\n",
      "2022-04-29 19:10:39,168 - stpipe.Image2Pipeline.photom - INFO - Using photom reference file: /grp/crds/cache/references/jwst/jwst_nircam_photom_0095.fits\n",
      "2022-04-29 19:10:39,169 - stpipe.Image2Pipeline.photom - INFO - Using area reference file: /grp/crds/cache/references/jwst/jwst_nircam_area_0026.fits\n",
      "2022-04-29 19:10:39,280 - stpipe.Image2Pipeline.photom - INFO - Using instrument: NIRCAM\n",
      "2022-04-29 19:10:39,281 - stpipe.Image2Pipeline.photom - INFO -  detector: NRCBLONG\n",
      "2022-04-29 19:10:39,281 - stpipe.Image2Pipeline.photom - INFO -  exp_type: NRC_IMAGE\n",
      "2022-04-29 19:10:39,281 - stpipe.Image2Pipeline.photom - INFO -  filter: F356W\n",
      "2022-04-29 19:10:39,282 - stpipe.Image2Pipeline.photom - INFO -  pupil: CLEAR\n",
      "2022-04-29 19:10:41,178 - stpipe.Image2Pipeline.photom - INFO - Pixel area map copied to output.\n",
      "2022-04-29 19:10:41,181 - stpipe.Image2Pipeline.photom - INFO - PHOTMJSR value: 0.434522\n",
      "2022-04-29 19:10:41,207 - stpipe.Image2Pipeline.photom - INFO - Step photom done\n",
      "2022-04-29 19:10:41,303 - stpipe.Image2Pipeline.resample - INFO - Step resample running with args (<ImageModel(2048, 2048) from jw01144030001_02101_00004_nrcblong_cal.fits>,).\n",
      "2022-04-29 19:10:41,305 - stpipe.Image2Pipeline.resample - INFO - Step resample parameters are: {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': False, 'output_use_index': True, 'save_results': True, 'skip': False, 'suffix': 'i2d', 'search_output_file': True, 'input_dir': 'v1', 'pixfrac': 1.0, 'kernel': 'square', 'fillval': 'INDEF', 'weight_type': 'ivm', 'output_shape': None, 'crpix': None, 'crval': None, 'rotation': None, 'pixel_scale_ratio': 1.0, 'pixel_scale': None, 'single': False, 'blendheaders': True, 'allowed_memory': None}\n",
      "2022-04-29 19:10:41,361 - stpipe.Image2Pipeline.resample - INFO - Drizpars reference file: /grp/crds/cache/references/jwst/jwst_nircam_drizpars_0001.fits\n",
      "2022-04-29 19:10:41,597 - stpipe.Image2Pipeline.resample - INFO - Resampling science data\n",
      "2022-04-29 19:10:43,816 - stpipe.Image2Pipeline.resample - INFO - Drizzling (2048, 2048) --> (2067, 2074)\n",
      "2022-04-29 19:10:45,199 - stpipe.Image2Pipeline.resample - INFO - Resampling var_rnoise\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-04-29 19:10:47,195 - stpipe.Image2Pipeline.resample - INFO - Drizzling (2048, 2048) --> (2067, 2074)\n",
      "2022-04-29 19:10:48,638 - stpipe.Image2Pipeline.resample - INFO - Resampling var_poisson\n",
      "2022-04-29 19:10:50,617 - stpipe.Image2Pipeline.resample - INFO - Drizzling (2048, 2048) --> (2067, 2074)\n",
      "2022-04-29 19:10:52,020 - stpipe.Image2Pipeline.resample - INFO - Resampling var_flat\n",
      "2022-04-29 19:10:53,959 - stpipe.Image2Pipeline.resample - INFO - Drizzling (2048, 2048) --> (2067, 2074)\n",
      "2022-04-29 19:10:55,369 - stpipe.Image2Pipeline.resample - INFO - Update S_REGION to POLYGON ICRS  80.481741540 -69.491103721 80.584990438 -69.492175465 80.587967530 -69.455879412 80.484893019 -69.454809479\n",
      "2022-04-29 19:10:55,965 - stpipe.Image2Pipeline.resample - INFO - Saved model in jw01144030001_02101_00004_nrcblong_i2d.fits\n",
      "2022-04-29 19:10:55,965 - stpipe.Image2Pipeline.resample - INFO - Step resample done\n",
      "2022-04-29 19:10:55,966 - stpipe.Image2Pipeline - INFO - Finished processing product v1/jw01144030001_02101_00004_nrcblong\n",
      "2022-04-29 19:10:55,968 - stpipe.Image2Pipeline - INFO - ... ending calwebb_image2\n",
      "2022-04-29 19:10:56,554 - stpipe.Image2Pipeline - INFO - Saved model in jw01144030001_02101_00004_nrcblong_cal.fits\n",
      "2022-04-29 19:10:56,555 - stpipe.Image2Pipeline - INFO - Step Image2Pipeline done\n"
     ]
    }
   ],
   "source": [
    "from jwst.pipeline import Image2Pipeline\n",
    "\n",
    "# filename = args.rate_image\n",
    "filename = 'v1/jw01144030001_02101_00004_nrcblong_cal.fits'\n",
    "\n",
    "\n",
    "# distortionfile = args.distortion_asdf\n",
    "distortionfile = '/grp/crds/cache/references/jwst/jwst_nircam_distortion_0133.asdf'\n",
    "\n",
    "print(f'reducing file {filename} with distortion file {distortionfile}')\n",
    "\n",
    "stage2_img = Image2Pipeline()\n",
    "stage2_img.resample.save_results = False\n",
    "stage2_img.save_results = True\n",
    "stage2_img.override_distortion = distortionfile\n",
    "\n",
    "level2_result = stage2_img.run(filename) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run only the alignment step of level-3 pipeline (```tweakreg```, which includes Gaia alignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Documentation: https://jwst-pipeline.readthedocs.io/en/latest/jwst/tweakreg/README.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, you need to edit the tweakreg code in the local JWST pipeline environment, because it's hardcoded to skip ```tweakreg``` if there is only one input image (it's expecting an association with multiple exposures). If you need help figuring out where the scripts are located, use (thanks, Bryan!):\n",
    "\n",
    "```python\n",
    "import jwst, tweakwcs\n",
    "print(jwst.tweakreg.__file__)\n",
    "print(tweakwcs.__file__)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Installation instructions are here: https://github.com/spacetelescope/jwst#installing-for-developers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Go to: ```/jwst/jwst/tweakreg/tweakreg_step.py``` and you'll see the code starting at **line 136**:\n",
    "\n",
    "```python\n",
    "        if len(grp_img) == 1:\n",
    "            self.log.info(\"* Images in GROUP 1:\")\n",
    "            for im in grp_img[0]:\n",
    "                self.log.info(\"     {}\".format(im.meta.filename))\n",
    "            self.log.info('')\n",
    "\n",
    "            # we need at least two exposures to perform image alignment\n",
    "            self.log.warning(\"At least two exposures are required for image \"\n",
    "                             \"alignment.\")\n",
    "            self.log.warning(\"Nothing to do. Skipping 'TweakRegStep'...\")\n",
    "            self.skip = True\n",
    "            for model in images:\n",
    "                model.meta.cal_step.tweakreg = \"SKIPPED\"\n",
    "                # Remove the attached catalogs\n",
    "                del model.catalog\n",
    "            return input\n",
    "```\n",
    "**Comment this whole section out.**\n",
    "\n",
    "\n",
    "\n",
    "Next, at **line 196**:\n",
    "\n",
    "```python\n",
    "except ValueError as e:\n",
    "    msg = e.args[0]\n",
    "    if (msg == \"Too few input images (or groups of images) with \"\n",
    "            \"non-empty catalogs.\"):\n",
    "        # we need at least two exposures to perform image alignment\n",
    "        self.log.warning(msg)\n",
    "        self.log.warning(\"At least two exposures are required for \"\n",
    "                         \"image alignment.\")\n",
    "        self.log.warning(\"Nothing to do. Skipping 'TweakRegStep'...\")\n",
    "        self.skip = True\n",
    "        for model in images:\n",
    "            model.meta.cal_step.tweakreg = \"SKIPPED\"\n",
    "        return images\n",
    "```    \n",
    "\n",
    "This section should be changed to: \n",
    "\n",
    "```python\n",
    "except ValueError as e:\n",
    "    raise e\n",
    "\n",
    "    # msg = e.args[0]\n",
    "    # if (msg == \"Too few input images (or groups of images) with \"\n",
    "    #         \"non-empty catalogs.\"):\n",
    "    #     # we need at least two exposures to perform image alignment\n",
    "    #     self.log.warning(msg)\n",
    "    #     self.log.warning(\"At least two exposures are required for \"\n",
    "    #                      \"image alignment.\")\n",
    "    #     self.log.warning(\"Nothing to do. Skipping 'TweakRegStep'...\")\n",
    "    #     self.skip = True\n",
    "    #     for model in images:\n",
    "    #         model.meta.cal_step.tweakreg = \"SKIPPED\"\n",
    "    #     return images\n",
    "    #\n",
    "    # else:\n",
    "    #     raise e   \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, go to ```tweakwcs/imalign.py``` and **line 601**:\n",
    "\n",
    "```python\n",
    "# check that we have enough input images:\n",
    "if (refcat is None and len(wcs_gcat) < 2) or len(wcs_gcat) == 0:\n",
    "    raise ValueError(\"Too few input images (or groups of images) with \"\n",
    "                     \"non-empty catalogs.\")\n",
    "```\n",
    "\n",
    "**Comment this whole section out.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from jwst.tweakreg import TweakRegStep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get all the available step parameters by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "        save_catalogs = boolean(default=False) # Write out catalogs?\n",
      "        catalog_format = string(default='ecsv') # Catalog output file format\n",
      "        kernel_fwhm = float(default=2.5) # Gaussian kernel FWHM in pixels\n",
      "        snr_threshold = float(default=10.0) # SNR threshold above the bkg\n",
      "        brightest = integer(default=1000) # Keep top ``brightest`` objects\n",
      "        peakmax = float(default=None) # Filter out objects with pixel values >= ``peakmax``\n",
      "        enforce_user_order = boolean(default=False) # Align images in user specified order?\n",
      "        expand_refcat = boolean(default=False) # Expand reference catalog with new sources?\n",
      "        minobj = integer(default=15) # Minimum number of objects acceptable for matching\n",
      "        searchrad = float(default=1.0) # The search radius in arcsec for a match\n",
      "        use2dhist = boolean(default=True) # Use 2d histogram to find initial offset?\n",
      "        separation = float(default=0.5) # Minimum object separation in arcsec\n",
      "        tolerance = float(default=1.0) # Matching tolerance for xyxymatch in arcsec\n",
      "        xoffset = float(default=0.0), # Initial guess for X offset in arcsec\n",
      "        yoffset = float(default=0.0) # Initial guess for Y offset in arcsec\n",
      "        fitgeometry = option('shift', 'rshift', 'rscale', 'general', default='general') # Fitting geometry\n",
      "        nclip = integer(min=0, default=3) # Number of clipping iterations in fit\n",
      "        sigma = float(min=0.0, default=3.0) # Clipping limit in sigma units\n",
      "        align_to_gaia = boolean(default=False)  # Align to GAIA catalog\n",
      "        gaia_catalog = option('GAIADR2', 'GAIADR1', default='GAIADR2')\n",
      "        min_gaia = integer(min=0, default=5) # Min number of GAIA sources needed\n",
      "        save_gaia_catalog = boolean(default=False)  # Write out GAIA catalog as a separate product\n",
      "        output_use_model = boolean(default=True)  # When saving use `DataModel.meta.filename`\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "from jwst.tweakreg import TweakRegStep \n",
    "print(TweakRegStep.spec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note that you need to use a search radius large enough that when it's divided by 10, sources are still found.**\n",
    "\n",
    "See https://jira.stsci.edu/browse/JP-2585 for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-04-29 19:15:27,356 - stpipe.TweakRegStep - INFO - TweakRegStep instance created.\n",
      "2022-04-29 19:15:27,478 - stpipe.TweakRegStep - INFO - Step TweakRegStep running with args ([<ImageModel(2048, 2048) from jw01144030001_02101_00004_nrcblong_cal.fits>],).\n",
      "2022-04-29 19:15:27,481 - stpipe.TweakRegStep - INFO - Step TweakRegStep parameters are: {'pre_hooks': [], 'post_hooks': [], 'output_file': None, 'output_dir': None, 'output_ext': '.fits', 'output_use_model': True, 'output_use_index': True, 'save_results': True, 'skip': False, 'suffix': None, 'search_output_file': True, 'input_dir': '', 'save_catalogs': False, 'catalog_format': 'ecsv', 'kernel_fwhm': 2.302, 'snr_threshold': 20.0, 'brightest': 1000, 'peakmax': None, 'enforce_user_order': False, 'expand_refcat': False, 'minobj': 15, 'searchrad': 2.0, 'use2dhist': True, 'separation': 0.5, 'tolerance': 1.0, 'xoffset': 0.0, 'yoffset': 0.0, 'fitgeometry': 'general', 'nclip': 3, 'sigma': 3.0, 'align_to_gaia': True, 'gaia_catalog': 'GAIADR2', 'min_gaia': 5, 'save_gaia_catalog': True}\n",
      "2022-04-29 19:15:32,085 - stpipe.TweakRegStep - INFO - Detected 1000 sources in jw01144030001_02101_00004_nrcblong_cal.fits.\n",
      "2022-04-29 19:15:32,103 - stpipe.TweakRegStep - INFO - \n",
      "2022-04-29 19:15:32,104 - stpipe.TweakRegStep - INFO - Number of image groups to be aligned: 1.\n",
      "2022-04-29 19:15:32,104 - stpipe.TweakRegStep - INFO - Image groups:\n",
      "2022-04-29 19:15:32,133 - stpipe.TweakRegStep - INFO - * Images in GROUP 'jw01144030001_02101_00004_nrcblong_cal':\n",
      "2022-04-29 19:15:32,134 - stpipe.TweakRegStep - INFO -      jw01144030001_02101_00004_nrcblong_cal\n",
      "2022-04-29 19:15:32,134 - stpipe.TweakRegStep - INFO - \n",
      "2022-04-29 19:15:32,135 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:32,136 - stpipe.TweakRegStep - INFO - ***** tweakwcs.imalign.align_wcs() started on 2022-04-29 19:15:32.135119\n",
      "2022-04-29 19:15:32,136 - stpipe.TweakRegStep - INFO -       Version 0.7.4\n",
      "2022-04-29 19:15:32,137 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:32,231 - stpipe.TweakRegStep - INFO - Selected image 'GROUP ID: jw01144030001_02101_00004_nrcblong_cal' as reference image\n",
      "2022-04-29 19:15:32,243 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:32,244 - stpipe.TweakRegStep - INFO - ***** tweakwcs.imalign.align_wcs() ended on 2022-04-29 19:15:32.243306\n",
      "2022-04-29 19:15:32,244 - stpipe.TweakRegStep - INFO - ***** tweakwcs.imalign.align_wcs() TOTAL RUN TIME: 0:00:00.108187\n",
      "2022-04-29 19:15:32,245 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:37,164 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:37,164 - stpipe.TweakRegStep - INFO - ***** tweakwcs.imalign.align_wcs() started on 2022-04-29 19:15:37.164349\n",
      "2022-04-29 19:15:37,165 - stpipe.TweakRegStep - INFO -       Version 0.7.4\n",
      "2022-04-29 19:15:37,165 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:37,253 - stpipe.TweakRegStep - INFO - Aligning image catalog 'GROUP ID: 987654' to the reference catalog.\n",
      "2022-04-29 19:15:37,401 - stpipe.TweakRegStep - INFO - Matching sources from 'jw01144030001_02101_00004_nrcblong_cal' catalog with sources from the reference 'Unnamed' catalog.\n",
      "2022-04-29 19:15:37,401 - stpipe.TweakRegStep - INFO - Computing initial guess for X and Y shifts...\n",
      "2022-04-29 19:15:37,422 - stpipe.TweakRegStep - INFO - Found initial X and Y shifts of 0.02073, 0.01188 with significance of 100.8 and 3081 matches.\n",
      "2022-04-29 19:15:37,425 - stpipe.TweakRegStep - INFO - Found 969 matches for 'GROUP ID: 987654'...\n",
      "2022-04-29 19:15:37,427 - stpipe.TweakRegStep - INFO - Performing 'general' fit\n",
      "2022-04-29 19:15:37,435 - stpipe.TweakRegStep - INFO - Computed 'general' fit for GROUP ID: 987654:\n",
      "2022-04-29 19:15:37,435 - stpipe.TweakRegStep - INFO - XSH: 0.0616053  YSH: -0.0636483    PROPER ROT: 0.021207    \n",
      "2022-04-29 19:15:37,436 - stpipe.TweakRegStep - INFO - <ROT>: 0.021207  SKEW: -0.00159983    ROT_X: 0.022007  ROT_Y: 0.0204071\n",
      "2022-04-29 19:15:37,436 - stpipe.TweakRegStep - INFO - <SCALE>: 0.999798  SCALE_X: 0.999797  SCALE_Y: 0.999798\n",
      "2022-04-29 19:15:37,437 - stpipe.TweakRegStep - INFO - \n",
      "2022-04-29 19:15:37,437 - stpipe.TweakRegStep - INFO - FIT RMSE: 0.0500056   FIT MAE: 0.0177509\n",
      "2022-04-29 19:15:37,438 - stpipe.TweakRegStep - INFO - Final solution based on 922 objects.\n",
      "2022-04-29 19:15:37,492 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:37,492 - stpipe.TweakRegStep - INFO - ***** tweakwcs.imalign.align_wcs() ended on 2022-04-29 19:15:37.492099\n",
      "2022-04-29 19:15:37,493 - stpipe.TweakRegStep - INFO - ***** tweakwcs.imalign.align_wcs() TOTAL RUN TIME: 0:00:00.327750\n",
      "2022-04-29 19:15:37,493 - stpipe.TweakRegStep - INFO -  \n",
      "2022-04-29 19:15:38,111 - stpipe.TweakRegStep - INFO - Saved model in jw01144030001_02101_00004_nrcblong_tweakregstep.fits\n",
      "2022-04-29 19:15:38,111 - stpipe.TweakRegStep - INFO - Step TweakRegStep done\n"
     ]
    }
   ],
   "source": [
    "tweakreg = TweakRegStep()\n",
    "tweakreg.kernel_fwhm = 2.302   \n",
    "tweakreg.snr_threshold = 20.0  \n",
    "tweakreg.brightest = 1000       \n",
    "tweakreg.save_results = True\n",
    "tweakreg.searchrad = 2.0\n",
    "tweakreg.align_to_gaia = True\n",
    "tweakreg.gaia_catalog = 'GAIADR2'\n",
    "tweakreg.save_gaia_catalog = True\n",
    "\n",
    "tweak = tweakreg.run(level2_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It runs! Check some stuff using Bryan's JWebbinar notebook snippets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal_file = 'jw01144001001_02101_00005_nrcalong_cal.fits'\n",
    "tweak_file = tweak[0].meta.filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from jwst import datamodels\n",
    "\n",
    "# Open the files using datamodels\n",
    "cal_data = datamodels.open(cal_file)\n",
    "tweak_data = datamodels.open(tweak_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal_detector_to_world = cal_data.meta.wcs.get_transform('detector', 'world')\n",
    "run_detector_to_world = tweak_data.meta.wcs.get_transform('detector', 'world')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at RA, Dec in the center of the detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = (1024, 1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(80.49804757647831, -69.52285585022852)\n",
      "(80.49805064896346, -69.52280678198792)\n"
     ]
    }
   ],
   "source": [
    "print(cal_detector_to_world(x, y))\n",
    "print(run_detector_to_world(x, y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shift in RA, Dec is (0.0111, 0.1766) arcsec\n",
      "This is (0.178, 2.849) pixels.\n"
     ]
    }
   ],
   "source": [
    "# What is the shift in the WCS before/after this step?\n",
    "delta_ra = run_detector_to_world(x, y)[0] - cal_detector_to_world(x, y)[0]\n",
    "delta_dec = run_detector_to_world(x, y)[1] - cal_detector_to_world(x, y)[1]\n",
    "print('Shift in RA, Dec is ({:.4f}, {:.4f}) arcsec'.format(delta_ra * 3600., delta_dec * 3600.))\n",
    "print('This is ({:.3f}, {:.3f}) pixels.'.format(delta_ra * 3600. / .062, delta_dec * 3600. / 0.062))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full GWCS model is contained in the ASDF extension of the file, and can be seen through the `meta` property. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tweak_data.meta.wcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several world coordinate systems available in the file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['detector', 'v2v3', 'v2v3vacorr', 'v2v3corr', 'world']"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# What coordinate frames are available?\n",
    "tweak_data.meta.wcs.available_frames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Frame2D(name=\"detector\", unit=(Unit(\"pix\"), Unit(\"pix\")), axes_names=('x', 'y'), axes_order=(0, 1))>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# What is the input frame of the WCS object?\n",
    "tweak_data.meta.wcs.input_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<CelestialFrame(name=\"world\", unit=(Unit(\"deg\"), Unit(\"deg\")), axes_names=('lon', 'lat'), axes_order=(0, 1), reference_frame=<ICRS Frame>)>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# What is the output frame of the WCS object?\n",
    "tweak_data.meta.wcs.output_frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='transforms'></a>\n",
    "#### Do the coordinate transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a transformation function to go from detector pixels to location on the sky."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(80.46244087451473, -69.54361944812793)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tweak_data.meta.wcs(50, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(80.46244087451473, -69.54361944812793)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get the transform to go from detector to world coordinates\n",
    "detector_to_world = tweak_data.meta.wcs.get_transform('detector', 'world')\n",
    "pix_ra, pix_dec = detector_to_world(50, 50)\n",
    "pix_ra, pix_dec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a function for the inverse transformation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50.01744404024873, 50.00591127827067)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "world_to_detector = tweak_data.meta.wcs.get_transform('world', 'detector')\n",
    "pix_x, pix_y = world_to_detector(pix_ra, pix_dec)\n",
    "pix_x, pix_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='find_coords'></a>\n",
    "#### Determine positions using transforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the transformation functions we defined above, we can now easily determine the RA and Dec of some randomly-chosen sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "sources_x = [925., 945., 940.]\n",
    "sources_y = [1045, 1183.2, 1120.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call the transform function\n",
    "sources_ra, sources_dec = detector_to_world(sources_x, sources_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([80.50026231, 80.50673304, 80.50374956])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sources_ra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-69.52439028, -69.52346722, -69.52381691])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sources_dec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now the opposite case: My target is at a given RA and Dec, so where is it in this image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "targ_ra = 80.51\n",
    "targ_dec = -69.53"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call the inverse transform function\n",
    "targ_x, targ_y = world_to_detector(targ_ra, targ_dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'targ_x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[0;32mIn [2]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTarget located at (x, y) = (\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m, \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m)\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;241m.\u001b[39mformat(\u001b[43mtarg_x\u001b[49m, targ_y))\n",
      "\u001b[0;31mNameError\u001b[0m: name 'targ_x' is not defined"
     ]
    }
   ],
   "source": [
    "print('Target located at (x, y) = ({}, {})'.format(targ_x, targ_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/arest/anaconda/envs/jwst/lib/python3.8/site-packages/jwst/tweakreg/__init__.py\n",
      "/Users/arest/anaconda/envs/jwst/lib/python3.8/site-packages/tweakwcs/__init__.py\n"
     ]
    }
   ],
   "source": [
    "import jwst, tweakwcs\n",
    "print(jwst.tweakreg.__file__)\n",
    "print(tweakwcs.__file__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}