In [37]:
# Standard python imports
%matplotlib inline
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
from os.path import join
import os
import glob
import inspect
from importlib import reload
import warnings
import re
from collections import OrderedDict
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import time
import itertools
import matplotlib.image as mpimg
import fitsio
import astropy.cosmology
from astropy.visualization import astropy_mpl_style
plt.style.use(astropy_mpl_style)
from astropy.io import fits
from astropy.visualization import make_lupton_rgb
In [2]:
basedir = "/home/s1/ezaborow/projects/strong_lensing/analysis/schechter_lenses"
In [10]:
coordfile = os.path.join(basedir, "coords.txt")
coords = pd.read_csv(coordfile, names=["name", "ra", "dec"])
coords
Out[10]:
In [86]:
bands = ['g', 'r', 'i', 'z']
for j in np.arange(len(coords)):
# Load image data
g = None
r = None
i = None
z = None
imfile = coords.loc[j, "name"] + "_cutouts.fits"
filepath = os.path.join(basedir, "images", imfile)
with fits.open(filepath) as hdul:
# hdul.info()
# hdul[0].header
imbands = []
imdata = []
for k in np.arange(len(hdul)):
header = hdul[k].header
if (header["TYPE"] == "IMAGE") and (header["BAND"] == "g"):
g = hdul[k].data
if (header["TYPE"] == "IMAGE") and (header["BAND"] == "r"):
r = hdul[k].data
if (header["TYPE"] == "IMAGE") and (header["BAND"] == "i"):
i = hdul[k].data
if (header["TYPE"] == "IMAGE") and (header["BAND"] == "z"):
z = hdul[k].data
if g is not None:
imbands.append('g')
imdata.append(g)
if r is not None:
imbands.append('r')
imdata.append(r)
if i is not None:
imbands.append('i')
imdata.append(i)
if z is not None:
imbands.append('z')
imdata.append(z)
# Plot images
print(coords.loc[j, "name"])
print("RA: {}, DEC: {}".format(coords.loc[j, "ra"], coords.loc[j, "dec"]))
print("BANDS: {}".format(imbands))
if len(imbands) == 1:
fig, axs = plt.subplots(1, len(imbands) + 1, figsize=(5, 5))
else:
fig, axs = plt.subplots(1, len(imbands) + 1, figsize=(20, 20))
for m, band in enumerate(imbands):
_ = axs[m].imshow(imdata[m])
try:
outpath = os.path.join(basedir, "images", coords.loc[j, "name"] + ".png")
rgb = make_lupton_rgb(image_r=imdata[2],
image_g=imdata[1],
image_b=imdata[0],
Q=10,
stretch=1000,
filename=outpath
)
_ = axs[m+1].imshow(rgb)
except:
_ = axs[m+1].axis("off")
_ = plt.show()
In [ ]: