Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F83455525
mockimgs_sb_display_fits
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 17, 05:51
Size
13 KB
Mime Type
text/x-python
Expires
Thu, Sep 19, 05:51 (2 d)
Engine
blob
Format
Raw Data
Handle
20841072
Attached To
rARRAKIHS ARRAKIHS
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
Log In to Comment