In [1]:
%matplotlib inline 
%load_ext autoreload 
%autoreload 2

#---------------------------------------------------------------------------#
__author__ = 'Song Huang'
__email__ = 'shuang89@ucsc.edu'
__version__ = '170901A'

from __future__ import print_function, \
    division, \
    absolute_import

#---------------------------------------------------------------------------#
import os
import sys
import math
import glob
import copy
import warnings
import subprocess

import sep
import photutils

import numpy as np
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Numba related
try:
    import numba
    from numba import jit
    use_numba = True
except Exception: 
    use_numba = False
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Matplotlib related
import matplotlib as mpl
import matplotlib.pyplot as plt
# To avoid conflict with isophote.Ellipse
from matplotlib.patches import Ellipse as mpl_ellip
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
plt.rc('text', usetex=True)

# Astropy related
from astropy import wcs
from astropy.io import fits
from astropy.table import \
    Table, \
    Column, \
    vstack, \
    unique
from astropy import units as u
from astropy.units import Quantity
from astropy.coordinates import SkyCoord
from astropy.stats import sigma_clip
from astropy.nddata import Cutout2D
from astropy.stats import sigma_clip
from astropy.convolution import convolve, Box1DKernel
from astropy.visualization import ZScaleInterval, \
    PercentileInterval, \
    AsymmetricPercentileInterval

# Scipy related
from scipy.interpolate import LSQUnivariateSpline
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# About the Colormaps
IMG_CMAP = plt.get_cmap('viridis')
IMG_CMAP.set_bad(color='black')

SEG_CMAP = photutils.utils.random_cmap(ncolors=512, background_color=u'white')
SEG_CMAP.set_bad(color='white')
SEG_CMAP.set_under(color='white')
#---------------------------------------------------------------------------#

In [2]:
#---------------------------------------------------------------------------#
# User imports
from kungpao.isophote.ellipse import Ellipse
from kungpao.isophote.ellipse import Centerer
from kungpao.isophote.ellipse import Geometry
from kungpao.isophote.ellipse.model import build_model

from kungpao.galsbp import galSBP

import hsc_massive
from hsc_massive import \
    s16a_path, \
    sample_selection, \
    prepare_sed, \
    catalog_summary, \
    smhm, \
    plotting
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
#envir = s16a_path.set_env(version='astro4')
envir = s16a_path.set_env(version='kungpao')

ORG = plotting.ORG
BLK = plotting.BLK
BLU = plotting.BLU
GRN = plotting.GRN

# For Kungpao
x_images = '/Users/song/iraf/bin.macosx/x_images.e'
x_ttools = '/Users/song/iraf/extern/tables/bin.macosx/x_ttools.e'
x_isophote = '/Users/song/iraf/extern/stsdas/bin.macosx/x_isophote.e'
#---------------------------------------------------------------------------#

In [3]:
from pyraf import iraf

iraf.tables()
iraf.stsdas()
iraf.analysis()
iraf.isophote()

iraf.unlearn('ellipse')
iraf.unlearn('bmodel')



      +------------------------------------------------------------+
      |             Space Telescope Tables Package                 |    
      |                  TABLES Version 3.17                       |
      |                                                            |
      |   Space Telescope Science Institute, Baltimore, Maryland   |
      |   Copyright (C) 2014 Association of Universities for       |
      |            Research in Astronomy, Inc.(AURA)               |
      |       See stsdas$copyright.stsdas for terms of use.        |
      |         For help, send e-mail to help@stsci.edu            |
      +------------------------------------------------------------+
tables/:
 fitsio/        tbplot/         tobsolete/      ttools/


      +------------------------------------------------------------+
      |       Space Telescope Science Data Analysis System         |    
      |                   STSDAS Version 3.17                      |
      |                                                            |
      |   Space Telescope Science Institute, Baltimore, Maryland   |
      |   Copyright (C) 2014 Association of Universities for       |
      |            Research in Astronomy, Inc.(AURA)               |
      |       See stsdas$copyright.stsdas for terms of use.        |
      |         For help, send e-mail to help@stsci.edu            |
      |                                                            |
      +------------------------------------------------------------+
stsdas/:
 analysis/      examples        hst_calib/      sobsolete/
 contrib/       fitsio/         playpen/        toolbox/
 describe       graphics/       problems
isophote/:
 bmodel         geompar@        isomap          magpar@
 controlpar@    isoexam         isopall         samplepar@
 ellipse        isoimap         isoplot

User defined functions


