MeasurePSFres_deep.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Make a catalogue with real and model PSFs + magnitude etc, for PSF testing script\n",
    "\n",
    "#! /usr/bin/env python\n",
    "\n",
    "!jupyter nbconvert --to script MeasurePSFres_deep.ipynb\n",
    "\n",
    "from __future__ import print_function\n",
    "import os\n",
    "import numpy as np\n",
    "#from read_psf_cats import read_data, band_combinations\n",
    "import fitsio\n",
    "import treecorr\n",
    "import matplotlib\n",
    "import matplotlib\n",
    "matplotlib.use('Agg') # needs to be done before import pyplot\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table,join\n",
    "import h5py as h\n",
    "from astropy.wcs import WCS\n",
    "\n",
    "from __future__ import print_function\n",
    "import os\n",
    "import sys\n",
    "import shutil\n",
    "import logging\n",
    "import datetime\n",
    "import traceback\n",
    "import numpy as np\n",
    "import copy\n",
    "import glob\n",
    "import time\n",
    "import fitsio\n",
    "#import pixmappy\n",
    "import pandas\n",
    "import galsim\n",
    "import galsim.des\n",
    "#import piff\n",
    "import ngmix\n",
    "import wget\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#haven't edited or decided which bits of this cell I need\n",
    "def parse_args():\n",
    "    import argparse\n",
    "\n",
    "    parser = argparse.ArgumentParser(description='Run PSFEx on a set of exposures')\n",
    "\n",
    "    # Directory arguments\n",
    "    parser.add_argument('--sex_dir', default='/astro/u/mjarvis/bin/',\n",
    "                        help='location of sextrator executable')\n",
    "    parser.add_argument('--piff_exe', default='/astro/u/mjarvis/.conda/envs/py2.7/bin/piffify',\n",
    "                        help='location of piffify executable')\n",
    "    parser.add_argument('--findstars_dir', default='/astro/u/mjarvis/bin',\n",
    "                        help='location wl executables')\n",
    "    parser.add_argument('--work', default='/astro/u/mjarvis/work/y3_piff',\n",
    "                        help='location of intermediate outputs')\n",
    "    parser.add_argument('--scratch', default='/data/mjarvis/y3_piff',\n",
    "                        help='location of intermediate outputs')\n",
    "    parser.add_argument('--pixmappy_dir', default='/astro/u/mjarvis/work/y3_piff/astro',\n",
    "                        help='location of pixmappy astrometric solutions')\n",
    "    parser.add_argument('--tag', default=None,\n",
    "                        help='A version tag to add to the directory name')\n",
    "\n",
    "    # Exposure inputs\n",
    "    parser.add_argument('--base_exposures',\n",
    "                        default='/astro/u/mjarvis/work/y3_piff/exposures-ccds-Y3A1_COADD.fits',\n",
    "                        help='The base file with information about the DES exposures')\n",
    "    parser.add_argument('--file', default='',\n",
    "                        help='list of exposures (in lieu of separate exps)')\n",
    "    parser.add_argument('--exps', default='', nargs='+',\n",
    "                        help='list of exposures to run')\n",
    "\n",
    "    # Configuration files\n",
    "    parser.add_argument('--sex_config',\n",
    "                        default='/astro/u/mjarvis/rmjarvis/DESWL/psf/y3.sex',\n",
    "                        help='sextractor config file')\n",
    "    parser.add_argument('--piff_config',\n",
    "                        default='/astro/u/mjarvis/rmjarvis/DESWL/psf/piff.yaml',\n",
    "                        help='piff config file')\n",
    "    parser.add_argument('--findstars_config',\n",
    "                        default='/astro/u/mjarvis/rmjarvis/DESWL/psf/y3.config',\n",
    "                        help='findstars config file')\n",
    "    parser.add_argument('--sex_params',\n",
    "                        default='/astro/u/mjarvis/rmjarvis/DESWL/psf/sex.param_piff',\n",
    "                        help='sextractor param file')\n",
    "    parser.add_argument('--sex_filter',\n",
    "                        default='/astro/u/mjarvis/rmjarvis/DESWL/psf/sex.conv',\n",
    "                        help='name of sextractor filter file')\n",
    "    parser.add_argument('--sex_nnw',\n",
    "                        default='/astro/u/mjarvis/rmjarvis/DESWL/psf/sex.nnw',\n",
    "                        help='name of sextractor star file')\n",
    "    parser.add_argument('--tapebump_file',\n",
    "                        default='/astro/u/mjarvis/rmjarvis/DESWL/psf/mask_ccdnum.txt',\n",
    "                        help='name of tape bump file')\n",
    "    parser.add_argument('--make_symlinks', default=0, type=int,\n",
    "                        help='make symlinks in output dir, rather than move files')\n",
    "    parser.add_argument('--noweight', default=False, action='store_const', const=True,\n",
    "                        help='do not try to use a weight image.')\n",
    "\n",
    "\n",
    "    # Options\n",
    "    parser.add_argument('--clear_output', default=0, type=int,\n",
    "                        help='should the output directory be cleared before writing new files?')\n",
    "    parser.add_argument('--rm_files', default=1, type=int,\n",
    "                        help='remove unpacked files after finished')\n",
    "    parser.add_argument('--use_existing', default=0, type=int,\n",
    "                        help='use previously downloaded files if they exist')\n",
    "    parser.add_argument('--blacklist', default=1, type=int,\n",
    "                        help='add failed CCDs to the blacklist')\n",
    "    parser.add_argument('--run_piff', default=1, type=int,\n",
    "                        help='run piff on files')\n",
    "    parser.add_argument('--run_sextractor', default=1, type=int,\n",
    "                        help='run sextractor to remake input catalog')\n",
    "    parser.add_argument('--run_findstars', default=1, type=int,\n",
    "                        help='force a run of findstars to get input star catalog')\n",
    "    parser.add_argument('--mag_cut', default=-1, type=float,\n",
    "                        help='remove the top mags using mag_auto')\n",
    "    parser.add_argument('--min_mag', default=-1, type=float,\n",
    "                        help='remove stars brighter than this mag')\n",
    "    parser.add_argument('--nbright_stars', default=1, type=int,\n",
    "                        help='use median of this many brightest stars for min mag')\n",
    "    parser.add_argument('--max_mag', default=0, type=float,\n",
    "                        help='only use stars brighter than this mag')\n",
    "    parser.add_argument('--use_tapebumps', default=1, type=int,\n",
    "                        help='avoid stars in or near tape bumps')\n",
    "    parser.add_argument('--tapebump_extra', default=2, type=float,\n",
    "                        help='How much extra room around tape bumps to exclude stars in units of FWHM')\n",
    "    parser.add_argument('--single_ccd', default=0, type=int,\n",
    "                        help='Only do the specified ccd (used for debugging)')\n",
    "    parser.add_argument('--reserve', default=0, type=float,\n",
    "                        help='Reserve some fraction of the good stars for testing')\n",
    "    parser.add_argument('--get_psfex', default=False, action='store_const', const=True,\n",
    "                        help='Download the PSFEx files along the way')\n",
    "    parser.add_argument('--plot_fs', default=False, action='store_const', const=True,\n",
    "                        help='Make a size-magnitude plot of the findstars output')\n",
    "    parser.add_argument('--use_ngmix', default=False, action='store_const', const=True,\n",
    "                        help='Use ngmix rather than hsm for the measurements')\n",
    "\n",
    "    args = parser.parse_args()\n",
    "    return args"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ColDefs(\n",
      "    name = 'SOURCE_NUMBER'; format = '1J'; disp = 'I10'\n",
      "    name = 'EXTENSION'; format = '1I'; disp = 'I4'\n",
      "    name = 'CATALOG_NUMBER'; format = '1J'; disp = 'I7'\n",
      "    name = 'FLAGS_PSF'; format = '1J'; disp = 'I7'\n",
      "    name = 'VECTOR_CONTEXT'; format = '2D'; disp = 'G12.6'\n",
      "    name = 'X_IMAGE'; format = '1D'; unit = 'pixel'; disp = 'F11.4'\n",
      "    name = 'Y_IMAGE'; format = '1D'; unit = 'pixel'; disp = 'F11.4'\n",
      "    name = 'DELTAX_IMAGE'; format = '1E'; unit = 'pixel'; disp = 'F11.4'\n",
      "    name = 'DELTAY_IMAGE'; format = '1E'; unit = 'pixel'; disp = 'F11.4'\n",
      "    name = 'NORM_PSF'; format = '1E'; unit = 'count'; disp = 'G12.6'\n",
      "    name = 'FWHM_PSF'; format = '1E'; unit = 'pixel'; disp = 'F8.2'\n",
      "    name = 'ELLIP_PSF'; format = '1E'; disp = 'F7.5'\n",
      "    name = 'SNR_PSF'; format = '1E'; disp = 'G10.4'\n",
      "    name = 'CHI2_PSF'; format = '1E'; disp = 'G12.6'\n",
      "    name = 'RESI_PSF'; format = '1E'; disp = 'G12.6'\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "#read in list of stars made from Sextractor and PSFEx\n",
    "star_file= \"/global/homes/a/aamon/DES/DEStests/DEEP/deeppsfs/UltraVista/UVISTA_J_21_01_16_psfex-starlist.fits\"\n",
    "\n",
    "dat = fits.open(star_file)\n",
    "cols = dat[2].columns\n",
    "print(cols)\n",
    "\n",
    "  # This has the following columns:\n",
    "    # id: The original id from the SExtractor catalog\n",
    "    # x: The x position\n",
    "    # y: The y position\n",
    "    # sky: The local sky value\n",
    "    # noise: The estimated noise.  But these are all 0, so I think this isn't being calculated.\n",
    "    # size_flags: Error flags that occurred when estimating the size\n",
    "    # mag: The magnitude from SExtractor\n",
    "    # sg: SExtractor's star/galaxy estimate.  Currently SPREAD_MODEL.  (Actually, currently none)\n",
    "    # sigma0: The shapelet sigma that results in a b_11 = 0 shapelet parameter.\n",
    "    # star_flag: 1 if findstars thought this was a star, 0 otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wget( url, file):\n",
    "    full_file = os.path.join(url,file)\n",
    "    if not os.path.isfile(full_file):\n",
    "        # Sometimes this fails with an \"http protocol error, bad status line\".\n",
    "        # Maybe from too many requests at once or something.  So we retry up to 5 times.\n",
    "        nattempts = 5\n",
    "        cmd = 'wget -q --no-check-certificate %s'%(full_file)\n",
    "        for attempt in range(1,nattempts+1):\n",
    "            if os.path.exists(full_file):\n",
    "                break\n",
    "    return full_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the files we need. These files are \n",
    "#It looks to me that the image file is the full coadd, \n",
    "image_file = wget('ftp://ftp.star.ucl.ac.uk/whartley/ultraVISTA/','UVISTA_J_21_01_16_allpaw_skysub_015_dr3_rc_v5.fits.gz')\n",
    "\n",
    "#row['root'] = root\n",
    "#row['image_file'] = image_file\n",
    "\n",
    "#usually weight is in image file but in this case, it's a separate file\n",
    "weight_file = wget('ftp://ftp.star.ucl.ac.uk/whartley/ultraVISTA/','UVISTA_J_21_01_16_allpaw_skysub_015_dr3_rc_v5.weight.fits.gz')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Not sure this is necessary, but having this information might be useful for further tests\n",
    "\n",
    "def read_image_header(row, img_file):\n",
    "    \"\"\"Read some information from the image header and write into the df row.\n",
    "    \"\"\"\n",
    "    hdu = 0\n",
    "\n",
    "    # Note: The next line usually works, but fitsio doesn't support CONTINUE lines, which DES\n",
    "    #       image headers sometimes include.\n",
    "    #h = fitsio.read_header(img_file, hdu)\n",
    "    # I don't  care about any of the lines the sometimes use CONITNUE (e.g. OBSERVER), so I\n",
    "    # just remove them and make the header with the rest of the entries.\n",
    "    f = fitsio.FITS(img_file)\n",
    "    header_list = f[hdu].read_header_list()\n",
    "    header_list = [ d for d in header_list if 'CONTINUE' not in d['name'] ]\n",
    "    h = fitsio.FITSHDR(header_list)\n",
    "    try:\n",
    "        date = h['DATE-OBS']\n",
    "        date, time = date.strip().split('T',1)\n",
    "\n",
    "        filter = h['FILTER']\n",
    "        filter = filter.split()[0]\n",
    "\n",
    "        sat = h['SATURATE']\n",
    "        fwhm = h['FWHM']\n",
    "\n",
    "        ccdnum = int(h['CCDNUM'])\n",
    "        detpos = h['DETPOS'].strip()\n",
    "\n",
    "        telra = h['TELRA']\n",
    "        teldec = h['TELDEC']\n",
    "        telha = h['HA']\n",
    "        if galsim.__version__ >= '1.5.1':\n",
    "            telra = galsim.Angle.from_hms(telra) / galsim.degrees\n",
    "            teldec = galsim.Angle.from_dms(teldec) / galsim.degrees\n",
    "            telha = galsim.Angle.from_hms(telha) / galsim.degrees\n",
    "        else:\n",
    "            telra = galsim.HMS_Angle(telra) / galsim.degrees\n",
    "            teldec = galsim.DMS_Angle(teldec) / galsim.degrees\n",
    "            telha = galsim.HMS_Angle(telha) / galsim.degrees\n",
    "\n",
    "        airmass = float(h.get('AIRMASS',-999))\n",
    "        sky = float(h.get('SKYBRITE',-999))\n",
    "        sigsky = float(h.get('SKYSIGMA',-999))\n",
    "\n",
    "        tiling = int(h.get('TILING',0))\n",
    "        hex = int(h.get('HEX',0))\n",
    "\n",
    "    except Exception as e:\n",
    "        logger.info(\"Caught %s\",e)\n",
    "        logger.info(\"Cannot read header information from %s\", img_file)\n",
    "        raise\n",
    "\n",
    "    row['date'] = date\n",
    "    row['time'] = time\n",
    "    row['sat'] = sat\n",
    "    row['fits_filter'] = filter\n",
    "    row['fits_fwhm'] = fwhm\n",
    "    row['fits_ccdnum'] = ccdnum\n",
    "    row['telra'] = telra\n",
    "    row['teldec'] = teldec\n",
    "    row['telha'] = telha\n",
    "    row['airmass'] = airmass\n",
    "    row['sky'] = sky\n",
    "    row['sigsky'] = sigsky\n",
    "    row['tiling'] = tiling\n",
    "    row['hex'] = hex\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "usage: ipykernel_launcher.py [-h] [--sex_dir SEX_DIR] [--piff_exe PIFF_EXE]\n",
      "                             [--findstars_dir FINDSTARS_DIR] [--work WORK]\n",
      "                             [--scratch SCRATCH] [--pixmappy_dir PIXMAPPY_DIR]\n",
      "                             [--tag TAG] [--base_exposures BASE_EXPOSURES]\n",
      "                             [--file FILE] [--exps EXPS [EXPS ...]]\n",
      "                             [--sex_config SEX_CONFIG]\n",
      "                             [--piff_config PIFF_CONFIG]\n",
      "                             [--findstars_config FINDSTARS_CONFIG]\n",
      "                             [--sex_params SEX_PARAMS]\n",
      "                             [--sex_filter SEX_FILTER] [--sex_nnw SEX_NNW]\n",
      "                             [--tapebump_file TAPEBUMP_FILE]\n",
      "                             [--make_symlinks MAKE_SYMLINKS] [--noweight]\n",
      "                             [--clear_output CLEAR_OUTPUT]\n",
      "                             [--rm_files RM_FILES]\n",
      "                             [--use_existing USE_EXISTING]\n",
      "                             [--blacklist BLACKLIST] [--run_piff RUN_PIFF]\n",
      "                             [--run_sextractor RUN_SEXTRACTOR]\n",
      "                             [--run_findstars RUN_FINDSTARS]\n",
      "                             [--mag_cut MAG_CUT] [--min_mag MIN_MAG]\n",
      "                             [--nbright_stars NBRIGHT_STARS]\n",
      "                             [--max_mag MAX_MAG]\n",
      "                             [--use_tapebumps USE_TAPEBUMPS]\n",
      "                             [--tapebump_extra TAPEBUMP_EXTRA]\n",
      "                             [--single_ccd SINGLE_CCD] [--reserve RESERVE]\n",
      "                             [--get_psfex] [--plot_fs] [--use_ngmix]\n",
      "ipykernel_launcher.py: error: unrecognized arguments: -f /global/u2/a/aamon/.local/share/jupyter/runtime/kernel-eebb527c-9a9d-4dfa-8e39-be80983d5ebe.json\n"
     ]
    },
    {
     "ename": "SystemExit",
     "evalue": "2",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mSystemExit\u001b[0m\u001b[0;31m:\u001b[0m 2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/global/homes/a/aamon/.conda/envs/DESenv/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2886: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.\n",
      "  warn(\"To exit: use 'exit', 'quit', or Ctrl-D.\", stacklevel=1)\n"
     ]
    }
   ],
   "source": [
    "# Make the work directory if it does not exist yet.\n",
    "#This chunk of code doens't work\n",
    "args = parse_args()\n",
    "work = os.path.expanduser(args.work)\n",
    "try:\n",
    "    if not os.path.exists(work):\n",
    "        os.makedirs(work)\n",
    "except OSError:\n",
    "    if not os.path.exists(work): raise\n",
    "scratch = os.path.expanduser(args.scratch)\n",
    "logger.info('scratch dir = %s',scratch)\n",
    "try:\n",
    "    if not os.path.exists(scratch):\n",
    "        os.makedirs(scratch)\n",
    "\n",
    "except OSError:\n",
    "    if not os.path.exists(scratch): raise\n",
    "            \n",
    "# A listing Erin made of all the exposures in Y3 used in meds files\n",
    "all_exp = fitsio.read(args.base_exposures)\n",
    "# Switch to native endians, so pandas doesn't complain.\n",
    "all_exp = all_exp.astype(all_exp.dtype.newbyteorder('='))\n",
    "            \n",
    "row = pandas.DataFrame(info).iloc[0]\n",
    "read_image_header(row, image_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#put the stars data into a dataframe \n",
    "\n",
    "def read_findstars(star_file, img_file):\n",
    "    \"\"\"Read the findstars output file\n",
    "    \"\"\"\n",
    "    if not os.path.exists(star_file):\n",
    "        return None\n",
    "\n",
    "    # Read the output and make a DataFrome with the contents *********something buggy here\n",
    "    data = fitsio.read(star_file)\n",
    "    data = data.astype(data.dtype.newbyteorder('='))\n",
    "    print(data) \n",
    "    df = pandas.DataFrame(data)\n",
    "    print(df)\n",
    "    ntot = len(df)\n",
    "    ######nstars = df['star_flag'].sum()\n",
    "    #logger.info('   found %d stars',ntot)\n",
    "    print('   found %d stars',ntot)\n",
    "\n",
    "    #print('mag range = ',np.min(df['mag']), np.max(df['mag']))\n",
    "    #####is_star = df['star_flag'] == 1\n",
    "    #print('star mag range = ',np.min(df['mag'][is_star]), np.max(df['mag'][is_star]))\n",
    "    #print('zero point = ',magzp)\n",
    "    #####df['mag'] += magzp - 25.\n",
    "    #print('star mag range => ',np.min(df['mag'][is_star]), np.max(df['mag'][is_star]))\n",
    "\n",
    "    #Add on some extra information from the sextractor catalog\n",
    "    #INSTEAD I'LL USE THE WCS AND THE X,Y TO GET RA AND DEC.\n",
    "    image = galsim.fits.read(img_file)\n",
    "    wcs = image.wcs\n",
    "    w = wcs.WCS(hdulist[0].header)\n",
    "    world = w.wcs_pix2world((x,y))\n",
    "    print(world)\n",
    "    df['ra'] = world[:,0]\n",
    "    df['dec'] = world[:,1]\n",
    "    print(df)\n",
    "    return df\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ (['SIMPLE  =                    T / This is a FITS file                            ', 'BITPIX  =                    8 /                                                ', 'NAXIS   =                    2 / 2D data                                        ', 'NAXIS1  =                    1 / Number of rows                                 ', 'NAXIS2  =                    1 / Number of columns                              ', 'EXTEND  =                    T / This file may contain FITS extensions          ', 'END                                                                             '],)]\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "If using all scalar values, you must pass an index",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mValueError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-585a41c7b529>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdf\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0mread_findstars\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstar_file\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mimage_file\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-8-cba47c209a20>\u001b[0m in \u001b[0;36mread_findstars\u001b[0;34m(star_file, img_file)\u001b[0m\n\u001b[1;32m     11\u001b[0m     \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnewbyteorder\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'='\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m     \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m     \u001b[0mdf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpandas\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     14\u001b[0m     \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m     \u001b[0mntot\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/global/homes/a/aamon/.conda/envs/DESenv/lib/python2.7/site-packages/pandas/core/frame.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, data, index, columns, dtype, copy)\u001b[0m\n\u001b[1;32m    371\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mcolumns\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    372\u001b[0m                     \u001b[0mcolumns\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata_columns\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 373\u001b[0;31m                 \u001b[0mmgr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_init_dict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    374\u001b[0m             \u001b[0;32melif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'name'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    375\u001b[0m                 mgr = self._init_dict({data.name: data}, index, columns,\n",
      "\u001b[0;32m/global/homes/a/aamon/.conda/envs/DESenv/lib/python2.7/site-packages/pandas/core/frame.pyc\u001b[0m in \u001b[0;36m_init_dict\u001b[0;34m(self, data, index, columns, dtype)\u001b[0m\n\u001b[1;32m    437\u001b[0m                 \u001b[0;31m# GH10856\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    438\u001b[0m                 \u001b[0;31m# raise ValueError if only scalars in dict\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 439\u001b[0;31m                 \u001b[0mindex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mextract_index\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marrays\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m~\u001b[0m\u001b[0mmissing\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    440\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    441\u001b[0m                 \u001b[0mindex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_ensure_index\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/global/homes/a/aamon/.conda/envs/DESenv/lib/python2.7/site-packages/pandas/core/frame.pyc\u001b[0m in \u001b[0;36mextract_index\u001b[0;34m(data)\u001b[0m\n\u001b[1;32m   7391\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   7392\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mindexes\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mraw_lengths\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 7393\u001b[0;31m             raise ValueError('If using all scalar values, you must pass'\n\u001b[0m\u001b[1;32m   7394\u001b[0m                              ' an index')\n\u001b[1;32m   7395\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: If using all scalar values, you must pass an index"
     ]
    }
   ],
   "source": [
    "df= read_findstars(star_file,image_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read in psf model file\n",
    "\n",
    "psfex_file= \"/global/homes/a/aamon/DES/DEStests/deeppsfs/UltraVista/UVISTA_J_21_01_16_psfcat.psf\"\n",
    "dat = fits.open(psf_file)\n",
    "print(dat.info()) \n",
    "print(dat[1].header)\n",
    "data= dat[1].data\n",
    "\n",
    "#if args.get_psfex:\n",
    "#  if not (args.use_existing and os.path.exists(psfex_File)):\n",
    "#   psfex_file = wget(url_base, base_path + '/psf/', wdir, root + '_psfexcat.psf', logger)\n",
    "#    logger.info('psfex_file = %s',psfex_file)\n",
    "#    row['psfex_file'] = psfex_file\n",
    "#    keep_files.append(psfex_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#neither my starlist nor Mike's has an obs_flux in the starlist?\n",
    "\n",
    "def measure_psfex_shapes(df, psfex_file, image_file, noweight, wcs, fwhm): #, logger):\n",
    "    \"\"\"Measure shapes of the PSFEx solution at each location.\n",
    "    \"\"\"\n",
    "    #logger.info('Read in PSFEx file: %s',psfex_file)\n",
    "\n",
    "    #ignore fact that I have no star_file for now\n",
    "    ind = df.index[df] \n",
    "    #ind = df.index[df['star_flag'] == 1]\n",
    "    #logger.info('ind = %s',ind)\n",
    "    #n_psf = len(ind)\n",
    "    #logger.info('n_psf = %s',n_psf)\n",
    "\n",
    "    df['psfex_dx'] = [ -999. ] * len(df)\n",
    "    df['psfex_dy'] = [ -999. ] * len(df)\n",
    "    df['psfex_e1'] = [ -999. ] * len(df)\n",
    "    df['psfex_e2'] = [ -999. ] * len(df)\n",
    "    df['psfex_T'] = [ -999. ] * len(df)\n",
    "    df['psfex_flux'] = [ -999. ] * len(df)\n",
    "    df['psfex_flag'] = [ NOT_STAR ] * len(df)\n",
    "    df.loc[ind, 'psfex_flag'] = 0\n",
    "\n",
    "    full_image = galsim.fits.read(image_file, hdu=0)\n",
    "\n",
    "    if wcs is not None:\n",
    "        full_image.wcs = wcs\n",
    "\n",
    "    if not noweight:\n",
    "        print(\"I'm using a weight)\")\n",
    "        full_weight = galsim.fits.read(image_file, hdu=0)\n",
    "        full_weight.array[full_weight.array < 0] = 0.\n",
    "\n",
    "    stamp_size = 48\n",
    "\n",
    "    for i in ind:\n",
    "        x = df['X_IMAGE'].iloc[i]\n",
    "        y = df['Y_IMAGE'].iloc[i]\n",
    "        \n",
    "        #print('Measure PSFEx model shape at ',x,y)\n",
    "        image_pos = galsim.PositionD(x,y)\n",
    "        psf_i = psf.getPSF(image_pos)\n",
    "\n",
    "        b = galsim.BoundsI(int(x)-stamp_size/2, int(x)+stamp_size/2,\n",
    "                           int(y)-stamp_size/2, int(y)+stamp_size/2)\n",
    "        b = b & full_image.bounds\n",
    "        im = full_image[b]\n",
    "\n",
    "        im = psf_i.drawImage(image=im, method='no_pixel')\n",
    "        im *= df['obs_flux'].iloc[i]    \n",
    "\n",
    "        if noweight:\n",
    "            wt = None\n",
    "        else:\n",
    "            wt = full_weight[b]\n",
    "            var = wt.copy()\n",
    "            var.invertSelf()\n",
    "            im.addNoise(galsim.VariableGaussianNoise(rng, var))\n",
    "        \n",
    "        dx, dy, e1, e2, T, flux, flag = ngmix_fit(im, wt, fwhm, x, y, logger)\n",
    "        \n",
    "        if np.any(np.isnan([dx,dy,e1,e2,T,flux])):\n",
    "            logger.info(' *** NaN detected (%f,%f,%f,%f,%f,%f).',dx,dy,e1,e2,T,flux)\n",
    "            flag |= BAD_MEASUREMENT\n",
    "        else:\n",
    "            df.loc[i, 'psfex_dx'] = dx\n",
    "            df.loc[i, 'psfex_dy'] = dy\n",
    "            df.loc[i, 'psfex_e1'] = e1\n",
    "            df.loc[i, 'psfex_e2'] = e2\n",
    "            df.loc[i, 'psfex_T'] = T\n",
    "            df.loc[i, 'psfex_flux'] = flux\n",
    "        df.loc[i, 'psfex_flag'] |= flag\n",
    "    logger.info('final psfex_flag = %s',df['psfex_flag'][ind].values)\n",
    "    #print('df[ind] = ',df.loc[ind].describe())\n",
    "    flag_outliers(df, ind, 'psfex', 4., logger)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measure_psfex_shapes(df, psfex_file, image_file, noweight=False, wcs, fwhm) #, logger)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure_star_shapes(df, image_file, noweight, wcs, fwhm): #, logger):\n",
    "    \"\"\"Measure shapes of the raw stellar images at each location.\n",
    "    \"\"\"\n",
    "    #logger.info('Read in stars in file: %s',image_file)\n",
    "\n",
    "    ind = df.index[df] #['star_flag'] == 1]\n",
    "    #logger.info('ind = %s',ind)\n",
    "    n_psf = len(ind)\n",
    "    #logger.info('n_psf = %s',n_psf) #ignore logger for now\n",
    "    print('n_psf = %s',n_psf)\n",
    "\n",
    "    df['obs_dx'] = [ -999. ] * len(df)\n",
    "    df['obs_dy'] = [ -999. ] * len(df)\n",
    "    df['obs_e1'] = [ -999. ] * len(df)\n",
    "    df['obs_e2'] = [ -999. ] * len(df)\n",
    "    df['obs_T'] = [ -999. ] * len(df)\n",
    "    df['obs_flux'] = [ -999. ] * len(df)\n",
    "    df['obs_flag'] = [ NOT_STAR ] * len(df)\n",
    "    df.loc[ind, 'obs_flag'] = 0\n",
    "\n",
    "    full_image = galsim.fits.read(image_file, hdu=0)\n",
    "\n",
    "    if wcs is not None:\n",
    "        full_image.wcs = wcs\n",
    "\n",
    "    if not noweight:\n",
    "        full_weight = galsim.fits.read(image_file, hdu=2)\n",
    "        full_weight.array[full_weight.array < 0] = 0.\n",
    "\n",
    "    stamp_size = 48\n",
    "\n",
    "    for i in ind:\n",
    "        x = df['x'].iloc[i]\n",
    "        y = df['y'].iloc[i]\n",
    "\n",
    "        #print('Measure shape for star at ',x,y)\n",
    "        b = galsim.BoundsI(int(x)-stamp_size/2, int(x)+stamp_size/2,\n",
    "                           int(y)-stamp_size/2, int(y)+stamp_size/2)\n",
    "        b = b & full_image.bounds\n",
    "        im = full_image[b]\n",
    "\n",
    "        if noweight:\n",
    "            wt = None\n",
    "        else:\n",
    "            wt = full_weight[b]\n",
    "\n",
    "        \n",
    "        dx, dy, e1, e2, T, flux, flag = ngmix_fit(im, wt, fwhm, x, y, logger)\n",
    "\n",
    "        #logger.info('ngmix measurement: (%f,%f,%f,%f,%f,%f).',dx,dy,e1,e2,T,flux)\n",
    "        if np.any(np.isnan([dx,dy,e1,e2,T,flux])):\n",
    "            logger.info(' *** NaN detected (%f,%f,%f,%f,%f,%f).',dx,dy,e1,e2,T,flux)\n",
    "            flag |= BAD_MEASUREMENT\n",
    "        else:\n",
    "            df.loc[i, 'obs_dx'] = dx\n",
    "            df.loc[i, 'obs_dy'] = dy\n",
    "            df.loc[i, 'obs_e1'] = e1\n",
    "            df.loc[i, 'obs_e2'] = e2\n",
    "            df.loc[i, 'obs_T'] = T\n",
    "            df.loc[i, 'obs_flux'] = flux\n",
    "        df.loc[i, 'obs_flag'] |= flag\n",
    "    logger.info('final obs_flag = %s',df['obs_flag'][ind].values)\n",
    "    #print('df[ind] = ',df.loc[ind].describe())\n",
    "    flag_outliers(df, ind, 'obs', 4., logger)\n",
    "\n",
    "    # Any stars that weren't measurable here, don't use for PSF fitting.\n",
    "    df.loc[df['obs_flag']!=0, 'use'] = False\n",
    "    \n",
    "measure_star_shapes(df, image_file, noweight, wcs, fwhm, logger)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_cat_file = os.path.join(wdir, 'exp_psf_cat_%d.fits'%exp)\n",
    "        with fitsio.FITS(exp_cat_file,'rw',clobber=True) as f:\n",
    "            f.write_table(exp_stars_df.to_records(index=False), extname='stars')\n",
    "            f.write_table(exp_info_df.to_records(index=False), extname='info')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desenv",
   "language": "python",
   "name": "desenv"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}