Homec4science
Diffusion ARRAKIHS (master) Browse

arrakihs/Workspace/PAndAS/mockimgs_sb_display_fits
96a87b475cb8master

Last Change
Raw File

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("--colormap",
action="store",
default="Greys",
help='matplotlib colormap name (e.g. tab20c)')
parser.add_argument("--colormap_reverse",
action="store_true",
default=False,
help='reverse colormap')
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('--sbmaxcut',
action="store",
dest="sbmaxcut",
metavar='FLOAT',
type=float,
default=None,
help='cut pixels with sb > sbmaxcut')
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))
#cmap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors2)
cmap = mpl.colormaps[opt.colormap]
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.sbmaxcut!=None:
image = np.where(image>opt.sbmaxcut,100,image)
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}$"%(header["NAME"],header["FILTER"],header["OBJ_DIST"])
txt = r"%s : %s : D=%4.1f Mpc"%(header["NAME"],header["FILTER"],header["OBJ_DIST"])
fig.suptitle(txt,fontsize=20)
if opt.outputfilename:
plt.savefig(opt.outputfilename)
else:
plt.show()
c4science · Help