In [4]:
def display_single(img, pixel_scale=0.168,
                   xsize=8, ysize=8, ax=None,
                   stretch='arcsinh', scale='zscale', 
                   contrast=0.25, no_negative=False,
                   lower_percentile=1.0, upper_percentile=99.0,
                   cmap=IMG_CMAP, scale_bar=True,
                   scale_bar_length=20.0, scale_bar_fontsize=20,
                   scale_bar_y_offset=0.5
                   ):
    """
    Display a single image. 
    
    :param img: np 2-D array for image 
    
    :param xsize: int, default = 8
        Width of the image. 
    :param ysize: int, default = 8
        Height of the image. 
    """
    #---------------------------------------------------------------------------#
    if ax is None:
        fig = plt.figure(figsize=(xsize, ysize))
        ax1 = fig.add_subplot(111)
    else: 
        ax1 = ax
    #---------------------------------------------------------------------------#

    #---------------------------------------------------------------------------#
    # Stretch option
    if stretch.strip() is 'arcsinh':
        img_scale = np.arcsinh(img)
    elif stretch.strip() is 'log':
        if no_negative:
            img[img <= 0.0] = 1.0E-10
        img_scale = np.log(img)
    elif stretch.strip() is 'log10':
        if no_negative:
            img[img <= 0.0] = 1.0E-10
        img_scale = np.log10(img)
    elif stretch.strip() is 'linear':
        img_scale = img
    else: 
        raise Exception("# Wrong stretch option.")

    # Scale option
    if scale.strip() is 'zscale':
        zmin, zmax = ZScaleInterval(contrast=contrast).get_limits(img_scale)
    elif scale.strip() is 'percentile':
        zmin, zmax = AsymmetricPercentileInterval(
            lower_percentile=lower_percentile,
            upper_percentile=upper_percentile).get_limits(img_scale)

    ax1.imshow(img_scale, origin='lower', cmap=cmap,
               vmin=zmin, vmax=zmax)
    #---------------------------------------------------------------------------#

    #---------------------------------------------------------------------------#
    # Hide ticks and tick labels
    ax1.tick_params(labelbottom='off', labelleft='off', 
                    axis=u'both', which=u'both', length=0)    
    #---------------------------------------------------------------------------#

    #---------------------------------------------------------------------------#
    # Put scale bar on the image
    (img_size_x, img_size_y) = img.shape
    if scale_bar:
        scale_bar_x_0 = int(img_size_x * 0.95 - 
                            (scale_bar_length / pixel_scale))
        scale_bar_x_1 = int(img_size_x * 0.95)
        scale_bar_y = int(img_size_y * 0.10)
        scale_bar_text_x = (scale_bar_x_0 + scale_bar_x_1) / 2
        scale_bar_text_y = (scale_bar_y * scale_bar_y_offset)
        scale_bar_text = r'$%d^{\prime\prime}$' % int(scale_bar_length)
        scale_bar_text_size = scale_bar_fontsize

        ax1.plot([scale_bar_x_0, scale_bar_x_1], 
                 [scale_bar_y, scale_bar_y], 
                 linewidth=3, c='w', alpha=1.0)
        ax1.text(scale_bar_text_x, scale_bar_text_y, scale_bar_text, 
                 fontsize=scale_bar_text_size,
                 horizontalalignment='center', 
                 color='w')
    #---------------------------------------------------------------------------#
    
    #---------------------------------------------------------------------------#
    if ax is None:
        return fig
    else: 
        return ax1

Load Data


In [41]:
#---------------------------------------------------------------------------#
# Read the images and get the header (WCS)
gal_img = fits.open('redadd_529_HSC-I_full_img.fits')[0].data
gal_hdr = fits.open('redadd_529_HSC-I_full_img.fits')[0].header
gal_wcs = wcs.WCS(gal_hdr)

gal_bad = fits.open('redadd_529_HSC-I_full_bad.fits')[0].data
gal_sig = fits.open('redadd_529_HSC-I_full_sig.fits')[0].data
gal_psf = fits.open('redadd_529_HSC-I_full_psf.fits')[0].data

# The pixel values on PSF images are often too low, need to scale to a 
# larger value
psf_flux_scale = 3000.0
exposure_secs = (20.0 * 60.0)
gal_psf *= psf_flux_scale
gal_img_expcor = gal_img * exposure_secs 

