Page MenuHomec4science
Files F102801227

mockimgs_sb_display_fits
No OneTemporary
Actions

File Metadata

Created
Mon, Feb 24, 07:54

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')
####################################################################
# 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"]
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=mymap)
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)
if opt.outputfilename:
plt.savefig(opt.outputfilename)
else:
plt.show()

Event Timeline

c4science · Help