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

from __future__ import print_function, \
    division, \

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

import sep
import photutils

import numpy as np

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.gridspec as gridspec
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter, MaxNLocator, FormatStrFormatter
from matplotlib.collections import PatchCollection
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
plt.rc('text', usetex=True)

from astropy import wcs
from import fits
from astropy.visualization import ZScaleInterval, \
    PercentileInterval, \
from astropy.convolution import Gaussian2DKernel
from astropy.convolution import convolve

# Increase the pixel stack for large image

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

In [40]:
# Some useful functions from my code

from kungpao import io
from kungpao.galsbp import galSBP
from kungpao import utils
from kungpao import detection
from kungpao import imtools
from kungpao.display import display_single, IMG_CMAP, SEG_CMAP

In [4]:
# Convolution kernel I used for detecting objects on the image

kernel3 = np.asarray([[0.092163, 0.221178, 0.296069, 0.221178, 0.092163],
                      [0.221178, 0.530797, 0.710525, 0.530797, 0.221178],
                      [0.296069, 0.710525, 0.951108, 0.710525, 0.296069],
                      [0.221178, 0.530797, 0.710525, 0.530797, 0.221178],
                      [0.092163, 0.221178, 0.296069, 0.221178, 0.092163]])

# Directions to my IRAF executbles
ISO = '/Users/song/iraf/extern/stsdas/bin.macosx/x_isophote.e'
TBL = '/Users/song/iraf/extern/tables/bin.macosx/x_ttools.e'

Inspect the image

In [23]:
img_dir = '/Users/song/Dropbox/final-images-1/'

# One of my favorite galaxies of all time
hdr_test =, 'M104subsky12.fits'))[0].header
img_test =, 'M104subsky12.fits'))[0].data

# For running sep
img_test = img_test.byteswap().newbyteorder()

img_w, img_h = img_test.shape

In [19]:
# I just realize that these images do not have WCS information in the header.  
# It would be useful to have these information available 

wcs_test = wcs.WCS(hdr_test)

In [27]:
_ = display_single(img_test, contrast=0.2, scale_bar=None)

Image reduction

Measure and subtract a background to remove low-frequency structures

  • This will remove part of the extended halo. The purpose is to reduce the difficulty of source detection and deblending

In [30]:
bkg_local = sep.Background(img_test, bw=20, bh=20, fw=4, fh=4)

print("# Mean Sky / RMS Sky = %10.5f / %10.5f" % (bkg_local.globalback, bkg_local.globalrms))

_ = display_single(bkg_local.back(), contrast=0.5, scale_bar=None)

# The RMS map will be used as the image error for object detection 
_ = display_single(bkg_local.rms(), contrast=0.5, scale_bar=None)

# Mean Sky / RMS Sky =    1.72370 /   14.24700

In [31]:
# Subtract the very "local" sky background
img_subbkg = img_test - bkg_local.back()

# As you can see, most of the extended structures have been subtracted off.
_ = display_single(img_subbkg, contrast=0.3)