gal_img_swap = gal_img.byteswap().newbyteorder()
gal_sig_swap = gal_sig.byteswap().newbyteorder()
gal_bad_swap = gal_bad.byteswap().newbyteorder()
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Basic information for HSC
hsc_pixel_scale = 0.168 # arcsec/ pixel
hsc_zeropoint = 27.0    # mag

print("# 10 arcsecs == %d pixels" % (10.0 / hsc_pixel_scale))
hsc_pixels_20asec = (20.0 / hsc_pixel_scale)
hsc_pixels_10asec = (10.0 / hsc_pixel_scale)
hsc_pixels_5asec = (5.0 / hsc_pixel_scale)
hsc_pixels_3asec = (3.0 / hsc_pixel_scale)

print("# 0.01/pixel == %6.3f mag/arcsec^2" % (
    -2.5 * np.log10(0.01 / (hsc_pixel_scale ** 2.0)) + hsc_zeropoint)
     )
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Basic information about the image
(img_size_x, img_size_y) = gal_img.shape
x_cen_pix, y_cen_pix = (img_size_x / 2), (img_size_y / 2)

(psf_size_x, psf_size_y) = gal_psf.shape
x_cen_psf, y_cen_psf = (psf_size_x / 2), (psf_size_y / 2)

ra_cen, dec_cen = gal_wcs.wcs_pix2world(x_cen_pix, y_cen_pix, 1)
ra_cen, dec_cen = float(ra_cen), float(dec_cen)

img_radii_arcsec = (x_cen_pix * hsc_pixel_scale) * math.sqrt(2.0)

print("# Central coordinate: %9.6f, %9.6f" % (ra_cen, dec_cen))
print("# Searching radius for the image: %6.1f arcsec" % img_radii_arcsec)
#---------------------------------------------------------------------------#


# 10 arcsecs == 59 pixels
# 0.01/pixel == 28.127 mag/arcsec^2
# Central coordinate: 140.147266, -0.117204
# Searching radius for the image:  142.7 arcsec

1-D profile of the PSF model


In [42]:
psf_fig = display_single(gal_psf, xsize=4, ysize=4, ax=None,
                         stretch='log', scale='percentile', no_negative=True,
                         lower_percentile=0.5, upper_percentile=99.0,
                         contrast=0.15, scale_bar=True, scale_bar_y_offset=0.20,
                         scale_bar_length=1, scale_bar_fontsize=15.0)



In [43]:
psf_isos_iraf, psf_bin_iraf = galSBP.galSBP('redadd_529_HSC-I_full_psf.fits', 
                                            galX=20.5, galY=21.5, galQ=0.95, galPA=0.0, 
                                            maxSma=22, iniSma=4.0, stage=1, intMode='bi-linear',
                                            isophote=x_isophote, xttools=x_ttools, ellipStep=0.05,
                                            savePng=False)


----------------------------------------------------------------------------------------------------
###      galX, galY :  20.5 21.5
###      galR :  20.0
###      iniSma, maxSma :  10.0 22
###      Stage :  1
###      Step :  0.05
----------------------------------------------------------------------------------------------------
##       Set up the Ellipse configuration
----------------------------------------------------------------------------------------------------

----------------------------------------------------------------------------------------------------
##       Start the Ellipse Run: Attempt  1
----------------------------------------------------------------------------------------------------
###      Origin Image  : redadd_529_HSC-I_full_psf.fits
###      Input Image   : temp_URTCN.fits
###      Output Binary : redadd_529_HSC-I_full_psf_ellip_1.bin
----------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------
###     Input background value   :  0.0
###     1-D SBP background value :  1.2611639e-06
###     Current outer background :  1.2611639e-06
----------------------------------------------------------------------------------------------------

In [44]:
try: 
    os.remove('redadd_529_HSC-I_full_psf_ellip_1.fits')
except Exception: 
    pass

iraf.bmodel(parent='redadd_529_HSC-I_full_psf.fits', 
            table='redadd_529_HSC-I_full_psf_ellip_1.bin',
            output='redadd_529_HSC-I_full_psf_ellip_1.fits',
            minsma=0.0,
            highar='no')

In [46]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 10))

ax1.imshow(np.arcsinh(gal_psf))
ax1.set_title(r"$\mathrm{Data}$")

psf_iraf = fits.open('redadd_529_HSC-I_full_psf_ellip_1.fits')[0].data
ax2.imshow(np.arcsinh(psf_iraf * psf_flux_scale))
ax2.set_title(r"$\mathrm{Model}$")

ax3.imshow(np.arcsinh(gal_psf - psf_iraf * psf_flux_scale))
ax3.set_title(r"$\mathrm{Residual}$")


Out[46]:
<matplotlib.text.Text at 0x119723950>

In [53]:
psf_geometry = Geometry(x_cen_psf, y_cen_psf, 
                        3.0, 0.01, (90.0 * np.pi / 180.0))
psf_ellip = Ellipse(gal_psf, geometry=psf_geometry, verbose=False)

psf_isos = psf_ellip.fit_image(sma0=4.0, minsma=0.0, maxsma=(x_cen_psf * 0.9), 
                               step=0.15, sclip=2.0, nclip=0, integrmode='bi-linear', 
                               conver=0.02, maxit=50, fflag=0.6, maxgerr=0.5,
                               verbose=True)

psf_model = build_model(gal_psf, psf_isos, high_harmonics=False)


#
# Semi-      Isophote         Ellipticity    Position     Grad.   Data  Flag Iter. Stop
# major        mean                           Angle        rel.                    code
# axis       intensity                                    error
#(pixel)                                     (degree)
#
   4.00       10.87 ( 0.09)  0.010 (0.001) 173.14 ( 3.7)  0.037    26     0   50     2
   4.60        6.84 ( 0.03)  0.010 (0.001) 148.42 ( 3.1)  0.027    29     0   50     2
   5.29        4.17 ( 0.02)  0.006 (0.001) 137.35 ( 5.0)  0.022    34     0   18     0
   6.08        2.52 ( 0.01)  0.010 (0.001)  10.00 ( 4.0)  0.023    39     0   24     0
   7.00        1.54 ( 0.01)  0.011 (0.001)  19.83 ( 2.3)  0.019    44     0   10     0
   8.05        0.94 ( 0.01)  0.011 (0.001)  56.01 ( 3.1)  0.029    51     0   50     2
   9.25        0.58 ( 0.01)  0.024 (0.002)  71.38 ( 2.2)  0.044    58     0   10     0
  10.64        0.34 ( 0.00)  0.005 (0.002) 122.60 ( 9.1)  0.044    67     0   19     0
  12.24        0.21 ( 0.00)  0.028 (0.002)  75.61 ( 2.3)  0.071    76     0   10     0
  14.07        0.12 ( 0.00)  0.045 (0.003)  57.24 ( 1.9)  0.098    87     0   10     0
  16.18        0.08 ( 0.00)  0.088 (0.003)  48.98 ( 1.1)  0.113    98     0   10     0
   3.48       17.21 ( 0.10)  0.026 (0.002) 134.69 ( 2.3)  0.030    22     0   50     2
   3.02       26.04 ( 0.18)  0.031 (0.002) 147.83 ( 2.4)  0.038    19     0   50     2
   2.63       37.85 ( 0.19)  0.032 (0.002) 136.44 ( 2.0)  0.039    17     0   12     0
   2.29       50.61 ( 0.29)  0.030 (0.003) 138.67 ( 3.3)  0.036    15     0   12     0
   1.99       65.98 ( 0.42)  0.030 (0.004) 146.25 ( 3.8)  0.059    13     0   50     2
   1.73       82.24 ( 0.39)  0.027 (0.004) 147.77 ( 4.1)  0.053    13     0   19     0
   1.50       98.15 ( 0.44)  0.034 (0.004) 145.36 ( 4.0)  0.057    13     0   50     2
   1.31      112.28 ( 0.52)  0.038 (0.006) 145.36 ( 4.8)  0.076    13     0   50     2
   1.14      124.23 ( 0.61)  0.037 (0.008) 145.36 ( 6.8)  0.102    13     0   50     2
   0.99      134.24 ( 0.82)  0.027 (0.012) 145.36 (14.3)  0.153    13     0   10     0
   0.86      142.82 ( 1.08)  0.023 (0.020) 144.84 (26.9)  0.249    13     0   18     0
   0.75      151.00 ( 1.20)  0.092 (0.029) 165.21 (10.3)  0.443    13     0   10     0
   0.65      156.03 ( 1.13)  0.086 (0.032) 167.06 (12.0)  0.491    13     0   10     0
   0.57      160.39 ( 1.02)  0.080 (0.033) 167.41 (13.0)  0.493    13     0   10     0
   0.00      182.20
Interpolating....Done
SMA= 16.1
Done

In [54]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 10))

ax1.imshow(np.arcsinh(gal_psf))
ax1.set_title(r"$\mathrm{Data}$")

ax2.imshow(np.arcsinh(psf_model))
ax2.set_title(r"$\mathrm{Model}$")

