Page MenuHomec4science
Files F92758207

mockimgs_sb_display_fits
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 23, 11:49

mockimgs_sb_display_fits
View Options

#!/usr/bin/python3
import numpy as np
from astropy import constants as cte
from astropy import units as u
import argparse
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib as mpl
from pNbody import *
from pNbody import ic
from pNbody.Mockimgs import lib as libMockimgs
import pickle
import os
from astropy.io import fits
####################################################################
# option parser
####################################################################
description="display surface brightness images from a fits file"
epilog ="""
Examples:
--------
mockimgs_sb_display_fits image.fits
mockimgs_sb_display_fits image1.fits image2.fits
mockimgs_sb_display_fits image1.fits image2.fits -o output.png
mockimgs_sb_display_fits image.fits --add_axes --ax_unit kpc --ax_max 3
mockimgs_sb_display_fits image.fits --add_axes --ax_unit arcsec --ax_max 10
mockimgs_sb_display_fits image.fits --add_axes --ax_unit arcmin --ax_max 60
mockimgs_sb_display_fits image.fits --sbmin 25 --sbmax 32
mockimgs_sb_display_fits image.fits --sbcontours 28 30.5
mockimgs_sb_display_fits image.fits --colorbar
"""
parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument(action="store",
dest="files",
metavar='FILE',
type=str,
default=None,
nargs='*',
help='a list of fits files')
parser.add_argument("-o",
action="store",
type=str,
dest="outputfilename",
default=None,
help="Name of the output file")
parser.add_argument("--add_axes",
action="store_true",
default=False,
help='add axes to the figure')
parser.add_argument("--colorbar",
action="store_true",
default=False,
help='add colorbar to the figure')
parser.add_argument('--fig_size',
action="store",
metavar='FLOAT',
type=float,
default=6,
help='size of a single figure in inches')
parser.add_argument('--sbmin',
action="store",
dest="sbmin",
metavar='FLOAT',
type=float,
default=None,
help='surface brightness minimum')
parser.add_argument('--sbmax',
action="store",
dest="sbmax",
metavar='FLOAT',
type=float,
default=None,
help='surface brightness maximum')
parser.add_argument('--sbcontours',
action="store",
dest="sbcontours",
metavar='FLOAT',
type=float,
default=None,
nargs="*",
help='surface brightness contours')
parser.add_argument("--ax_unit",
action="store",
metavar='STRING',
type=str,
default='pixels',
help='axes units (kpc, arcsec, pixels)')
parser.add_argument('--ax_max',
action="store",
metavar='FLOAT',
type=float,
default=20,
help='extention of the image in the axes units')
parser.add_argument("--colormap",
action="store",
default="tab20c",
help='matplotlib colormap name (e.g. tab20c)')
parser.add_argument("--colormap_reverse",
action="store_true",
default=False,
help='reverse colormap')
def mycmap():
import matplotlib.colors as colors
carray = np.array([[0.19215686274509805, 0.5098039215686274, 0.7411764705882353, 1.0],
[0.4196078431372549, 0.6823529411764706, 0.8392156862745098, 1.0],
[0.6196078431372549, 0.792156862745098, 0.8823529411764706, 1.0],
[0.7764705882352941, 0.8588235294117647, 0.9372549019607843, 1.0],
[0.9019607843137255, 0.3333333333333333, 0.050980392156862744, 1.0],
[0.9921568627450981, 0.5529411764705883, 0.23529411764705882, 1.0],
[0.9921568627450981, 0.6823529411764706, 0.4196078431372549, 1.0],
[0.9921568627450981, 0.8156862745098039, 0.6352941176470588, 1.0],
[0.19215686274509805, 0.6392156862745098, 0.32941176470588235, 1.0],
[0.4549019607843137, 0.7686274509803922, 0.4627450980392157, 1.0],
[0.6313725490196078, 0.8509803921568627, 0.6078431372549019, 1.0],
[0.7803921568627451, 0.9137254901960784, 0.7529411764705882, 1.0],
[0.4588235294117647, 0.4196078431372549, 0.6941176470588235, 1.0],
[0.6196078431372549, 0.6039215686274509, 0.7843137254901961, 1.0],
[0.7372549019607844, 0.7411764705882353, 0.8627450980392157, 1.0],
[0.8549019607843137, 0.8549019607843137, 0.9215686274509803, 1.0],
[0.38823529411764707, 0.38823529411764707, 0.38823529411764707, 1.0],
[0.5882352941176471, 0.5882352941176471, 0.5882352941176471, 1.0],
[0.7411764705882353, 0.7411764705882353, 0.7411764705882353, 1.0],
[0.8509803921568627, 0.8509803921568627, 0.8509803921568627, 1.0]]
)
carray = np.array([[0.19215686274509805, 0.6392156862745098, 0.32941176470588235, 1.0],
[0.4549019607843137, 0.7686274509803922, 0.4627450980392157, 1.0],
[0.6313725490196078, 0.8509803921568627, 0.6078431372549019, 1.0],
[0.7803921568627451, 0.9137254901960784, 0.7529411764705882, 1.0],
[0.19215686274509805, 0.5098039215686274, 0.7411764705882353, 1.0],
[0.4196078431372549, 0.6823529411764706, 0.8392156862745098, 1.0],
[0.6196078431372549, 0.792156862745098, 0.8823529411764706, 1.0],
[0.7764705882352941, 0.8588235294117647, 0.9372549019607843, 1.0],
[0.9019607843137255, 0.3333333333333333, 0.050980392156862744, 1.0],
[0.9921568627450981, 0.5529411764705883, 0.23529411764705882, 1.0],
[0.9921568627450981, 0.6823529411764706, 0.4196078431372549, 1.0],
[0.9921568627450981, 0.8156862745098039, 0.6352941176470588, 1.0],
[1.0, 1.0, 1.0, 1.0]]
)
new_cmap = colors.LinearSegmentedColormap.from_list('mymap',carray,N=13)
return new_cmap
def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
import matplotlib.colors as colors
new_cmap = colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),cmap(np.linspace(minval, maxval, n)))
return new_cmap
####################################################################
# main
####################################################################
if __name__ == '__main__':
params = {
"axes.labelsize": 14,
"axes.titlesize": 18,
"font.size": 12,
"legend.fontsize": 12,
"xtick.labelsize": 14,
"ytick.labelsize": 14,
"text.usetex": False,
"figure.subplot.left": 0.0,
"figure.subplot.right": 1.0,
"figure.subplot.bottom": 0.0,
"figure.subplot.top": 1,
"figure.subplot.wspace": 0.02,
"figure.subplot.hspace": 0.02,
"figure.figsize" : (15, 4),
"lines.markersize": 6,
"lines.linewidth": 3.0,
}
plt.rcParams.update(params)
opt = parser.parse_args()
# number of images
n = len(opt.files)
fig, ax = plt.subplots(1,n)
# plot
fig = plt.gcf()
#fig.set_size_inches(opt.fig_size*n,opt.fig_size)
fig.set_size_inches(12,10)
if opt.add_axes:
fig.subplots_adjust(left=0.088)
fig.subplots_adjust(right=0.985)
fig.subplots_adjust(bottom=0.054)
fig.subplots_adjust(top=0.974)
fig.subplots_adjust(wspace=0.0)
fig.subplots_adjust(hspace=0.0)
else:
fig.subplots_adjust(left=0.0)
fig.subplots_adjust(right=0.995)
fig.subplots_adjust(bottom=0.0)
fig.subplots_adjust(top=1)
fig.subplots_adjust(wspace=0.02)
fig.subplots_adjust(hspace=0.02)
#colors2 = plt.cm.gist_heat(np.linspace(0.3, 1.0, 255))
#mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors2)
#mymap = mpl.colormaps["tab20c"]
#cmap = mpl.colormaps[opt.colormap]
# truncate
#cmap = truncate_colormap(cmap, 0.0, 0.59,int(512*0.59))
cmap = mycmap()
if opt.colormap_reverse:
cmap = cmap.reversed()
for i in range(n):
# open image
hdul = fits.open(opt.files[i])
image = hdul[0].data
header = hdul[0].header
# find the image extend
if opt.ax_unit == 'pixels':
NX = header["NX"]
NY = header["NY"]
DX = opt.ax_max
imgext_xmin = 0
imgext_xmax = NX
imgext_ymin = 0
imgext_ymax = NY
xmin = NX//2 - DX
xmax = NX//2 + DX
ymin = NY//2 - DX
ymax = NY//2 + DX
xlabel = r"x [pixel]"
ylabel = r"y [pixel]"
elif opt.ax_unit == 'kpc':
imgext_xmin = header["XMINKPC"]
imgext_xmax = header["XMAXKPC"]
imgext_ymin = header["YMINKPC"]
imgext_ymax = header["YMAXKPC"]
xmin = opt.ax_max
xmax = -opt.ax_max
ymin = opt.ax_max
ymax = -opt.ax_max
xlabel = r"x [kpc]"
ylabel = r"y [kpc]"
elif opt.ax_unit == 'arcsec':
imgext_xmin = header["XMIN"]
imgext_xmax = header["XMAX"]
imgext_ymin = header["YMIN"]
imgext_ymax = header["YMAX"]
xmin = opt.ax_max
xmax = -opt.ax_max
ymin = opt.ax_max
ymax = -opt.ax_max
xlabel = r"x [arcsec]"
ylabel = r"y [arcsec]"
elif opt.ax_unit == 'arcmin':
imgext_xmin = header["XMIN"]/60
imgext_xmax = header["XMAX"]/60
imgext_ymin = header["YMIN"]/60
imgext_ymax = header["YMAX"]/60
xmin = opt.ax_max
xmax = -opt.ax_max
ymin = opt.ax_max
ymax = -opt.ax_max
xlabel = r"x [arcmin]"
ylabel = r"y [arcmin]"
# set ax
if n==1:
axii = ax
else:
axii = ax[i]
axii.set_aspect('equal')
if opt.sbmax==None or opt.sbmin==None:
# remove value set to 100
tmp = np.compress(image.ravel()<100,image.ravel())
if opt.sbmax==None:
opt.sbmax = tmp.max()
if opt.sbmin==None:
opt.sbmin = tmp.min()
# plot image
im = axii.imshow(image,aspect='equal',extent=(imgext_xmin,imgext_xmax,imgext_ymin,imgext_ymax),interpolation=None,vmin=opt.sbmin,vmax=opt.sbmax,cmap=cmap)
if opt.sbcontours is not None:
axii.contour(image,levels=opt.sbcontours,colors='k',linewidths=1,linestyles='solid',aspect='equal',extent=(imgext_xmin,imgext_xmax,imgext_ymin,imgext_ymax),origin='upper')
# labels
axii.set_xlabel(xlabel)
axii.set_ylabel(ylabel)
# limits
axii.set_xlim(xmin,xmax)
axii.set_ylim(ymin,ymax)
# add axes
if not opt.add_axes:
axii.get_xaxis().set_visible(False)
axii.get_yaxis().set_visible(False)
axii.get_xaxis().set_ticks([])
axii.get_yaxis().set_ticks([])
# add colorbar
if opt.colorbar:
plt.colorbar(im,label="surface brightness [mag/arcsec]",ax=axii,location='right')
#txt = r"$\rm{%s}\,\,:\rm{%s}\,\,:\,\,\rm{D=%4.1f}\,\rm{Mpc}$"%(params["name"],params["filter_name"],params["object_distance"].to(u.Mpc).value)
#txt = r"%s : %s : D=%4.1f Mpc"%(params["name"],params["filter_name"],params["object_distance"].to(u.Mpc).value)
#fig.suptitle(txt,fontsize=20)
#txt = r"$\rm{%s}\,\,:\rm{%s}\,\,:\,\,\rm{D=%4.1f}\,\rm{Mpc}$"%(header["NAME"],header["FILTER"],header["OBJ_DIST"])
los = header["OBJ_LOS"]
los = los[1:-1].split(',')
los = np.array(los,float)
#los = "({0:3.1f},{1:3.1f},{2:3.1f})".format(los[0],los[1],los[2])
los = "los=({0:1.0f},{1:1.0f},{2:1.0f})".format(los[0],los[1],los[2])
txt = r"%s : %s : D=%4.1f Mpc %s"%(header["NAME"],header["FILTER"],header["OBJ_DIST"],los)
plt.title(txt,fontsize=20,loc="left")
# add time
if "OBJ_TNOW" in header:
ax = plt.gca()
tnow = "{0:5.2f} [Gyr]".format(float(header["OBJ_TNOW"]))
ax.text(0,ymin*0.75,tnow,ha='center',va='center',fontsize=20)
if opt.outputfilename:
plt.savefig(opt.outputfilename)
else:
plt.show()

Event Timeline

c4science · Help