ax3.imshow(np.arcsinh(gal_psf - psf_model))
ax3.set_title(r"$\mathrm{Residual}$")


Out[54]:
<matplotlib.text.Text at 0x11ff9a410>

In [50]:
plt.scatter(psf_isos.sma ** 0.25, psf_isos.intens, s=15)
plt.scatter(psf_isos_iraf['sma'] ** 0.25, psf_isos_iraf['intens'] * psf_flux_scale, s=6)


Out[50]:
<matplotlib.collections.PathCollection at 0x11ff1cb10>

Search for stars on the image


In [30]:
#---------------------------------------------------------------------------#
# Search for bright stars in Gaia in time
from astroquery.gaia import Gaia
from astropy.coordinates import ICRS, FK5

img_cen_ra_dec = SkyCoord(ra_cen, dec_cen, unit=('deg', 'deg'), 
                          frame='icrs')

img_search_x = Quantity(hsc_pixel_scale * (gal_img.shape)[0], u.arcsec)
img_search_y = Quantity(hsc_pixel_scale * (gal_img.shape)[1], u.arcsec)

gaia_results = Gaia.query_object_async(coordinate=img_cen_ra_dec, 
                                       width=img_search_x, 
                                       height=img_search_y,
                                       verbose=False)

ra_gaia, dec_gaia = np.asarray(gaia_results['ra']), np.asarray(gaia_results['dec'])
rmask_gaia_arcsec = 694.7 * np.exp(-gaia_results['phot_g_mean_mag'] / 4.04)

x_gaia, y_gaia = gal_wcs.wcs_world2pix(ra_gaia, dec_gaia, 1)

gaia_results.add_column(Column(data=x_gaia, name='x_pix'))
gaia_results.add_column(Column(data=y_gaia, name='y_pix'))
gaia_results.add_column(Column(data=rmask_gaia_arcsec, name='rmask_arcsec'))

# Keep the bright ones
gaia_gmag_limit = 17.5
gaia_bright = gaia_results[gaia_results['phot_g_mean_mag'] <= gaia_gmag_limit]

print("# Gaia finds %d stars and %d with g < %4.2f mag" % (len(gaia_results),
                                                          len(gaia_bright),
                                                          gaia_gmag_limit))
#---------------------------------------------------------------------------#


Launched query: 'SELECT DISTANCE(POINT('ICRS',ra,dec),                 POINT('ICRS',140.147272588,-0.117209164971)) AS dist, *                 FROM gaiadr1.gaia_source WHERE CONTAINS(                POINT('ICRS',ra,dec),                BOX('ICRS',140.147272588,-0.117209164971, 0.0587155555556, 0.0587155555556))=1                 ORDER BY dist ASC'
Retrieving async. results...
Query finished.
# Gaia finds 9 stars and 5 with g < 17.50 mag

Step by Step test


In [34]:
fig = plt.figure(figsize=(8, 8))
ax1 = fig.add_subplot(111)

gal_img_zscale = np.arcsinh(gal_img)

zmin, zmax = ZScaleInterval(contrast=0.25).get_limits(gal_img_zscale)

ax1.imshow(gal_img_zscale, origin='lower', cmap=IMG_CMAP,
           vmin=zmin, vmax=zmax)

#---------------------------------------------------------------------------#
# Hide ticks and tick labels
ax1.tick_params(labelbottom='off', labelleft='off', 
                axis=u'both', which=u'both', length=0)    
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Put scale bar on the image
scale_bar_x_0 = int(img_size_x * 0.95 - hsc_pixels_20asec)
scale_bar_x_1 = int(img_size_x * 0.95)
scale_bar_y = int(img_size_y * 0.10)
scale_bar_text_x = (scale_bar_x_0 + scale_bar_x_1) / 2
scale_bar_text_y = (scale_bar_y * 0.50)
scale_bar_text = r'$20^{\prime\prime}$'
scale_bar_text_size = 20

ax1.plot([scale_bar_x_0, scale_bar_x_1], [scale_bar_y, scale_bar_y], 
         linewidth=3, c='w', alpha=1.0)
ax1.text(scale_bar_text_x, scale_bar_text_y, scale_bar_text, 
         fontsize=scale_bar_text_size,
         horizontalalignment='center', 
         color='w')
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Show stars
ax1.scatter(gaia_results['x_pix'], 
            gaia_results['y_pix'], c=ORG(0.8), 
            s=70, alpha=0.7, marker='+')

ax1.scatter(gaia_bright['x_pix'] - 1, 
            gaia_bright['y_pix'], c=ORG(0.9), 
            s=70, alpha=0.8)

# Plot an ellipse for each object
for star in gaia_results:
    smask = mpl_ellip(xy=(star['x_pix'], 
                          star['y_pix']),
                      width=(2.0 * star['rmask_arcsec'] / hsc_pixel_scale),
                      height=(2.0 * star['rmask_arcsec'] / hsc_pixel_scale),
                      angle=0.0)
    smask.set_facecolor(ORG(0.3))
    smask.set_edgecolor(ORG(0.9))
    smask.set_alpha(0.3)
    ax1.add_artist(smask)
#---------------------------------------------------------------------------#



In [28]:
star_use = gaia_bright[0]
star_mod_step = 0.2

star_r_ini, star_e_ini, star_pa_ini = 50.0, 0.00, 0.0
star_r_min, star_r_max = 0.0, (star_use['rmask_arcsec'] / hsc_pixel_scale * 1.1)
star_geom = Geometry(star_use['x_pix'] - 1, 
                     star_use['y_pix'], 
                     star_r_ini, star_e_ini, star_pa_ini)
star_ellip = Ellipse(gal_img, geometry=star_geom, verbose=False)

try:
    star_isos = star_ellip.fit_image(sma0=star_r_ini, 
                                     minsma=star_r_min, 
                                     maxsma=star_r_max, 
                                     step=0.10, linear=False, 
                                     integrmode='median', 
                                     sclip=2.0, nclip=3, 
                                     conver=0.07, maxit=80, fflag=0.8,
                                     fixgeom=False, verbose=True)
except Exception: 
    star_isos = star_ellip.fit_image(sma0=star_r_ini, 
                                     minsma=star_r_min, 
                                     maxsma=star_r_max, 
                                     step=0.10, linear=False, 
                                     integrmode='median', 
                                     sclip=2.0, nclip=2, 
                                     conver=0.07, maxit=80, fflag=0.8,
                                     fixgeom=True, verbose=True)
    
star_mod = build_model(gal_img, star_isos, 
                       high_harmonics=False, 
                       step=star_mod_step, 
                       verbose=True)

plt.scatter(star_isos.sma, star_isos.intens)


#
# Semi-      Isophote         Ellipticity    Position     Grad.   Data  Flag Iter. Stop
# major        mean                           Angle        rel.                    code
# axis       intensity                                    error
#(pixel)                                     (degree)
#
#
# Semi-      Isophote         Ellipticity    Position     Grad.   Data  Flag Iter. Stop
# major        mean                           Angle        rel.                    code
# axis       intensity                                    error
#(pixel)                                     (degree)
#
  50.00        3.13 ( 0.08)  0.000 (0.004)   0.00 ( 0.0)  0.115    85    15    0     4
  55.00        2.36 ( 0.05)  0.000 (0.003)   0.00 ( 0.0)  0.092   103    18    0     4
  60.50        1.77 ( 0.02)  0.000 (0.002)   0.00 ( 0.0)  0.064   106    21    0     4
  66.55        1.30 ( 0.02)  0.000 (0.002)   0.00 ( 0.0)  0.056   111    16    0     4
  73.21        0.95 ( 0.01)  0.000 (0.002)   0.00 ( 0.0)  0.056   110    17    0     4
  80.53        0.71 ( 0.01)  0.000 (0.002)   0.00 ( 0.0)  0.062   112    15    0     4
  88.58        0.54 ( 0.01)  0.000 (0.001)   0.00 ( 0.0)  None    114    13    0     4
  97.44        0.43 ( 0.01)  0.000 (0.000)   0.00 ( 0.0)  None    115    12    0     4
 107.18        0.33 ( 0.01)  0.000 (0.000)   0.00 ( 0.0)  None    115    12    0     4
 117.90        0.25 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None    112    15    0     4
 129.69        0.18 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None    118     9    0     4
 142.66        0.12 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     95    32    0     4
 156.92        0.06 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     93    34    0     4
 172.61        0.04 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     90    37    0     4
 189.87        0.03 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     88    39    0     4
 208.86        0.02 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     83    44    0     4
 229.75        0.01 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     79    48    0     4
 252.72        0.01 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     78    49    0     4
 278.00       -0.00 ( 0.00)  0.000 (0.000)   0.00 ( 0.0)  None     77    50    0     4
  45.45        4.25 ( 0.09)  0.000 (0.004)   0.00 ( 0.0)  0.109    68    15    0     4
  41.32        5.63 ( 0.11)  0.000 (0.004)   0.00 ( 0.0)  0.100    58    11    0     4
  37.57        7.22 ( 0.13)  0.000 (0.006)   0.00 ( 0.0)  0.107    45    12    0     4
  34.15        8.83 ( 0.29)  0.000 (0.013)   0.00 ( 0.0)  0.199    39     8    0     4
  31.05        9.69 ( 0.78)  0.000 (0.022)   0.00 ( 0.0)  0.559    39     0    0     4
  28.22       12.44 ( 1.26)  0.000 (0.023)   0.00 ( 0.0)  0.502    32     0    0     4
  25.66       15.83 ( 1.72)  0.000 (0.024)   0.00 ( 0.0)  0.585    32     0    0     4
  23.33       19.85 ( 2.31)  0.000 (0.021)   0.00 ( 0.0)  0.634    32     0    0     4
  21.20       25.14 ( 3.23)  0.000 (0.025)   0.00 ( 0.0)  0.789    32     0    0     4
  19.28       31.13 ( 4.51)  0.000 (0.031)   0.00 ( 0.0)  0.951    32     0    0     4
  17.52       29.90 ( 4.97)  0.000 (0.112)   0.00 ( 0.0)  5.690    31     1    0     4
  15.93       27.20 ( 5.15)  0.000 (2.322)   0.00 ( 0.0)  None     30     2    0     4
  14.48       29.30 ( 6.01)  0.000 (0.061)   0.00 ( 0.0)  2.456    31     1    0     4
  13.17       25.08 ( 5.58)  0.000 (3.207)   0.00 ( 0.0)  None     32     0    0     4
  11.97        8.25 ( 1.51)  0.000 (0.925)   0.00 ( 0.0)  None     68     9    0     4
  10.88        3.14 ( 0.68)  0.000 (0.614)   0.00 ( 0.0)  None     60    10    0     4
   9.89        0.93 ( 0.16)  0.000 (0.125)   0.00 ( 0.0)  None     54    10    0     4
   8.99        0.91 ( 0.14)  0.000 (0.091)   0.00 ( 0.0)  None     58     0    0     4
   8.18        0.28 ( 0.05)  0.000 (0.052)   0.00 ( 0.0)  None     45     8    0     4
   7.43        0.25 ( 0.03)  0.000 (0.023)   0.00 ( 0.0)  None     46     2    0     4
   6.76        0.21 ( 0.02)  0.000 (0.028)   0.00 ( 0.0)  None     42     2    0     4
   6.14        0.18 ( 0.02)  0.000 (0.037)   0.00 ( 0.0)  None     38     2    0     4
   5.58        0.15 ( 0.02)  0.000 (0.048)   0.00 ( 0.0)  None     34     3    0     4
   5.08        0.11 ( 0.02)  0.000 (0.058)   0.00 ( 0.0)  None     30     3    0     4
   4.61        0.08 ( 0.02)  0.000 (0.076)   0.00 ( 0.0)  None     29     1    0     4
   4.20        0.08 ( 0.02)  0.000 (0.085)   0.00 ( 0.0)  None     28     0    0     4
   3.81        0.08 ( 0.03)  0.000 (0.106)   0.00 ( 0.0)  None     25     0    0     4
   3.47        0.08 ( 0.03)  0.000 (0.122)   0.00 ( 0.0)  None     23     0    0     4
   3.15        0.07 ( 0.03)  0.000 (0.125)   0.00 ( 0.0)  None     21     0    0     4
   2.87        0.05 ( 0.03)  0.000 (0.113)   0.00 ( 0.0)  None     20     0    0     4
   2.60        0.04 ( 0.03)  0.000 (0.110)   0.00 ( 0.0)  None     18     0    0     4
   2.37        0.03 ( 0.03)  0.000 (0.130)   0.00 ( 0.0)  None     16     0    0     4
   2.15        0.01 ( 0.03)  0.000 (0.135)   0.00 ( 0.0)  None     15     0    0     4
   1.96        0.00 ( 0.04)  0.000 (0.130)   0.00 ( 0.0)  None     14     0    0     4
   1.78       -0.01 ( 0.04)  0.000 (0.139)   0.00 ( 0.0)  None     14     0    0     4
   1.62       -0.02 ( 0.04)  0.000 (0.159)   0.00 ( 0.0)  None     14     0    0     4
   1.47       -0.02 ( 0.03)  0.000 (0.183)   0.00 ( 0.0)  None     14     0    0     4
   1.34       -0.03 ( 0.03)  0.000 (0.199)   0.00 ( 0.0)  None     14     0    0     4
   1.22       -0.03 ( 0.02)  0.000 (0.214)   0.00 ( 0.0)  None     14     0    0     4
   1.10       -0.03 ( 0.02)  0.000 (0.228)   0.00 ( 0.0)  None     14     0    0     4
   1.00       -0.02 ( 0.01)  0.000 (0.158)   0.00 ( 0.0)  None     13     1    0     4
   0.91       -0.03 ( 0.01)  0.000 (0.186)   0.00 ( 0.0)  None     12     2    0     4
   0.83       -0.03 ( 0.01)  0.000 (0.098)   0.00 ( 0.0)  1.827    12     2    0     4
   0.75       -0.04 ( 0.01)  0.000 (0.112)   0.00 ( 0.0)  1.775    12     2    0     4
   0.69       -0.04 ( 0.00)  0.000 (0.092)   0.00 ( 0.0)  1.303    12     2    0     4
   0.62       -0.04 ( 0.00)  0.000 (0.080)   0.00 ( 0.0)  1.049    12     2    0     4
   0.57       -0.03 ( 0.00)  0.000 (0.106)   0.00 ( 0.0)  2.078    13     1    0     4
   0.52       -0.04 ( 0.00)  0.000 (0.015)   0.00 ( 0.0)  None     13     1    0     4
   0.00       -0.03
Interpolating....
Done
SMA=277.8
Done
Out[28]:
<matplotlib.collections.PathCollection at 0x117e20bd0>

In [35]:
star_mod = build_model(gal_img, star_isos, 
                       high_harmonics=False, 
                       step=star_mod_step, 
                       verbose=True)


Interpolating....
Done
/Users/song/Dropbox/work/project/kungpao/kungpao/isophote/ellipse/model.py:199: RuntimeWarning: divide by zero encountered in double_scalars
  phi = max((phi + 0.75 / r), PHI_MIN)
SMA=277.8
Done

In [36]:
fig = plt.figure(figsize=(8, 8))
ax1 = fig.add_subplot(111)

gal_img_zscale = np.arcsinh(gal_img - star_mod)

zmin, zmax = ZScaleInterval(contrast=0.25).get_limits(gal_img_zscale)

ax1.imshow(gal_img_zscale, origin='lower', cmap=IMG_CMAP,
           vmin=zmin, vmax=zmax)

#---------------------------------------------------------------------------#
# Hide ticks and tick labels
ax1.tick_params(labelbottom='off', labelleft='off', 
                axis=u'both', which=u'both', length=0)    
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Put scale bar on the image
scale_bar_x_0 = int(img_size_x * 0.95 - hsc_pixels_20asec)
scale_bar_x_1 = int(img_size_x * 0.95)
scale_bar_y = int(img_size_y * 0.10)
scale_bar_text_x = (scale_bar_x_0 + scale_bar_x_1) / 2
scale_bar_text_y = (scale_bar_y * 0.50)
scale_bar_text = r'$20^{\prime\prime}$'
scale_bar_text_size = 20

ax1.plot([scale_bar_x_0, scale_bar_x_1], [scale_bar_y, scale_bar_y], 
         linewidth=3, c='w', alpha=1.0)
ax1.text(scale_bar_text_x, scale_bar_text_y, scale_bar_text, 
         fontsize=scale_bar_text_size,
         horizontalalignment='center', 
         color='w')
#---------------------------------------------------------------------------#

#---------------------------------------------------------------------------#
# Show stars
ax1.scatter(gaia_results['x_pix'], 
            gaia_results['y_pix'], c=ORG(0.8), 
            s=70, alpha=0.7, marker='+')

ax1.scatter(gaia_bright['x_pix'], 
            gaia_bright['y_pix'], c=ORG(0.9), 
            s=70, alpha=0.8)

# Plot an ellipse for each object
for star in gaia_results:
    smask = mpl_ellip(xy=(star['x_pix'], 
                          star['y_pix']),
                      width=(2.0 * star['rmask_arcsec'] / hsc_pixel_scale),
                      height=(2.0 * star['rmask_arcsec'] / hsc_pixel_scale),
                      angle=0.0)
    smask.set_facecolor('none')
    smask.set_edgecolor(ORG(0.9))
    smask.set_alpha(0.3)
    ax1.add_artist(smask)
#---------------------------------------------------------------------------#