# coding=utf8
# Copyright (c) MPE/IR-Submm Group. See LICENSE.rst for license information.
#
# Functions for plotting DYSMALPY kinematic model fit results
from __future__ import (absolute_import, division, print_function,
unicode_literals)
# Standard library
import logging
import copy
from collections import OrderedDict
import os
import datetime
# Third party imports
import numpy as np
import astropy.units as u
import matplotlib as mpl
# Check if there is a display for plotting, or if there is an SSH/TMUX session.
# If no display, or if SSH/TMUX, use the matplotlib "agg" backend for plotting.
havedisplay = "DISPLAY" in os.environ
if havedisplay:
exitval = os.system('python -c "import matplotlib.pyplot as plt; plt.figure()"')
skipconds = (("SSH_CLIENT" in os.environ) | ("TMUX" in os.environ) | ("SSH_CONNECTION" in os.environ) | (os.environ["TERM"].lower().strip()=='screen') | (exitval != 0))
havedisplay = not skipconds
if not havedisplay:
mpl.use('agg')
import astropy.modeling as apy_mod
import astropy.io.fits as fits
from astropy.wcs import WCS
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.ticker import MultipleLocator
import matplotlib.colors as mplcolors
from mpl_toolkits.axes_grid1 import ImageGrid, AxesGrid
from matplotlib import colorbar
import corner
from spectral_cube import SpectralCube, BooleanArrayMask
# Package imports
from .utils import calc_pix_position, apply_smoothing_3D, gaus_fit_sp_opt_leastsq, gaus_fit_apy_mod_fitter
from .utils_io import create_vel_profile_files_obs1d, create_vel_profile_files_intrinsic, \
read_bestfit_1d_obs_file, read_model_intrinsic_profile
from .observation import Observation
from .aperture_classes import CircApertures
from .data_classes import Data1D, Data2D
from .extern.altered_colormaps import new_diverging_cmap
#from .config import Config_create_model_data, Config_simulate_cube
try:
from dysmalpy.utils_least_chi_squares_1d_fitter import LeastChiSquares1D
_loaded_LeastChiSquares1D = True
except:
_loaded_LeastChiSquares1D = False
# New colormap:
new_diverging_cmap('RdBu_r', diverge = 0.5,
gamma_lower=1.5, gamma_upper=1.5,
name_new='RdBu_r_stretch')
# Base colormaps:
try:
cmap_viridis = mpl.colormap['viridis']
cmap_spectral_r = mpl.colormap['Spectral_r']
cmap_greys = mpl.colormaps['Greys']
cmap_plasma = mpl.colormaps['plasma']
cmap_rdbu_r = mpl.colormaps["RdBu_r_stretch"]
cmap_seismic = mpl.colormaps['seismic']
except:
import matplotlib.cm as cm
cmap_viridis = cm.viridis
cmap_spectral_r = mpl.colormaps["Spectral_r"]
cmap_greys = cm.Greys
cmap_plasma = cm.plasma
cmap_rdbu_r = mpl.colormaps["RdBu_r_stretch"]
cmap_seismic = cm.seismic
# Default settings for contours on 2D maps:
_kwargs_contour_defaults = { 'colors_cont': 'black',
'alpha_cont': 1.,
'ls_cont': 'solid',
'lw_cont': 0.75,
'delta_cont_v': 25.,
'delta_cont_disp': 25.,
'delta_cont_flux': 20., #5.,
####
'lw_cont_minor': 0.2,
'alpha_cont_minor': 1.,
'colors_cont_minor': 'black',
'ls_cont_minor': '-',
'delta_cont_v_minor': None,
'delta_cont_disp_minor': None,
'delta_cont_flux_minor': None,
####
'lw_cont_minor_resid': 0.15,
'colors_cont_minor_resid': 'grey',
'alpha_cont_minor_resid' : 1.,
'ls_cont_minor_resid' : ':',
'delta_cont_v_minor_resid': None,
'delta_cont_disp_minor_resid': None,
'delta_cont_flux_minor_resid': None
}
__all__ = ['plot_trace', 'plot_corner', 'plot_bestfit']
# LOGGER SETTINGS
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger('DysmalPy')
logger.setLevel(logging.INFO)
[docs]
def plot_bestfit(mcmcResults, gal,
show_1d_apers=False,
fileout=None,
vcrop=False,
vcrop_value=800.,
remove_shift=False,
overwrite=False,
fill_mask=False,
**plot_kwargs):
"""
Plot data, bestfit model, and residuals from the MCMC fitting.
"""
plot_data_model_comparison(gal, theta = mcmcResults.bestfit_parameters,
fileout=fileout,
vcrop=vcrop, vcrop_value=vcrop_value, show_1d_apers=show_1d_apers,
remove_shift=remove_shift, fill_mask=fill_mask,
overwrite=overwrite, **plot_kwargs)
return None
[docs]
def plot_trace(bayesianResults, fileout=None, overwrite=False):
"""
Plot trace of Bayesian samples
"""
if bayesianResults.fit_method.lower() == 'mcmc':
plot_trace_mcmc(bayesianResults, fileout=fileout, overwrite=overwrite)
elif bayesianResults.fit_method.lower() == 'nested':
plot_trace_nested(bayesianResults, fileout=fileout, overwrite=overwrite)
else:
raise ValueError("plot_trace() not supported for fit method: {}".format(bayesianResults.fit_method))
def plot_run(bayesianResults, fileout=None, overwrite=False):
if bayesianResults.fit_method.lower() == 'nested':
plot_run_nested(bayesianResults, fileout=fileout, overwrite=overwrite)
else:
raise ValueError("plot_run() not supported for fit method: {}".format(bayesianResults.fit_method))
[docs]
def plot_corner(bayesianResults, gal=None, fileout=None,
step_slice=None, blob_name=None, overwrite=False):
"""
Plot corner plot of Bayesian result posterior distributions.
Optional:
step slice: 2 element tuple/array with beginning and end step number to use
"""
if bayesianResults.fit_method.lower() == 'mcmc':
plot_corner_mcmc(bayesianResults, gal=gal, fileout=fileout,
step_slice=step_slice, blob_name=blob_name,
overwrite=overwrite)
elif bayesianResults.fit_method.lower() == 'nested':
plot_corner_nested(bayesianResults, fileout=fileout, overwrite=overwrite)
else:
raise ValueError("plot_corner() not supported for fit method: {}".format(bayesianResults.fit_method))
def plot_trace_mcmc(mcmcResults, fileout=None, overwrite=False):
"""
Plot trace of MCMC walkers
"""
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
names = make_clean_bayesian_plot_names(mcmcResults)
######################################
# Setup plot:
f = plt.figure()
scale = 1.75
nRows = len(names)
nWalkers = mcmcResults.sampler_results['nWalkers']
f.set_size_inches(4.*scale, nRows*scale)
gs = gridspec.GridSpec(nRows, 1, hspace=0.2)
axes = []
alpha = max(0.01, min(10./nWalkers, 1.))
# Define random color inds for tracking some walkers:
nTraceWalkers = 5
cmap = cmap_viridis
alphaTrace = 0.8
lwTrace = 1.5
trace_inds = np.random.randint(0,nWalkers, size=nTraceWalkers)
trace_colors = []
for i in range(nTraceWalkers):
trace_colors.append(cmap(1./float(nTraceWalkers)*i))
norm_inds = np.setdiff1d(range(nWalkers), trace_inds)
for k in range(nRows):
axes.append(plt.subplot(gs[k,0]))
axes[k].plot(mcmcResults.sampler_results['chain'][norm_inds,:,k].T, '-',
color='black', alpha=alpha, rasterized=True)
for j in range(nTraceWalkers):
axes[k].plot(mcmcResults.sampler_results['chain'][trace_inds[j],:,k].T, '-',
color=trace_colors[j], lw=lwTrace, alpha=alphaTrace)
axes[k].set_ylabel(names[k])
if k == nRows-1:
axes[k].set_xlabel('Step number')
else:
axes[k].set_xticks([])
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.show()
return None
def plot_trace_nested(bayesianResults, fileout=None, overwrite=False):
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
from dynesty import plotting as dyplot
names = make_clean_bayesian_plot_names(bayesianResults)
scale = 1.5
ndim = bayesianResults.nparams_free
figsize = (4*scale, ndim*scale)
# Plot dynesty trace plot:
f, axes = dyplot.traceplot(bayesianResults.sampler_results,
truths=bayesianResults.bestfit_parameters,
labels=names, show_titles=True,
trace_cmap='viridis',
title_kwargs={'fontsize': 15, 'y': 1.05},
quantiles=None,
fig=plt.subplots(ndim, 2, figsize=figsize))
f.tight_layout()
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.show()
return None
def plot_run_nested(bayesianResults, fileout=None, overwrite=False):
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
from dynesty import plotting as dyplot
f, axes = dyplot.runplot(bayesianResults.sampler_results)
scale = 1.5
nRows = bayesianResults.nparams_free
f.set_size_inches(4*scale, nRows*scale)
f.tight_layout()
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.show()
return None
def plot_corner_mcmc(mcmcResults, gal=None, fileout=None, step_slice=None,
blob_name=None, overwrite=False):
"""
Plot corner plot of MCMC result posterior distributions.
Optional:
step slice: 2 element tuple/array with beginning and end step number to use
"""
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
names = make_clean_bayesian_plot_names(mcmcResults)
if step_slice is None:
sampler_chain = mcmcResults.sampler_results['flatchain']
else:
sampler_chain = mcmcResults.sampler_results['chain'][:,step_slice[0]:step_slice[1],:].reshape((-1, mcmcResults.sampler_results['nParam']))
truths = mcmcResults.bestfit_parameters
truths_l68 = mcmcResults.bestfit_parameters_l68_err
truths_u68 = mcmcResults.bestfit_parameters_u68_err
try:
truths_l68_percentile = mcmcResults.bestfit_parameters_l68_err_percentile
truths_u68_percentile = mcmcResults.bestfit_parameters_u68_err_percentile
except:
truths_l68_percentile = None
truths_u68_percentile = None
if gal is not None:
# Get prior locations:
priors = []
pfree_dict = gal.model.get_free_parameter_keys()
comps_names = pfree_dict.keys()
for compn in comps_names:
comp = gal.model.components.__getitem__(compn)
params_names = pfree_dict[compn].keys()
for paramn in params_names:
if pfree_dict[compn][paramn] >= 0:
# Free parameter:
# check if uniform or prior
if 'center' in comp.__getattribute__(paramn).prior.__dict__.keys():
priors.append(comp.__getattribute__(paramn).prior.center)
else:
priors.append(None)
else:
priors = None
###############
if blob_name is not None:
names_nice = names[:]
if isinstance(blob_name, str):
blob_names = [blob_name]
else:
blob_names = blob_name[:]
for blobn in blob_names:
blob_true = mcmcResults.__dict__['bestfit_{}'.format(blobn)]
blob_l68_err = mcmcResults.__dict__['bestfit_{}_l68_err'.format(blobn)]
blob_u68_err = mcmcResults.__dict__['bestfit_{}_u68_err'.format(blobn)]
if blobn.lower() == 'fdm':
names.append('Blob: fDM(RE)')
names_nice.append(r'$f_{\mathrm{DM}}(R_E)$')
elif blobn.lower() == 'alpha':
names.append('Blob: alpha')
names_nice.append(r'$\alpha$')
elif blobn.lower() == 'mvirial':
names.append('Blob: Mvirial')
names_nice.append(r'$\log_{10}(M_{\rm vir})$')
elif blobn.lower() == 'rb':
names.append('Blob: rB')
names_nice.append(r'$R_B$')
else:
names.append(blobn)
names_nice.append(blobn)
if isinstance(blob_name, str):
if step_slice is None:
blobs = mcmcResults.sampler_results['flatblobs']
else:
blobs = mcmcResults.sampler_results['blobs'][step_slice[0]:step_slice[1],:,:].reshape((-1, 1))
else:
indv = blob_names.index(blobn)
if step_slice is None:
blobs = mcmcResults.sampler_results['flatblobs'][:,indv]
else:
blobs = mcmcResults.sampler_results['blobs'][step_slice[0]:step_slice[1],:,:].reshape((-1, mcmcResults.sampler_results['blobs'].shape[2]))[:,indv]
sampler_chain = np.concatenate( (sampler_chain, np.array([blobs]).T ), axis=1)
truths = np.append(truths, blob_true)
truths_l68 = np.append(truths_l68, blob_l68_err)
truths_u68 = np.append(truths_u68, blob_u68_err )
if priors is not None:
priors.append(None)
if truths_l68_percentile is not None:
truths_l68_percentile = np.append(truths_l68_percentile, mcmcResults.__dict__['bestfit_{}_l68_err_percentile'.format(blobn)])
truths_u68_percentile = np.append(truths_u68_percentile, mcmcResults.__dict__['bestfit_{}_u68_err_percentile'.format(blobn)])
title_kwargs = {'horizontalalignment': 'left', 'x': 0.}
fig = corner.corner(sampler_chain,
labels=names,
title_quantiles=[0.16,0.5,0.84], # Required by corner, will be overwritten later with self-calculated truth/l/u68
quantiles= [.02275, 0.15865, 0.84135, .97725], # Plot raw upper, lower 1 and 2 sigma quantiles
truths=truths,
plot_datapoints=False,
show_titles=True,
bins=40,
plot_contours=True,
verbose=False,
title_kwargs=title_kwargs)
axes = fig.axes
nFreeParam = len(truths)
for i in range(nFreeParam):
ax = axes[i*nFreeParam + i]
# Format the quantile display.
best = truths[i]
q_m = truths_l68[i]
q_p = truths_u68[i]
title_fmt=".2f"
fmt = "{{0:{0}}}".format(title_fmt).format
title = r"${{{0}}}_{{-{1}}}^{{+{2}}}$"
title = title.format(fmt(best), fmt(q_m), fmt(q_p))
# Add in the column name if it's given.
if names is not None:
title = "{0} = {1}".format(names[i], title)
ax.set_title(title, **title_kwargs)
xlim = ax.get_xlim()
ylim = ax.get_ylim()
if truths_l68_percentile is not None:
if (truths_l68_percentile[i] != q_m) | (truths_u68_percentile[i] != q_p):
ax.axvline(best-q_m, ls='--', color='#9467bd') # purple
ax.axvline(best+q_p, ls='--', color='#9467bd') # purple
if priors is not None:
if priors[i] is not None:
ax.axvline(priors[i], ls=':', color='#ff7f0e') # orange
ax.set_xlim(xlim)
ax.set_ylim(ylim)
if priors is not None:
for i in range(nFreeParam):
for j in range(nFreeParam):
# need the off-diagonals:
if j >= i:
pass
else:
ax = axes[i*nFreeParam + j]
xlim = ax.get_xlim()
ylim = ax.get_ylim()
if priors[i] is not None:
ax.axhline(priors[i], ls=':', color='#ff7f0e') # orange
if priors[j] is not None:
ax.axvline(priors[j], ls=':', color='#ff7f0e') # orange
if (priors[i] is not None) & (priors[j] is not None):
ax.scatter([priors[j]], [priors[i]], marker='s', edgecolor='#ff7f0e', facecolor='None') # orange
#
ax.set_xlim(xlim)
ax.set_ylim(ylim)
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight')
plt.close(fig)
else:
plt.show()
return None
def plot_corner_nested(bayesianResults, fileout=None, overwrite=False):
from dynesty import plotting as dyplot
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
MAP_vals = bayesianResults.bestfit_parameters
ndim = bayesianResults.nparams_free
names = make_clean_bayesian_plot_names(bayesianResults, short=True)
# Plot dynesty corner plot:
# initialize figure
scale = 1.5
f, axes = plt.subplots(ndim, ndim, figsize=(scale*ndim,scale*ndim))
axes = axes.reshape((ndim, ndim))
fg, ax = dyplot.cornerplot(bayesianResults.sampler_results, color='blue',
truths=MAP_vals,
#span=[(0,1), (0,1), (0,1), (0,1)],
labels=names,
show_titles=True, title_kwargs={'y': 1.05},
quantiles=None, fig=(f, axes[:, :]))
title_fmt = '.2f'
for i in range(ndim):
axi = ax[i,i]
qm = MAP_vals[i]
ql = bayesianResults.bestfit_parameters_l68[i]
qh = bayesianResults.bestfit_parameters_u68[i]
q_minus, q_plus = qm - ql, qh - qm
fmt = "{{0:{0}}}".format(title_fmt).format
title = r"${{{0}}}_{{-{1}}}^{{+{2}}}$"
title = title.format(fmt(qm), fmt(q_minus), fmt(q_plus))
title = "{0} = {1}".format(names[i], title)
axi.set_title(title)
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
return None
# -------------------------------------------------------------
def plot_data_model_comparison(gal,theta=None,
fileout=None,
vcrop=False,
show_1d_apers=False,
vcrop_value=800.,
remove_shift=False,
overwrite=False,
fill_mask=False,
show_contours=False,
show_ruler=True,
ruler_loc='lowerleft',
**plot_kwargs):
dummy_gal = copy.deepcopy(gal)
if remove_shift:
dummy_gal.model.geometries[obs.name].xshift = 0
dummy_gal.model.geometries[obs.name].yshift = 0
# Update model parameters and create new set of model data if
# theta array provided
if theta is not None:
dummy_gal.model.update_parameters(theta)
dummy_gal.create_model_data()
# Plot data and model for each observation
for obs_name in dummy_gal.observations:
obs = dummy_gal.observations[obs_name]
plot_single_obs_data_model_comparison(obs, dummy_gal.model,
fileout=fileout,
dscale=dummy_gal.dscale,
vcrop=vcrop,
show_1d_apers=show_1d_apers,
vcrop_value=vcrop_value,
overwrite=overwrite,
fill_mask=fill_mask,
show_contours=show_contours,
show_ruler=show_ruler,
ruler_loc=ruler_loc,
**plot_kwargs)
def plot_single_obs_data_model_comparison(obs, model, theta = None,
fileout=None,
dscale=None,
vcrop=False,
show_1d_apers=False,
vcrop_value=800.,
remove_shift=False,
overwrite=False,
fill_mask=False,
show_contours=False,
show_ruler=True,
ruler_loc='lowerleft',
**plot_kwargs):
"""
Plot data, model, and residuals between the data and this model.
"""
dummy_obs = copy.deepcopy(obs)
dummy_model = copy.deepcopy(model)
if remove_shift:
dummy_model.geometries[obs.name].xshift = 0
dummy_model.geometries[obs.name].yshift = 0
if fileout is not None:
fileout_in = fileout
f_r = fileout.split('.')
f_base = ".".join(f_r[:-1])
fileout = "{}_{}.{}".format(f_base, obs.name, f_r[-1])
if (not overwrite) & os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
if dummy_obs.instrument.ndim == 1:
plot_data_model_comparison_1D(dummy_obs, fileout=fileout,
overwrite=overwrite,
**plot_kwargs)
elif dummy_obs.instrument.ndim == 2:
plot_data_model_comparison_2D(dummy_obs, dummy_model,
fileout=fileout,
show_contours=show_contours,
show_ruler=show_ruler,
ruler_loc=ruler_loc,
overwrite=overwrite,
**plot_kwargs)
elif dummy_obs.instrument.ndim == 3:
plot_data_model_comparison_3D(dummy_obs, dummy_model,
show_1d_apers=show_1d_apers,
fileout=fileout,
dscale=dscale,
vcrop=vcrop,
vcrop_value=vcrop_value,
overwrite=overwrite,
fill_mask=fill_mask,
show_contours=show_contours,
show_ruler=show_ruler,
ruler_loc=ruler_loc,
**plot_kwargs)
elif dummy_obs.instrument.ndim == 0:
plot_data_model_comparison_0D(dummy_obs, fileout=fileout,
overwrite=overwrite,
**plot_kwargs)
else:
logger.warning("nDim="+str(dummy_obs.instrument.ndim)+" not supported!")
raise ValueError("nDim="+str(dummy_obs.instrument.ndim)+" not supported!")
return None
def plot_data_model_comparison_0D(obs, fileout=None,
overwrite=False):
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
data = obs.data
model_data = obs.model_data
######################################
# Setup plot:
scale = 3.5
f = plt.figure(figsize=(2.2 * scale, scale))
ax = f.add_subplot(111)
# Plot the observed spectrum with error shading
ax.plot(data.x, data.data, color='black', lw=1.5)
ax.fill_between(data.x, data.data - data.error, data.data + data.error, color='black', alpha=0.2)
# Plot the model spectrum
ax.plot(model_data.x, model_data.data, color='red', lw=1.5)
# Plot the residuals
ax.plot(model_data.x, data.data - model_data.data, color='blue', lw=1.0)
ax.set_ylabel('Normalized Flux')
ax.set_xlabel(obs.instrument.spec_type.capitalize() + ' [' + obs.instrument.spec_step.unit.name + ']')
f.suptitle(obs.name)
# Save to file:
if fileout is not None:
f.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close(f)
else:
plt.show()
def plot_data_model_comparison_1D(obs, fileout=None, overwrite=False):
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
# Default: fit in 1D, compare to 1D data:
data = obs.data
model_data = obs.model_data
# Correct model for instrument dispersion if the data is instrument corrected:
if 'inst_corr' in data.data.keys():
inst_corr = data.data['inst_corr']
if inst_corr:
model_data.data['dispersion'] = \
np.sqrt( model_data.data['dispersion']**2 - \
obs.instrument.lsf.dispersion.to(u.km/u.s).value**2 )
######################################
# Setup plot:
f = plt.figure()
scale = 3.5
ncols = 0
for cond in [obs.fit_options.fit_flux, obs.fit_options.fit_velocity, obs.fit_options.fit_dispersion]:
if cond:
ncols += 1
nrows = 2
f.set_size_inches(1.1*ncols*scale, nrows*scale)
gs = gridspec.GridSpec(nrows, ncols, wspace=0.35, hspace=0.2)
keyxtitle = r'$r$ [arcsec]'
keyyarr, keyytitlearr, keyyresidtitlearr = ([] for _ in range(3))
if obs.fit_options.fit_flux:
keyyarr.append('flux')
keyytitlearr.append(r'Flux [arb]')
keyyresidtitlearr.append(r'$\mathrm{Flux_{data} - Flux_{model}}$ [arb]')
if obs.fit_options.fit_velocity:
keyyarr.append('velocity')
keyytitlearr.append(r'$V$ [km/s]')
keyyresidtitlearr.append(r'$V_{\mathrm{data}} - V_{\mathrm{model}}$ [km/s]')
if obs.fit_options.fit_dispersion:
keyyarr.append('dispersion')
keyytitlearr.append(r'$\sigma$ [km/s]')
keyyresidtitlearr.append(r'$\sigma_{\mathrm{data}} - \sigma_{\mathrm{model}}$ [km/s]')
errbar_lw = 0.5
errbar_cap = 1.5
axes = []
k = -1
for j in range(ncols):
# Comparison:
axes.append(plt.subplot(gs[0,j]))
k += 1
if keyyarr[j] == 'velocity':
if hasattr(obs.data, 'mask_velocity'):
if obs.data.mask_velocity is not None:
msk = obs.data.mask_velocity
else:
msk = obs.data.mask
else:
msk = obs.data.mask
elif keyyarr[j] == 'dispersion':
if hasattr(obs.data, 'mask_vel_disp'):
if obs.data.mask_vel_disp is not None:
msk = obs.data.mask_vel_disp
else:
msk = obs.data.mask
else:
msk = obs.data.mask
elif keyyarr[j] == 'flux':
msk = obs.data.mask
else:
msk = np.array(np.ones(len(data.rarr)), dtype=bool)
# Masked points
if data.error[keyyarr[j]] is not None:
axes[k].errorbar( data.rarr[~msk], data.data[keyyarr[j]][~msk],
xerr=None, yerr=data.error[keyyarr[j]][~msk],
marker=None, ls='None', ecolor='darkgrey', zorder=-1.,
lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None, alpha=0.5 )
if np.any(~msk):
lbl_mask = 'Masked data'
else:
lbl_mask = None
axes[k].scatter( data.rarr[~msk], data.data[keyyarr[j]][~msk],
c='darkgrey', marker='o', s=25, lw=1, label=lbl_mask, alpha=0.5)
# Unmasked points
if data.error[keyyarr[j]] is not None:
axes[k].errorbar( data.rarr[msk], data.data[keyyarr[j]][msk],
xerr=None, yerr=data.error[keyyarr[j]][msk],
marker=None, ls='None', ecolor='k', zorder=-1.,
lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None )
# Weights: effective errorbars show in blue, for reference
if hasattr(data, 'weight'):
if obs.data.weight is not None:
wgt = obs.data.weight
axes[k].errorbar( data.rarr[msk], data.data[keyyarr[j]][msk],
xerr=None, yerr=data.error[keyyarr[j]][msk]/np.sqrt(wgt[msk]),
marker=None, ls='None', ecolor='blue', zorder=-1.,
lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None )
axes[k].scatter( data.rarr[msk], data.data[keyyarr[j]][msk],
c='black', marker='o', s=25, lw=1, label='Data')
# Masked points
axes[k].scatter( model_data.rarr[~msk], model_data.data[keyyarr[j]][~msk],
c='lightsalmon', marker='s', s=25, lw=1, label=None, alpha=0.3)
# Unmasked points
axes[k].scatter( model_data.rarr[msk], model_data.data[keyyarr[j]][msk],
c='red', marker='s', s=25, lw=1, label='Model')
axes[k].set_xlabel(keyxtitle)
axes[k].set_ylabel(keyytitlearr[j])
axes[k].axhline(y=0, ls='--', color='k', zorder=-10.)
if k == 0:
handles, lbls = axes[k].get_legend_handles_labels()
if len(lbls) > 2:
ind_reord = [1,0,2]
labels_arr = []
handles_arr = []
for ir in ind_reord:
labels_arr.append(lbls[ir])
handles_arr.append(handles[ir])
else:
labels_arr = lbls
handles_arr = handles
frameon = True
borderpad = 0.25
markerscale = 0.8 #1.
labelspacing= 0.25
handletextpad = 0.2
handlelength = 1.
fontsize_leg= 7.5
legend = axes[k].legend(handles_arr, labels_arr,
labelspacing=labelspacing, borderpad=borderpad,
markerscale=markerscale,
handletextpad=handletextpad,
handlelength=handlelength,
loc='lower right',
frameon=frameon, numpoints=1,
scatterpoints=1,
fontsize=fontsize_leg)
# Residuals:
axes.append(plt.subplot(gs[1,j]))
k += 1
# Masked points
if data.error[keyyarr[j]] is not None:
axes[k].errorbar( data.rarr[~msk], data.data[keyyarr[j]][~msk]-model_data.data[keyyarr[j]][~msk],
xerr=None, yerr = data.error[keyyarr[j]][~msk],
marker=None, ls='None', ecolor='darkgrey', zorder=-1.,
lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None, alpha=0.5 )
axes[k].scatter( data.rarr[~msk], data.data[keyyarr[j]][~msk]-model_data.data[keyyarr[j]][~msk],
c='lightsalmon', marker='s', s=25, lw=1, label=None, alpha=0.3)
# Unmasked points:
if data.error[keyyarr[j]] is not None:
axes[k].errorbar( data.rarr[msk], data.data[keyyarr[j]][msk]-model_data.data[keyyarr[j]][msk],
xerr=None, yerr = data.error[keyyarr[j]][msk],
marker=None, ls='None', ecolor='k', zorder=-1.,
lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None )
# Weights: effective errorbars show in blue, for reference
if hasattr(data, 'weight'):
if obs.data.weight is not None:
wgt = obs.data.weight
axes[k].errorbar( data.rarr[msk], data.data[keyyarr[j]][msk]-model_data.data[keyyarr[j]][msk],
xerr=None, yerr = data.error[keyyarr[j]][msk]/np.sqrt(wgt[msk]),
marker=None, ls='None', ecolor='blue', zorder=-1.,
lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None )
#
axes[k].scatter( data.rarr[msk], data.data[keyyarr[j]][msk]-model_data.data[keyyarr[j]][msk],
c='red', marker='s', s=25, lw=1, label=None)
axes[k].axhline(y=0, ls='--', color='k', zorder=-10.)
axes[k].set_xlabel(keyxtitle)
axes[k].set_ylabel(keyyresidtitlearr[j])
# Figure title is the observation name
f.suptitle(obs.name)
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.show()
def plot_data_model_comparison_2D(obs, model,
fileout=None,
symmetric_residuals=True,
max_residual=100.,
max_residual_flux=None,
overwrite=False,
show_contours=False,
show_ruler=True,
ruler_loc='lowerleft',
**plot_kwargs):
if show_contours:
# Set contour defaults, if not specifed:
for key in _kwargs_contour_defaults.keys():
if key not in plot_kwargs.keys():
plot_kwargs[key] = _kwargs_contour_defaults[key]
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
try:
if 'inst_corr' in obs.data.data.keys():
inst_corr = obs.data.data['inst_corr']
except:
inst_corr = False
######################################
# Setup plot:
f = plt.figure(figsize=(9.5, 6))
scale = 3.5
# Plot settings:
nrows_ncols=(1, 3)
direction="row"
axes_pad=0.5
label_mode="1"
share_all=True
cbar_location="right"
cbar_mode="each"
cbar_size="5%"
cbar_pad="1%"
nrows = 0
for cond in [obs.fit_options.fit_flux, obs.fit_options.fit_velocity, obs.fit_options.fit_dispersion]:
if cond:
nrows += 1
cntr = 0
if obs.fit_options.fit_flux:
cntr += 1
grid_flux = ImageGrid(f, 100*nrows+10+cntr,
nrows_ncols=nrows_ncols,
direction=direction,
axes_pad=axes_pad,
label_mode=label_mode,
share_all=share_all,
cbar_location=cbar_location,
cbar_mode=cbar_mode,
cbar_size=cbar_size,
cbar_pad=cbar_pad)
if obs.fit_options.fit_velocity:
cntr += 1
grid_vel = ImageGrid(f, 100*nrows+10+cntr,
nrows_ncols=nrows_ncols,
direction=direction,
axes_pad=axes_pad,
label_mode=label_mode,
share_all=share_all,
cbar_location=cbar_location,
cbar_mode=cbar_mode,
cbar_size=cbar_size,
cbar_pad=cbar_pad)
if obs.fit_options.fit_dispersion:
cntr += 1
grid_disp = ImageGrid(f, 100*nrows+10+cntr,
nrows_ncols=nrows_ncols,
direction=direction,
axes_pad=axes_pad,
label_mode=label_mode,
share_all=share_all,
cbar_location=cbar_location,
cbar_mode=cbar_mode,
cbar_size=cbar_size,
cbar_pad=cbar_pad)
keyxarr = ['data', 'model', 'residual']
keyxtitlearr = ['Data', 'Model', 'Residual']
keyyarr, keyytitlearr, grid_arr, annstr_arr = ([] for _ in range(4))
if obs.fit_options.fit_flux:
keyyarr.append('flux')
keyytitlearr.append(r'Flux')
grid_arr.append(grid_flux)
annstr_arr.append('f')
if obs.data is not None:
flux_vmin = obs.data.data['flux'][obs.data.mask].min()
flux_vmax = obs.data.data['flux'][obs.data.mask].max()
if flux_vmin == flux_vmax:
flux_vmin = obs.model_data.data['flux'].min()
flux_vmax = obs.model_data.data['flux'].max()
else:
flux_vmin = obs.model_data.data['flux'].min()
flux_vmax = obs.model_data.data['flux'].max()
if max_residual_flux is None:
max_residual_flux = np.max(np.abs(obs.data.data['flux'][obs.data.mask]))
if obs.fit_options.fit_velocity:
keyyarr.append('velocity')
keyytitlearr.append(r'$V$')
grid_arr.append(grid_vel)
annstr_arr.append('V')
vel_vmin = obs.data.data['velocity'][obs.data.mask].min()
vel_vmax = obs.data.data['velocity'][obs.data.mask].max()
try:
vel_shift = model.geometries[obs.name].vel_shift.value
except:
vel_shift = 0
vel_vmin -= vel_shift
vel_vmax -= vel_shift
if obs.fit_options.fit_dispersion:
keyyarr.append('dispersion')
keyytitlearr.append(r'$\sigma$')
grid_arr.append(grid_disp)
annstr_arr.append('\sigma')
if obs.data is not None:
disp_vmin = obs.data.data['dispersion'][obs.data.mask].min()
disp_vmax = obs.data.data['dispersion'][obs.data.mask].max()
else:
disp_vmin = obs.model_data.data['dispersion'].min()
disp_vmax = obs.model_data.data['dispersion'].max()
int_mode = "nearest"
origin = 'lower'
cmap = copy.copy(cmap_spectral_r)
# color_bad = 'black'
# color_annotate = 'white'
color_bad = 'white'
color_annotate = 'black'
cmap.set_bad(color=color_bad)
cmap_resid = copy.copy(cmap_rdbu_r)
cmap_resid.set_bad(color=color_bad)
for j in range(len(keyyarr)):
grid = grid_arr[j]
for ax, k, xt in zip(grid, keyxarr, keyxtitlearr):
if k == 'data':
im = obs.data.data[keyyarr[j]].copy()
if keyyarr[j] == 'velocity':
im -= vel_shift
vmin = vel_vmin
vmax = vel_vmax
elif keyyarr[j] == 'dispersion':
vmin = disp_vmin
vmax = disp_vmax
elif keyyarr[j] == 'flux':
vmin = flux_vmin
vmax = flux_vmax
cmaptmp = cmap
elif k == 'model':
im = obs.model_data.data[keyyarr[j]].copy()
if keyyarr[j] == 'velocity':
im -= vel_shift
vmin = vel_vmin
vmax = vel_vmax
elif keyyarr[j] == 'dispersion':
if inst_corr:
im = np.sqrt(im ** 2 - obs.instrument.lsf.dispersion.to(
u.km / u.s).value ** 2)
vmin = disp_vmin
vmax = disp_vmax
elif keyyarr[j] == 'flux':
vmin = flux_vmin
vmax = flux_vmax
cmaptmp = cmap
elif k == 'residual':
im_data = obs.data.data[keyyarr[j]].copy()
im_model = obs.model_data.data[keyyarr[j]].copy()
if keyyarr[j] == 'dispersion':
if inst_corr:
im_model = np.sqrt(im_model ** 2 -
obs.instrument.lsf.dispersion.to( u.km / u.s).value ** 2)
im = im_data - im_model
if symmetric_residuals:
if keyyarr[j] == 'flux':
vmin = -max_residual_flux
vmax = max_residual_flux
else:
vmin = -max_residual
vmax = max_residual
cmaptmp = cmap_resid
else:
raise ValueError("key not supported.")
# Mask image:
im[~obs.data.mask] = np.nan
imax = ax.imshow(im, cmap=cmaptmp, interpolation=int_mode,
vmin=vmin, vmax=vmax, origin=origin)
if len(model.geometries) > 0:
ax = plot_major_minor_axes_2D(ax, obs, model, im, obs.data.mask)
if show_ruler:
pixscale = obs.instrument.pixscale.value
ax = plot_ruler_arcsec_2D(ax, pixscale, len_arcsec=1.,
ruler_loc=ruler_loc, color=color_annotate)
if show_contours:
ax = plot_contours_2D_multitype(im, ax=ax, mapname=keyyarr[j], plottype=k,
vmin=vmin, vmax=vmax, kwargs=plot_kwargs)
if k == 'data':
ax.set_ylabel(keyytitlearr[j])
for pos in ['top', 'bottom', 'left', 'right']:
ax.spines[pos].set_visible(False)
ax.set_xticks([])
ax.set_yticks([])
if j == 0:
ax.set_title(xt)
####
if k == 'residual':
med = np.median(im[obs.data.mask])
rms = np.std(im[obs.data.mask])
median_str = r"${}".format(annstr_arr[j])+r"_{med}="+r"{:0.1f}".format(med)+r"$"
scatter_str = r"${}".format(annstr_arr[j])+r"_{rms}="+r"{:0.1f}".format(rms)+r"$"
ax.annotate(median_str,
(0.01,-0.05), xycoords='axes fraction',
ha='left', va='top', fontsize=8)
ax.annotate(scatter_str,
(0.99,-0.05), xycoords='axes fraction',
ha='right', va='top', fontsize=8)
cbar = ax.cax.colorbar(imax)
cbar.ax.tick_params(labelsize=8)
# Figure title is the observation name
f.suptitle(obs.name)
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
def plot_data_model_comparison_3D(obs, model,
theta = None,
fileout=None,
dscale=None,
symmetric_residuals=True,
show_1d_apers = False,
max_residual=100.,
inst_corr = True,
vcrop=False,
vcrop_value=800.,
overwrite=False,
moment=False,
remove_shift=False,
fill_mask=False,
show_contours=False,
show_ruler=True,
ruler_loc='lowerleft',
**plot_kwargs):
if fileout is not None:
fileout_in = fileout
f_r = fileout.split('.')
f_base = ".".join(f_r[:-1])
fileout_aperture = "{}_apertures.{}".format(f_base, f_r[-1])
fileout_spaxel = "{}_spaxels.{}".format(f_base, f_r[-1])
fileout_channel = "{}_channels.{}".format(f_base, f_r[-1])
else:
fileout_aperture = fileout_spaxel = fileout_channel = None
# Check for existing file:
if (not overwrite) and (fileout is not None):
for f in [fileout, fileout_aperture, fileout_spaxel, fileout_channel]:
if os.path.isfile(f):
logger.warning("overwrite={} & File already exists! Will not save file: {} \n ".format(overwrite, f))
return None
plot_3D_extracted_to_1D_2D(obs, model, fileout=fileout, dscale=dscale,
symmetric_residuals=symmetric_residuals,
max_residual=max_residual, show_1d_apers=show_1d_apers,
inst_corr=True, remove_shift=False,
vcrop=vcrop, vcrop_value=vcrop_value,
overwrite=overwrite, fill_mask=fill_mask,
show_ruler=show_ruler, ruler_loc=ruler_loc,
show_contours=show_contours, **plot_kwargs)
plot_spaxel_compare_3D_cubes(obs, fileout=fileout_spaxel, typ='all',
show_model=True, overwrite=overwrite)
plot_aperture_compare_3D_cubes(obs, model, fileout=fileout_aperture,
fill_mask=fill_mask, overwrite=overwrite)
plot_channel_maps_3D_cube(obs, model, show_data=True, show_model=True,
show_residual=True, fileout=fileout_channel,
vbounds = [-450., 450.], delv = 100.,
vbounds_shift=True, cmap=cmap_greys,
overwrite=overwrite)
return None
#############################################################
#############################################################
def plot_3D_extracted_to_1D_2D(obs_in, model_in,
fileout=None,
dscale=None,
symmetric_residuals=True, max_residual=100.,
show_1d_apers=False, inst_corr=None,
xshift = None,
yshift = None,
vcrop=False,
vcrop_value=800.,
remove_shift=False,
overwrite=False,
fill_mask=False,
show_contours=False,
show_ruler=True,
ruler_loc='lowerleft',
**plot_kwargs):
obs = copy.deepcopy(obs_in)
model = copy.deepcopy(model_in)
print("plot_model_multid: ndim=3: moment={}".format(obs.instrument.moment))
obs_extract, model = extract_1D_2D_obs_from_cube(obs, model, inst_corr=True,
fill_mask=fill_mask)
if fileout is not None:
f_r = fileout.split('.')
f_base = ".".join(f_r[:-1])
fileout_1D = "{}_{}.{}".format(f_base, "extract_1D", f_r[-1])
fileout_2D = "{}_{}.{}".format(f_base, "extract_2D", f_r[-1])
else:
fileout_1D = fileout_2D = None
# Data haven't actually been corrected for instrument LSF yet
# (Note: 1D/2D *models* will be corrected for LSF during plotting,
# based on the data['inst_corr'] setting)
if obs_extract['extract_1D'].data.data['inst_corr']:
inst_corr_sigma = obs_extract['extract_1D'].instrument.lsf.dispersion.to(u.km/u.s).value
disp_prof_1D = np.sqrt(obs_extract['extract_1D'].data.data['dispersion']**2 - inst_corr_sigma**2 )
disp_prof_1D[~np.isfinite(disp_prof_1D)] = 0.
obs_extract['extract_1D'].data.data['dispersion'] = disp_prof_1D
if 'filled_mask_data' in obs_extract['extract_1D'].data.__dict__.keys():
disp_prof_1D = np.sqrt(obs_extract['extract_1D'].data.filled_mask_data.data['dispersion']**2 - inst_corr_sigma**2 )
disp_prof_1D[~np.isfinite(disp_prof_1D)] = 0.
obs_extract['extract_1D'].data.filled_mask_data.data['dispersion'] = disp_prof_1D
if obs_extract['extract_2D'].data.data['inst_corr']:
inst_corr_sigma = obs_extract['extract_2D'].instrument.lsf.dispersion.to(u.km/u.s).value
im = obs_extract['extract_2D'].data.data['dispersion'].copy()
im = np.sqrt(im ** 2 - inst_corr_sigma ** 2)
im[~np.isfinite(im)] = 0.
obs_extract['extract_2D'].data.data['dispersion'] = im
plot_data_model_comparison_1D(obs_extract['extract_1D'], fileout=fileout_1D,
overwrite=overwrite, **plot_kwargs)
plot_data_model_comparison_2D(obs_extract['extract_2D'], model,
fileout=fileout_2D, show_contours=show_contours,
show_ruler=show_ruler, ruler_loc=ruler_loc,
overwrite=overwrite, **plot_kwargs)
return None
#############################################################
def plot_aperture_compare_3D_cubes(obs, model, datacube=None, errcube=None,
modelcube=None, mask=None,
fileout=None,
slit_width=None, slit_pa=None,
aper_dist=None, fill_mask=False,
skip_fits=True, overwrite=False):
#############################################################
if datacube is None:
datacube = obs.data.data
if errcube is None:
errcube = obs.data.error
if modelcube is None:
modelcube = obs.model_data.data
if mask is None:
mask = obs.data.mask
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
######################################
if slit_width is None:
try:
slit_width = obs.instrument.beam.major.to(u.arcsec).value
except:
slit_width = obs.instrument.beam.major_fwhm.to(u.arcsec).value
if slit_pa is None:
slit_pa = model.geometries[obs.name].pa.value
if mask is None:
mask = obs.data.mask.copy()
pixscale = obs.instrument.pixscale.value
rpix = slit_width/pixscale/2.
# Aper centers: pick roughly number fitting into size:
nx = datacube.shape[2]
ny = datacube.shape[1]
try:
center_pixel = (obs.mod_options.xcenter + model.geometries[obs.name].xshift.value,
obs.mod_options.ycenter + model.geometries[obs.name].yshift.value)
except:
#center_pixel = (int(nx / 2.) + model.geometries[obs.name].geometry.xshift,
# int(ny / 2.) + model.geometries[obs.name].yshift)
center_pixel = (int(nx / 2.) + model.geometries[obs.name].xshift.value,
int(ny / 2.) + model.geometries[obs.name].yshift.value)
aper_centers_arcsec = aper_centers_arcsec_from_cube(datacube, obs, model,
mask=mask, slit_width=slit_width, slit_pa=slit_pa,
aper_dist=aper_dist, fill_mask=fill_mask)
#############################################################
specarr = datacube.spectral_axis.to(u.km/u.s).value
apertures = CircApertures(rarr=aper_centers_arcsec, slit_PA=slit_pa, rpix=rpix,
nx=nx, ny=ny, center_pixel=center_pixel, pixscale=pixscale, moment=False)
if not skip_fits:
data_scaled = datacube.unmasked_data[:].value
if errcube is not None:
ecube = errcube.unmasked_data[:].value * mask
else:
ecube = None
model_scaled = modelcube.unmasked_data[:].value
aper_centers, flux1d, vel1d, disp1d = apertures.extract_1d_kinematics(spec_arr=specarr,
cube=data_scaled, mask=mask, err=ecube,
center_pixel = center_pixel, pixscale=pixscale)
#####
aper_centers_mod, flux1d_mod, vel1d_mod, disp1d_mod = apertures.extract_1d_kinematics(spec_arr=specarr,
cube=model_scaled, mask=mask, err=None,
center_pixel = center_pixel, pixscale=pixscale)
apertures_mom = CircApertures(rarr=aper_centers_arcsec, slit_PA=slit_pa, rpix=rpix,
nx=nx, ny=ny, center_pixel=center_pixel, pixscale=pixscale, moment=True)
aper_centers_mod_mom, flux1d_mod_mom, vel1d_mod_mom, disp1d_mod_mom = apertures_mom.extract_1d_kinematics(spec_arr=specarr,
cube=model_scaled, mask=mask, err=None,
center_pixel = center_pixel, pixscale=pixscale)
aper_centers2, flux1d2, vel1d2, disp1d2 = apertures_mom.extract_1d_kinematics(spec_arr=specarr,
cube=data_scaled, mask=mask, err=ecube,
center_pixel = center_pixel, pixscale=pixscale)
######################################
# Setup plot:
nrows = int(np.round(np.sqrt(len(aper_centers_arcsec))))
ncols = int(np.ceil(len(aper_centers_arcsec)/(1.*nrows)))
padx = 0.25
pady = 0.25
xextra = 0.15
yextra = 0.
scale = 2.5
f = plt.figure()
f.set_size_inches((ncols+(ncols-1)*padx+xextra)*scale, (nrows+pady+yextra)*scale)
suptitle = '{}: ndim={}'.format(obs.name,obs.instrument.ndim)
gs = gridspec.GridSpec(nrows, ncols, wspace=padx, hspace=pady)
axes = []
for i in range(nrows):
for j in range(ncols):
axes.append(plt.subplot(gs[i,j]))
#############################################################
# for each ap:
for k in range(len(axes)):
ax = axes[k]
if k < len(aper_centers_arcsec):
_, datarr, errarr = apertures.apertures[k].extract_aper_spec(spec_arr=specarr,
cube=datacube, mask=mask, err=errcube, skip_specmask=True)
_, modarr, _ = apertures.apertures[k].extract_aper_spec(spec_arr=specarr,
cube=modelcube, mask=mask, skip_specmask=True)
_, maskarr, _ = apertures.apertures[k].extract_aper_spec(spec_arr=specarr,
cube=mask, skip_specmask=True)
maskarr[maskarr>0] = 1.
if not skip_fits:
gmod_flux2=flux1d_mod_mom[k]
gmod_vel2=vel1d_mod_mom[k]
gmod_disp2=disp1d_mod_mom[k]
gmod_flux2 = gmod_vel2 = gmod_disp2 = None
gdata_flux2=flux1d2[k]
gdata_vel2=vel1d2[k]
gdata_disp2=disp1d2[k]
gdata_flux=flux1d[k]
gdata_vel=vel1d[k]
gdata_disp=disp1d[k]
gmod_flux=flux1d_mod[k]
gmod_vel=vel1d_mod[k]
gmod_disp=disp1d_mod[k]
else:
gdata_flux = gdata_vel = gdata_disp = None
gdata_flux2 = gdata_vel2 = gdata_disp2 = None
gmod_flux = gmod_vel = gmod_disp = None
gmod_flux2 = gmod_vel2 = gmod_disp2 = None
ax = plot_spaxel_fit(specarr, datarr, maskarr, err=errarr,
gdata_flux=gdata_flux, gdata_vel=gdata_vel, gdata_disp=gdata_disp,
gdata_flux2=gdata_flux2, gdata_vel2=gdata_vel2, gdata_disp2=gdata_disp2,
model=modarr, gmod_flux=gmod_flux, gmod_vel=gmod_vel, gmod_disp=gmod_disp,
gmod_flux2=gmod_flux2, gmod_vel2=gmod_vel2, gmod_disp2=gmod_disp2,
ax=ax)
ax.annotate('Ap {}'.format(k),
(0.02,0.98), xycoords='axes fraction',
ha='left', va='top', fontsize=8)
else:
ax.set_axis_off()
#############################################################
# Save to file:
f.suptitle(suptitle, fontsize=16, y=0.925)
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
def plot_spaxel_compare_3D_cubes(obs, datacube=None, errcube=None,
modelcube=None, mask=None,
fileout=None,
typ='all',
show_model=True,
skip_masked=False,
skip_fits=True,
overwrite=False):
if typ.strip().lower() not in ['all', 'diag']:
raise ValueError("typ={} not recognized for `plot_spaxel_compare_3D_cubes`!".format(typ))
# Clean up:
typ = typ.strip().lower()
if datacube is None:
datacube = obs.data.data
if errcube is None:
errcube = obs.data.error
if mask is None:
mask = np.array(obs.data.mask, dtype=float)
if show_model:
if modelcube is None:
try:
modelcube = obs.model_data.data
except:
# Skip showing model
show_model = False
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
# Aper centers: pick roughly number fitting into size:
nx = datacube.shape[2]
ny = datacube.shape[1]
npix = np.max([nx,ny])
specarr = datacube.spectral_axis.to(u.km/u.s).value
######################################
# Setup plot:
if typ == 'all':
rowinds = np.where(np.sum(np.sum(mask,axis=0),axis=1)>0)[0]
colinds = np.where(np.sum(np.sum(mask,axis=0),axis=0)>0)[0]
nrows = len(rowinds)
ncols = len(colinds)
elif typ == 'diag':
nrows = int(np.round(np.sqrt(npix)))
ncols = int(np.ceil(npix/(1.*nrows)))
rowinds = np.arange(nrows*ncols)
colinds = np.arange(nrows*ncols)
padx = 0.25
pady = 0.25
xextra = 0.15
yextra = 0.
scale = 2.5
f = plt.figure()
figsize = ((ncols+(ncols-1)*padx+xextra)*scale,
(nrows+(nrows-1)*pady+yextra)*scale)
f.set_size_inches((ncols+(ncols-1)*padx+xextra)*scale,
(nrows+(nrows-1)*pady+yextra)*scale)
suptitle = '{}: ndim={}'.format(obs.name, obs.instrument.ndim)
gs = gridspec.GridSpec(nrows, ncols, wspace=padx, hspace=pady)
axes = []
for i in range(nrows):
for j in range(ncols):
if typ == 'all':
# invert rows:
ii = nrows - 1 - i
else:
ii = i
axes.append(plt.subplot(gs[ii,j]))
#############################################################
if not skip_fits:
if show_model:
mom0_mod = modelcube.moment0().to(u.km/u.s).value
mom1_mod = modelcube.moment1().to(u.km/u.s).value
mom2_mod = modelcube.linewidth_sigma().to(u.km/u.s).value
maskbool = np.array(mask, dtype=bool)
datacube_masked = datacube.with_mask(maskbool)
mom0_dat = datacube_masked.moment0().to(u.km/u.s).value
mom1_dat = datacube_masked.moment1().to(u.km/u.s).value
mom2_dat = datacube_masked.linewidth_sigma().to(u.km/u.s).value
# for each spax:
k = -1
for i in rowinds:
for j in colinds:
skip = False
if typ == 'diag':
if (i == j):
k += 1
else:
skip = True
elif typ == 'all':
k += 1
if (k < len(axes)) & (not skip) & (i < npix):
ax = axes[k]
datarr = datacube[:,i,j].value
maskarr = mask[:,i,j]
errarr = errcube[:,i,j].value
if show_model:
modarr = modelcube[:,i,j].value
else:
modarr = None
do_plot = True
if skip_masked:
if (maskarr.max() <= 0):
do_plot = False
if ((typ == 'all') & (do_plot)) | (typ == 'diag'):
if not skip_fits:
if show_model:
flux1d_mod_mom = mom0_mod[i,j]
vel1d_mod_mom = mom1_mod[i,j]
disp1d_mod_mom = mom2_mod[i,j]
best_fit = gaus_fit_sp_opt_leastsq(specarr, modarr, mom0_mod[i,j],
mom1_mod[i,j], mom2_mod[i,j])
flux1d_mod = best_fit[0] * np.sqrt(2 * np.pi) * best_fit[2]
vel1d_mod = best_fit[1]
disp1d_mod = best_fit[2]
else:
flux1d_mod_mom = None
vel1d_mod_mom = None
disp1d_mod_mom = None
flux1d_mod = None
vel1d_mod = None
disp1d_mod = None
maskarr_bool = np.array(maskarr, dtype=bool)
flux1d2 = mom0_dat[i,j]
vel1d2 = mom1_dat[i,j]
disp1d2 = mom2_dat[i,j]
try:
best_fit = gaus_fit_apy_mod_fitter(specarr[maskarr_bool], datarr[maskarr_bool],
mom0_dat[i,j], mom1_dat[i,j], mom2_dat[i,j], yerr=errarr[maskarr_bool])
except:
best_fit = [np.NaN, np.NaN, np.NaN]
flux1d = best_fit[0] * np.sqrt(2 * np.pi) * best_fit[2]
vel1d = best_fit[1]
disp1d = best_fit[2]
gmod_flux2=flux1d_mod_mom
gmod_vel2=vel1d_mod_mom
gmod_disp2=disp1d_mod_mom
else:
flux1d = vel1d = disp1d = None
flux1d2 = vel1d2 = disp1d2 = None
flux1d_mod = vel1d_mod = disp1d_mod = None
gmod_flux2 = gmod_vel2 = gmod_disp2 = None
ax = plot_spaxel_fit(specarr, datarr, maskarr, err=errarr,
gdata_flux=flux1d, gdata_vel=vel1d, gdata_disp=disp1d,
gdata_flux2=flux1d2, gdata_vel2=vel1d2, gdata_disp2=disp1d2,
model=modarr, gmod_flux=flux1d_mod, gmod_vel=vel1d_mod,
gmod_disp=disp1d_mod,
gmod_flux2=gmod_flux2, gmod_vel2=gmod_vel2,
gmod_disp2=gmod_disp2,
ax=ax)
ax.annotate('Pix ({},{})'.format(j,i),
(0.02,0.98), xycoords='axes fraction',
ha='left', va='top', fontsize=8)
if (maskarr.max() <= 0):
ax.set_facecolor('#f0f0f0')
else:
ax.set_axis_off()
elif (k < len(axes)) & (not skip) & (k >= npix):
ax = axes[k]
ax.set_axis_off()
#############################################################
# Save to file:
if suptitle is not None:
yoff = 0.105 - 0.01*(36.875-f.get_size_inches()[1])/18.75
ytitlepos = 1.-yoff
if f.get_size_inches()[1] < 20.:
ytitlefontsize = 20
else:
ytitlefontsize = 30
f.suptitle(suptitle, fontsize=ytitlefontsize, y=ytitlepos)
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
#############################################################
def plot_channel_maps_3D_cube(obs, model, show_data=True,
show_model=True, show_residual=True,
vbounds = [-450., 450.], delv = 100.,
vbounds_shift=True,
cmap=cmap_greys, cmap_resid=cmap_seismic,
fileout=None, overwrite=False):
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
if (not show_data) | (not show_model):
show_residual = False
vbounds = np.array(vbounds)
if vbounds_shift:
vbounds += model.geometries[obs.name].vel_shift.value
v_slice_lims_arr = np.arange(vbounds[0], vbounds[1]+delv, delv)
#################################################
# center slice: flux limits:
ind = int(np.round((len(v_slice_lims_arr)-2)/2.))
v_slice_lims = v_slice_lims_arr[ind:ind+2]
if show_data:
subcube = obs.data.data.spectral_slab(v_slice_lims[0]*u.km/u.s, v_slice_lims[1]*u.km/u.s)
im = subcube.moment0().value * obs.data.mask
else:
subcube = obs.model_data.data.spectral_slab(v_slice_lims[0]*u.km/u.s, v_slice_lims[1]*u.km/u.s)
im = subcube.moment0().value * obs.model_data.mask
flux_lims = [im.min(), im.max()]
fac = 1.
immax = np.max(np.abs(im))
flux_lims_resid = [-fac*immax, fac*immax]
#################################################
f = plt.figure()
scale = 3.
show_multi = True
if show_residual:
n_cols = 3
n_rows = len(v_slice_lims_arr)-1
ind_dat = 0
ind_mod = 1
ind_resid = 2
elif (show_data) & (show_model):
n_cols = 2
n_rows = len(v_slice_lims_arr)-1
ind_dat = 0
ind_mod = 1
ind_resid = None
else:
n_cols = int(np.ceil(np.sqrt(len(v_slice_lims_arr)-1)))
n_rows = int(np.ceil((len(v_slice_lims_arr)-1.)/(1.*n_cols)))
show_multi = False
if show_data:
ind_dat = 0
ind_mod = ind_resid = None
else:
ind_mod = 0
ind_dat = ind_resid = None
wspace_outer = 0.
wspace = hspace = padfac = 0.1
fac = 1.
f.set_size_inches(fac*scale*n_cols+(n_cols-1)*scale*padfac,scale*n_rows+(n_rows-1)*scale*padfac)
gs = gridspec.GridSpec(n_rows,n_cols, wspace=wspace, hspace=hspace)
if show_multi:
if show_data:
axes_dat = []
if show_model:
axes_mod = []
if show_residual:
axes_resid = []
for j in range(n_rows):
if show_data:
axes_dat.append(plt.subplot(gs[j,ind_dat]))
if show_model:
axes_mod.append(plt.subplot(gs[j,ind_mod]))
if show_residual:
axes_resid.append(plt.subplot(gs[j,ind_resid]))
else:
if show_data:
axes_dat = []
if show_model:
axes_mod = []
for j in range(n_rows):
for i in range(n_cols):
if show_data:
axes_dat.append(plt.subplot(gs[j,i]))
if show_model:
axes_mod.append(plt.subplot(gs[j,i]))
types = []
axes_stack = []
if show_data:
types.append('data')
axes_stack.append(axes_dat)
if show_model:
types.append('model')
axes_stack.append(axes_mod)
if show_residual:
types.append('residual')
axes_stack.append(axes_resid)
center = np.array([(obs.model_data.data.shape[2]-1.)/2., (obs.model_data.data.shape[1]-1.)/2.])
center[0] += model.geometries[obs.name].xshift.value
center[1] += model.geometries[obs.name].yshift.value
k = -1
for ii in range(n_rows):
if show_multi:
k += 1
v_slice_lims = v_slice_lims_arr[k:k+2]
for j in range(n_cols):
if show_multi:
ind_stack = j
else:
k += 1
ind_stack = 0
v_slice_lims = v_slice_lims_arr[k:k+2]
##
ax = axes_stack[ind_stack][k]
typ = types[ind_stack]
flims = flux_lims
cmap_tmp = cmap
residual=False
if typ == 'data':
cube = obs.data.data*obs.data.mask
color_contours='blue'
elif typ == 'model':
cube = obs.model_data.data*obs.model_data.mask
color_contours='red'
elif typ == 'residual':
cube = obs.data.data*obs.data.mask - obs.model_data.data*obs.model_data.mask
flims = flux_lims_resid
cmap_tmp = cmap_resid
color_contours='black'
residual=True
ax = plot_channel_slice(ax=ax,speccube=cube, center=center,
v_slice_lims=v_slice_lims, flux_lims=flims,
cmap=cmap_tmp, color_contours=color_contours,
residual=residual)
###########################
label_str = None
if (show_multi & (ii == 0)):
if (j==0):
label_str = "{}: {}".format(obs.name, typ.capitalize())
else:
label_str = typ.capitalize()
elif ((not show_multi) & (k == 0)):
label_str = "{}: {}".format(obs.name, typ.capitalize())
if label_str is not None:
ax.set_title(label_str, fontsize=14)
###########################
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
#############################################################
def plot_spaxels_cube(cube=None, hdr=None, mask=None,
fname=None,
typ='all',
skip_masked=False,
fileout=None,
overwrite=False):
if typ.strip().lower() not in ['all', 'diag']:
raise ValueError("typ={} not recognized for `plot_spaxels_cube`!".format(typ))
# Clean up:
typ = typ.strip().lower()
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
if ((cube is None) | (hdr is None)) and (fname is None):
raise ValueError("Either 'fname' must be specified,"
"or both 'cube' and 'hdr' must be set!")
import astropy.units as u
from astropy.wcs import WCS
from spectral_cube import SpectralCube
if ((cube is None) | (hdr is None)):
datcube = SpectralCube.read(fname)
else:
if (hdr['SPEC_TYPE'] == 'VOPT'):
spec_unit = u.km/u.s
elif (hdr['SPEC_TYPE'] == 'WAVE'):
spec_unit = u.Angstrom
w = WCS(header=hdr)
datcube = SpectralCube(data=cube, wcs=w, mask=mask).with_spectral_unit(spec_unit)
#################################################
# Aper centers: pick roughly number fitting into size:
nx = datcube.shape[2]
ny = datcube.shape[1]
npix = np.max([nx,ny])
specarr = datcube.spectral_axis.to(u.km/u.s).value
######################################
# Setup plot:
if typ == 'all':
if mask is not None:
rowinds = np.where(np.sum(np.sum(mask,axis=0),axis=1)>0)[0]
colinds = np.where(np.sum(np.sum(mask,axis=0),axis=0)>0)[0]
else:
rowinds = np.where(np.sum(np.sum(np.abs(datcube),axis=0),axis=1)>0)[0]
colinds = np.where(np.sum(np.sum(np.abs(datcube),axis=0),axis=0)>0)[0]
nrows = len(rowinds)
ncols = len(colinds)
elif typ == 'diag':
nrows = int(np.round(np.sqrt(npix)))
ncols = int(np.ceil(npix/(1.*nrows)))
rowinds = np.arange(nrows*ncols)
colinds = np.arange(nrows*ncols)
padx = 0.25
pady = 0.25
xextra = 0.15
yextra = 0.
scale = 2.5
f = plt.figure()
figsize = ((ncols+(ncols-1)*padx+xextra)*scale,
(nrows+(nrows-1)*pady+yextra)*scale)
f.set_size_inches((ncols+(ncols-1)*padx+xextra)*scale,
(nrows+(nrows-1)*pady+yextra)*scale)
suptitle = None
gs = gridspec.GridSpec(nrows, ncols, wspace=padx, hspace=pady)
axes = []
for i in range(nrows):
for j in range(ncols):
if typ == 'all':
# invert rows:
ii = nrows - 1 - i
else:
ii = i
axes.append(plt.subplot(gs[ii,j]))
#############################################################
# for each spax:
k = -1
for i in rowinds:
for j in colinds:
skip = False
if typ == 'diag':
if (i == j):
k += 1
else:
skip = True
elif typ == 'all':
k += 1
if (k < len(axes)) & (not skip) & (i < npix):
ax = axes[k]
datarr = datcube[:,i,j].value
if mask is not None:
maskarr = mask[:,i,j]
else:
maskarr = specarr * 0. + 1.
do_plot = True
if skip_masked:
if (maskarr.max() <= 0):
do_plot = False
if ((typ == 'all') & (do_plot)) | (typ == 'diag'):
ax.plot(specarr, datarr, color='black', ls='-', lw=1., alpha=0.5, zorder=1.)
ax.plot(specarr, datarr*maskarr, color='black',lw=1.5, zorder=1.)
ax.axhline(y=0., ls='--', color='grey', alpha=0.5, zorder=-20.)
xlim = ax.get_xlim()
xrange = xlim[1]-xlim[0]
if xrange >= 1000.:
xmajloc = 500.
xminloc = 100.
elif xrange >= 500.:
xmajloc = 200.
xminloc = 50.
elif xrange >= 250.:
xmajloc = 100.
xminloc = 20.
elif xrange >= 100.:
xmajloc = 50.
xminloc = 10.
elif xrange >= 50.:
xmajloc = 10.
xminloc = 2.
elif xrange >= 10.:
xmajloc = 2.
xminloc = 0.5
else:
xmajloc = None
xminloc = None
if xmajloc is not None:
ax.xaxis.set_major_locator(MultipleLocator(xmajloc))
ax.xaxis.set_minor_locator(MultipleLocator(xminloc))
ax.annotate('Pix ({},{})'.format(j,i),
(0.02,0.98), xycoords='axes fraction',
ha='left', va='top', fontsize=8)
if (maskarr.max() <= 0):
ax.set_facecolor('#f0f0f0')
else:
ax.set_axis_off()
elif (k < len(axes)) & (not skip) & (k >= npix):
ax = axes[k]
ax.set_axis_off()
#############################################################
# Save to file:
if suptitle is not None:
yoff = 0.105 - 0.01*(36.875-f.get_size_inches()[1])/18.75
ytitlepos = 1.-yoff
if f.get_size_inches()[1] < 20.:
ytitlefontsize = 20
else:
ytitlefontsize = 30
f.suptitle(suptitle, fontsize=ytitlefontsize, y=ytitlepos)
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
#############################################################
def plot_spaxel_fit(specarr, data, mask, err=None,
gdata_flux=None, gdata_vel=None, gdata_disp=None,
gdata_flux2=None, gdata_vel2=None, gdata_disp2=None,
model=None,
gmod_flux=None, gmod_vel=None, gmod_disp=None,
gmod_flux2=None, gmod_vel2=None, gmod_disp2=None,
ax=None):
returnax = True
if ax is None:
returnax = False
ax = plt.subplot(111)
ax.plot(specarr, data*mask, color='black', marker='o', ms=4., zorder=1.)
if (mask.max() > 0):
ylim = ax.get_ylim()
ax.plot(specarr, data, color='black', marker='o', ms=4., mfc='None', ls='None', alpha=0.5, zorder=0.)
if (mask.max() > 0):
ax.set_ylim(ylim)
if gdata_flux is not None:
try:
gdata_A = gdata_flux / ( np.sqrt(2 * np.pi) * gdata_disp)
ax.plot(specarr, gdata_A*np.exp(-((specarr-gdata_vel)**2/(2.*gdata_disp**2))),
color='turquoise',zorder=10., lw=0.5)
except:
pass
if gdata_flux2 is not None:
try:
gdata_A2 = gdata_flux2 / ( np.sqrt(2 * np.pi) * gdata_disp2)
ax.plot(specarr, gdata_A2*np.exp(-((specarr-gdata_vel2)**2/(2.*gdata_disp2**2))),
color='tab:green', ls='--', zorder=5., lw=0.5)
except:
pass
if model is not None:
ax.plot(specarr, model, color='red', ls='-', lw=1., alpha=0.5, zorder=1.)
ax.plot(specarr, model*mask, color='red',lw=1.5, zorder=1.)
if gmod_flux is not None:
try:
gmod_A = gmod_flux / ( np.sqrt(2 * np.pi) * gmod_disp)
ax.plot(specarr, gmod_A*np.exp(-((specarr-gmod_vel)**2/(2.*gmod_disp**2))),
color='orange', ls='--', zorder=10., lw=0.5)
except:
pass
if gmod_flux2 is not None:
try:
gmod_A2 = gmod_flux2 / ( np.sqrt(2 * np.pi) * gmod_disp2)
ax.plot(specarr, gmod_A2*np.exp(-((specarr-gmod_vel2)**2/(2.*gmod_disp2**2))),
color='purple', ls=':', zorder=10., lw=0.5)
except:
pass
if err is not None:
ylim = ax.get_ylim()
ax.errorbar(specarr, data, xerr=None, yerr=err, alpha=0.25, capsize=0.,
marker=None, ls='None', ecolor='k', zorder=-1.)
ax.set_ylim(ylim)
ax.axhline(y=0., ls='--', color='grey', alpha=0.5, zorder=-20.)
xlim = ax.get_xlim()
xrange = xlim[1]-xlim[0]
if xrange >= 1000.:
xmajloc = 500.
xminloc = 100.
elif xrange >= 500.:
xmajloc = 200.
xminloc = 50.
elif xrange >= 250.:
xmajloc = 100.
xminloc = 20.
elif xrange >= 100.:
xmajloc = 50.
xminloc = 10.
elif xrange >= 50.:
xmajloc = 10.
xminloc = 2.
elif xrange >= 10.:
xmajloc = 2.
xminloc = 0.5
else:
xmajloc = None
xminloc = None
if xmajloc is not None:
ax.xaxis.set_major_locator(MultipleLocator(xmajloc))
ax.xaxis.set_minor_locator(MultipleLocator(xminloc))
if returnax:
return ax
else:
return None
#############################################################
def plot_channel_maps_cube(cube=None, hdr=None, mask=None,
fname=None,
vbounds = [-450., 450.],
delv = 100.,
cmap=cmap_greys,
fileout=None,
overwrite=False):
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
if ((cube is None) | (hdr is None)) and (fname is None):
raise ValueError("Either 'fname' must be specified,"
"or both 'cube' and 'hdr' must be set!")
import astropy.units as u
from astropy.wcs import WCS
from spectral_cube import SpectralCube
if ((cube is None) | (hdr is None)):
datcube = SpectralCube.read(fname)
else:
if (hdr['SPEC_TYPE'] == 'VOPT'):
spec_unit = u.km/u.s
elif (hdr['SPEC_TYPE'] == 'WAVE'):
spec_unit = u.Angstrom
w = WCS(header=hdr)
datcube = SpectralCube(data=cube, wcs=w, mask=mask).with_spectral_unit(spec_unit)
vbounds = np.array(vbounds)
v_slice_lims_arr = np.arange(vbounds[0], vbounds[1]+delv, delv)
#################################################
# center slice: flux limits:
ind = int(np.round((len(v_slice_lims_arr)-2)/2.))
v_slice_lims = v_slice_lims_arr[ind:ind+2]
subcube = datcube.spectral_slab(v_slice_lims[0]*u.km/u.s, v_slice_lims[1]*u.km/u.s)
im = subcube.moment0().value
flux_lims = [im.min(), im.max()]
fac = 1.
immax = np.max(np.abs(im))
flux_lims_resid = [-fac*immax, fac*immax]
#################################################
f = plt.figure()
scale = 3.
n_cols = int(np.ceil(np.sqrt(len(v_slice_lims_arr)-1)))
n_rows = int(np.ceil((len(v_slice_lims_arr)-1.)/(1.*n_cols)))
wspace_outer = 0.
wspace = hspace = padfac = 0.1
fac = 1.
f.set_size_inches(fac*scale*n_cols+(n_cols-1)*scale*padfac,scale*n_rows+(n_rows-1)*scale*padfac)
gs = gridspec.GridSpec(n_rows,n_cols, wspace=wspace, hspace=hspace)
axes = []
for j in range(n_rows):
for i in range(n_cols):
axes.append(plt.subplot(gs[j,i]))
center = np.array([(datcube.shape[2]-1.)/2., (datcube.shape[1]-1.)/2.])
k = -1
for ii in range(n_rows):
for j in range(n_cols):
k += 1
ax = axes[k]
if k > (len(v_slice_lims_arr)-2):
ax.set_axis_off()
else:
v_slice_lims = v_slice_lims_arr[k:k+2]
flims = flux_lims
cmap_tmp = cmap
color_contours='red'
ax = plot_channel_slice(ax=ax,speccube=datcube, center=center,
v_slice_lims=v_slice_lims, flux_lims=flims,
cmap=cmap_tmp, color_contours=color_contours,
residual=False)
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
def plot_channel_slice(ax=None, speccube=None, v_slice_lims=None, flux_lims=None,
center=None, show_pix_coords=False,
cmap=cmap_greys,
color_contours='red', color_center='cyan',
residual=False):
if ax is None:
ax = plt.gca()
if v_slice_lims is None:
v_slice_lims = [-50., 50.]
subcube = speccube.spectral_slab(v_slice_lims[0]*u.km/u.s, v_slice_lims[1]*u.km/u.s)
im = subcube.moment0().value
if center is None:
center = ((im.shape[1]-1.)/2., (im.shape[0]-1.)/2.)
int_mode = "nearest"
origin = 'lower'
if flux_lims is not None:
vmin = flux_lims[0]
vmax = flux_lims[1]
else:
vmin = None
vmax = None
imax = ax.imshow(im, cmap=cmap, interpolation=int_mode,
vmin=vmin, vmax=vmax, origin=origin)
sigma_levels = [1.,2.,3.]
levels = 1.0 - np.exp(-0.5 * np.array(sigma_levels) ** 2)
lw_contour = 1.5
rgba_color = mplcolors.colorConverter.to_rgba(color_contours)
color_levels = [list(rgba_color) for l in levels]
lw_arr = []
for ii, l in enumerate(levels):
color_levels[ii][-1] *= float(ii+1) / (len(levels))
lw_arr.append(lw_contour * float(ii+1+len(levels)*0.5) / (len(levels)*1.5))
contour_kwargs = {'colors': color_levels, 'linewidths': lw_arr, 'linestyles': '-'}
#################################################
# Syntax taken from corner/core.py
imflat = im.flatten()
imtmp = im.copy()
if residual:
imflat = np.abs(imflat)
imtmp = np.abs(imtmp)
inds = np.argsort(imflat)[::-1]
imflat = imflat[inds]
sm = np.cumsum(imflat)
sm /= sm[-1]
contour_levels = np.empty(len(levels))
for i, v0 in enumerate(levels):
try:
contour_levels[i] = imflat[sm<=v0][-1]
except:
contour_levels[i] = imflat[0]
contour_levels.sort()
m = np.diff(contour_levels) == 0
if np.any(m):
print("Too few points to create valid contours")
while np.any(m):
contour_levels[np.where(m)[0][0]] *= 1.0 - 1e-4
m = np.diff(contour_levels) == 0
contour_levels.sort()
ax.contour(imtmp, contour_levels, **contour_kwargs)
#################################################
ax.plot(center[0], center[1], '+', mew=1., ms=10., c=color_center)
xlim = ax.get_xlim()
ylim = ax.get_ylim()
ybase_offset = 0.035
x_base = xlim[0] + (xlim[1]-xlim[0])*0.075
y_base = ylim[0] + (ylim[1]-ylim[0])*(ybase_offset+0.075)
string = r'$[{:0.1f},{:0.1f}]$'.format(v_slice_lims[0], v_slice_lims[1])
ax.annotate(string, xy=(x_base, y_base),
xycoords="data",
xytext=(0,0),
color='black',
textcoords="offset points", ha="left", va="center",
fontsize=8)
if not show_pix_coords:
ax.set_xticks([])
ax.set_yticks([])
return ax
#############################################################
def plot_channel_maps_cube_overlay(obs, model,
vbounds = [-450., 450.], delv = 100.,
vbounds_shift=True,
cmap=cmap_greys, cmap_resid=cmap_seismic,
fileout=None, overwrite=False):
# Check for existing file:
if (not overwrite) and (fileout is not None):
if os.path.isfile(fileout):
logger.warning("overwrite={} & File already exists! Will not save file. \n {}".format(overwrite, fileout))
return None
vbounds = np.array(vbounds)
if vbounds_shift:
vbounds += model.geometries[obs.name].vel_shift.value
v_slice_lims_arr = np.arange(vbounds[0], vbounds[1]+delv, delv)
#################################################
# center slice: flux limits:
ind = int(np.round((len(v_slice_lims_arr)-2)/2.))
v_slice_lims = v_slice_lims_arr[ind:ind+2]
subcube = obs.data.data.spectral_slab(v_slice_lims[0]*u.km/u.s, v_slice_lims[1]*u.km/u.s)
im = subcube.moment0().value * obs.data.mask
flux_lims = [im.min(), im.max()]
fac = 1.
immax = np.max(np.abs(im))
flux_lims_resid = [-fac*immax, fac*immax]
#################################################
f = plt.figure()
scale = 2. #3.
n_cols = int(np.ceil(np.sqrt(len(v_slice_lims_arr)-1)))
n_rows = int(np.ceil((len(v_slice_lims_arr)-1.)/(1.*n_cols)))
wspace_outer = 0.
wspace = hspace = padfac = 0.1
fac = 1.
f.set_size_inches(fac*scale*n_cols+(n_cols-1)*scale*padfac,scale*n_rows+(n_rows-1)*scale*padfac)
gs = gridspec.GridSpec(n_rows,n_cols, wspace=wspace, hspace=hspace)
axes = []
for j in range(n_rows):
for i in range(n_cols):
axes.append(plt.subplot(gs[j,i]))
center = np.array([(obs.model_data.data.shape[2]-1.)/2., (obs.model_data.data.shape[1]-1.)/2.])
center[0] += model.geometries[obs.name].xshift.value
center[1] += model.geometries[obs.name].yshift.value
k = -1
for ii in range(n_rows):
for j in range(n_cols):
k += 1
v_slice_lims = v_slice_lims_arr[k:k+2]
##
ax = axes[k]
flims = flux_lims
cmap_tmp = cmap
cube_dat = obs.data.data*obs.data.mask
cube_mod = obs.model_data.data*obs.model_data.mask
color_contours='red'
ax = plot_channel_slice_diff_contours(ax=ax,speccube=cube_dat,
modcube=cube_mod, center=center,
v_slice_lims=v_slice_lims, flux_lims=flims,
cmap=cmap_tmp, color_contours=color_contours)
# ###########################
# label_str = None
# if (k == 0):
# label_str = "{}: {}".format(obs.name, typ.capitalize())
# if label_str is not None:
# ax.set_title(label_str, fontsize=14)
###########################
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
def plot_channel_slice_diff_contours(ax=None, speccube=None, modcube=None,
v_slice_lims=None, flux_lims=None,
center=None, show_pix_coords=False,
cmap=cmap_greys,
color_contours='red', color_center='cyan'):
if ax is None:
ax = plt.gca()
if v_slice_lims is None:
v_slice_lims = [-50., 50.]
subcube = speccube.spectral_slab(v_slice_lims[0]*u.km/u.s, v_slice_lims[1]*u.km/u.s)
modsubcube = modcube.spectral_slab(v_slice_lims[0]*u.km/u.s, v_slice_lims[1]*u.km/u.s)
im = subcube.moment0().value
immod = modsubcube.moment0().value
if center is None:
center = ((im.shape[1]-1.)/2., (im.shape[0]-1.)/2.)
int_mode = "nearest"
origin = 'lower'
if flux_lims is not None:
vmin = flux_lims[0]
vmax = flux_lims[1]
else:
vmin = None
vmax = None
imax = ax.imshow(im, cmap=cmap, interpolation=int_mode,
vmin=vmin, vmax=vmax, origin=origin)
sigma_levels = [1.,2.,3.]
levels = 1.0 - np.exp(-0.5 * np.array(sigma_levels) ** 2)
lw_contour = 1.5
rgba_color = mplcolors.colorConverter.to_rgba(color_contours)
color_levels = [list(rgba_color) for l in levels]
lw_arr = []
for ii, l in enumerate(levels):
color_levels[ii][-1] *= float(ii+1) / (len(levels))
lw_arr.append(lw_contour * float(ii+1+len(levels)*0.5) / (len(levels)*1.5))
contour_kwargs = {'colors': color_levels, 'linewidths': lw_arr, 'linestyles': '-'}
#################################################
# Syntax taken from corner/core.py
imflat = immod.flatten()
imtmp = immod.copy()
inds = np.argsort(imflat)[::-1]
imflat = imflat[inds]
sm = np.cumsum(imflat)
sm /= sm[-1]
contour_levels = np.empty(len(levels))
for i, v0 in enumerate(levels):
try:
contour_levels[i] = imflat[sm<=v0][-1]
except:
contour_levels[i] = imflat[0]
contour_levels.sort()
m = np.diff(contour_levels) == 0
if np.any(m):
print("Too few points to create valid contours")
while np.any(m):
contour_levels[np.where(m)[0][0]] *= 1.0 - 1e-4
m = np.diff(contour_levels) == 0
contour_levels.sort()
ax.contour(imtmp, contour_levels, **contour_kwargs)
#################################################
ax.plot(center[0], center[1], '+', mew=1., ms=10., c=color_center)
xlim = ax.get_xlim()
ylim = ax.get_ylim()
ybase_offset = 0.035
x_base = xlim[0] + (xlim[1]-xlim[0])*0.075
y_base = ylim[0] + (ylim[1]-ylim[0])*(ybase_offset+0.075)
string = r'$[{:0.1f},{:0.1f}]$'.format(v_slice_lims[0], v_slice_lims[1])
ax.annotate(string, xy=(x_base, y_base),
xycoords="data",
xytext=(0,0),
color='black',
textcoords="offset points", ha="left", va="center",
fontsize=8)
if not show_pix_coords:
ax.set_xticks([])
ax.set_yticks([])
return ax
def plot_3D_data_automask_info(obs, mask_dict, axes=None):
if axes is None:
# Setup plotting
return_axes = False
nrows = 1; ncols = 5
padx = pady = 0.2
xextra = 0.15; yextra = 0.
scale = 2.5
fig = plt.figure()
fig.set_size_inches((ncols+(ncols-1)*padx+xextra)*scale,
(nrows+(nrows-1)*pady+yextra)*scale )
gs = gridspec.GridSpec(nrows, ncols, hspace=pady, wspace=padx)
axes = []
for jj in range(nrows):
for mm in range(ncols):
axes.append(plt.subplot(gs[jj,mm]))
else:
return_axes = True
int_mode = "nearest"; origin = 'lower'; cmap=cmap_viridis
xcenter = obs.mod_options.xcenter
ycenter = obs.mod_options.ycenter
if xcenter is None:
xcenter =(obs.data.data.shape[2]-1)/2.
if ycenter is None:
ycenter =(obs.data.data.shape[1]-1)/2.
titles = ['Collapsed flux', 'Collapsed err', 'Segm', 'Mask2D', 'Masked flux']
for i, im in enumerate([mask_dict['fmap_cube_sn'],
mask_dict['emap_cube_sn'],
mask_dict['segm'], mask_dict['mask2D'],
mask_dict['mask2D']*mask_dict['fmap_cube_sn']]):
ax = axes[i]
imax = ax.imshow(im, cmap=cmap, interpolation=int_mode,
origin=origin, vmin=None, vmax=None)
ax.scatter(xcenter, ycenter, color='magenta', marker='+', s=30)
cax, kw = colorbar.make_axes_gridspec(ax, pad=0.01,
fraction=5./101., aspect=20.)
cbar = plt.colorbar(imax, cax=cax)
cbar.ax.tick_params(labelsize=8)
ax.set_title(titles[i], fontsize=12)
axes[i] = ax
if return_axes:
return axes
else:
return None
#############################################################
def plot_model_1D(gal,
best_dispersion=None,
inst_corr=True,
fileout_base=None):
for obs_name in gal.observations:
obs = gal.observations[obs_name]
if obs.instrument.ndim == 1:
if fileout_base is not None:
fileout_obs= "{}_{}.pdf".format(fileout_base, obs.name)
else:
fileout_obs = None
plot_single_obs_model_1D(obs,
best_dispersion=best_dispersion,
inst_corr=inst_corr,
fileout=fileout_obs)
def plot_model_2D(gal,
fileout_base=None,
symmetric_residuals=True,
max_residual=100.,
inst_corr=True,
show_contours=True,
show_ruler=True,
len_ruler=None,
ruler_unit='arcsec',
apply_mask=True,
ruler_loc='lowerleft',
color_annotate='black',
color_bad='white',
**plot_kwargs):
for obs_name in gal.observations:
obs = gal.observations[obs_name]
if obs.instrument.ndim == 2:
if fileout_base is not None:
fileout_obs= "{}_{}.pdf".format(fileout_base, obs.name)
else:
fileout_obs = None
plot_single_obs_model_2D(obs, gal.model,
dscale=gal.dscale,
fileout=fileout_obs,
symmetric_residuals=symmetric_residuals,
max_residual=max_residual,
inst_corr=inst_corr,
show_contours=show_contours,
show_ruler=show_ruler,
len_ruler=len_ruler,
ruler_unit=ruler_unit,
apply_mask=apply_mask,
ruler_loc=ruler_loc,
color_annotate=color_annotate,
color_bad=color_bad,
**plot_kwargs)
#############################################################
def plot_single_obs_model_1D(obs,
# fitvelocity=True,
# fitdispersion=True,
# fitflux=False,
best_dispersion=None,
inst_corr=True,
fileout=None):
######################################
# Setup data/model comparison: if this isn't the fit dimension
# data/model comparison (eg, fit in 2D, showing 1D comparison)
model_data = copy.deepcopy(obs.model_data)
if inst_corr:
model_data.data['dispersion'] = \
np.sqrt( model_data.data['dispersion']**2 - \
obs.instrument.lsf.dispersion.to(u.km/u.s).value**2 )
######################################
# Setup plot:
f = plt.figure()
scale = 3.5
ncols = 0
nrows = 1
keyxtitle = r'$r$ [arcsec]'
keyyarr = []
keyytitlearr = []
keyyresidtitlearr = []
if obs.fit_options.fit_flux:
ncols += 1
keyyarr.append('flux')
keyytitlearr.append('Flux [arb]')
keyyresidtitlearr.append(r'$\mathrm{Flux_{data} - Flux_{model}}$ [arb]')
if obs.fit_options.fit_velocity:
ncols += 1
keyyarr.append('velocity')
keyytitlearr.append(r'$V$ [km/s]')
keyyresidtitlearr.append(r'$V_{\mathrm{data}} - V_{\mathrm{model}}$ [km/s]')
if obs.fit_options.fit_dispersion:
ncols += 1
keyyarr.append('dispersion')
keyytitlearr.append(r'$\sigma$ [km/s]')
keyyresidtitlearr.append(r'$\sigma_{\mathrm{data}} - \sigma_{\mathrm{model}}$ [km/s]')
f.set_size_inches(1.1*ncols*scale, nrows*scale)
gs = gridspec.GridSpec(nrows, ncols, wspace=0.35, hspace=0.2)
errbar_lw = 0.5
errbar_cap = 1.5
axes = []
k = -1
for j in range(ncols):
# Comparison:
axes.append(plt.subplot(gs[0,j]))
k += 1
axes[k].scatter( model_data.rarr, model_data.data[keyyarr[j]],
c='red', marker='s', s=25, lw=1, label=None)
if keyyarr[j] == 'dispersion':
if best_dispersion:
axes[k].axhline(y=best_dispersion, ls='--', color='red', zorder=-1.)
axes[k].set_xlabel(keyxtitle)
axes[k].set_ylabel(keyytitlearr[j])
axes[k].axhline(y=0, ls='--', color='k', zorder=-10.)
f.suptitle(obs.name)
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.show()
def plot_single_obs_model_2D(obs, model,
dscale=None,
fileout=None,
symmetric_residuals=True,
max_residual=100.,
inst_corr=True,
show_contours=True,
show_ruler=True,
len_ruler=None,
ruler_unit='arcsec',
apply_mask=True,
ruler_loc='lowerleft',
color_annotate='black',
color_bad='white',
**plot_kwargs):
if show_contours:
# Set contour defaults, if not specifed:
for key in _kwargs_contour_defaults.keys():
if key not in plot_kwargs.keys():
plot_kwargs[key] = _kwargs_contour_defaults[key]
######################################
# Setup plot:
# f = plt.figure(figsize=(9.5, 6))
# scale = 3.5
f = plt.figure()
scale = 3.5
nrows = 1
ncols = 0
for cond in [obs.fit_options.fit_flux, obs.fit_options.fit_velocity,
obs.fit_options.fit_dispersion]:
if cond:
ncols += 1
f.set_size_inches(1.1*ncols*scale, nrows*scale)
cntr = 0
if obs.fit_options.fit_flux:
cntr += 1
grid_flux = ImageGrid(f, 100+ncols*10+cntr,
nrows_ncols=(1, 1),
direction="row",
axes_pad=0.5,
label_mode="1",
share_all=True,
cbar_location="right",
cbar_mode="each",
cbar_size="5%",
cbar_pad="1%",
)
if obs.fit_options.fit_velocity:
cntr += 1
grid_vel = ImageGrid(f, 100+ncols*10+cntr,
nrows_ncols=(1, 1),
direction="row",
axes_pad=0.5,
label_mode="1",
share_all=True,
cbar_location="right",
cbar_mode="each",
cbar_size="5%",
cbar_pad="1%",
)
if obs.fit_options.fit_dispersion:
cntr += 1
grid_disp = ImageGrid(f, 100+ncols*10+cntr,
nrows_ncols=(1, 1),
direction="row",
axes_pad=0.5,
label_mode="1",
share_all=True,
cbar_location="right",
cbar_mode="each",
cbar_size="5%",
cbar_pad="1%",
)
keyxarr = ['model']
keyxtitlearr = ['Model']
keyyarr, keyytitlearr, grid_arr = ([] for _ in range(3))
if obs.fit_options.fit_flux:
keyyarr.append('flux')
keyytitlearr.append(r'Flux')
grid_arr.append(grid_flux)
msk = np.isfinite(obs.model_data.data['flux'])
flux_vmin = obs.model_data.data['flux'][msk].min()
flux_vmax = obs.model_data.data['flux'][msk].max()
if obs.fit_options.fit_velocity:
keyyarr.append('velocity')
keyytitlearr.append(r'$V$')
grid_arr.append(grid_vel)
msk = np.isfinite(obs.model_data.data['velocity'])
vel_vmin = obs.model_data.data['velocity'][msk].min()
vel_vmax = obs.model_data.data['velocity'][msk].max()
if np.abs(vel_vmax) > 400.:
vel_vmax = 400.
if np.abs(vel_vmin) > 400.:
vel_vmin = -400.
if obs.fit_options.fit_dispersion:
keyyarr.append('dispersion')
keyytitlearr.append(r'$\sigma$')
grid_arr.append(grid_disp)
msk = np.isfinite(obs.model_data.data['dispersion'])
disp_vmin = obs.model_data.data['dispersion'][msk].min()
disp_vmax = obs.model_data.data['dispersion'][msk].max()
if np.abs(disp_vmax) > 500:
disp_vmax = 500.
if np.abs(disp_vmin) > 500:
disp_vmin = 0.
int_mode = "nearest"
origin = 'lower'
cmap = copy.copy(cmap_spectral_r)
cmap.set_bad(color=color_bad)
for j in range(len(keyyarr)):
msk = np.isfinite(obs.model_data.data[keyyarr[j]])
# Also use mask if defined:
msk[~obs.model_data.mask] = False
grid = grid_arr[j]
for ax, k in zip(grid, keyxarr):
im = copy.deepcopy(obs.model_data.data[keyyarr[j]])
if apply_mask:
im[~msk] = np.NaN
if keyyarr[j] == 'flux':
vmin = flux_vmin
vmax = flux_vmax
elif keyyarr[j] == 'velocity':
vel_shift = model.geometries[obs.name].vel_shift.value
im -= vel_shift
vel_vmin -= vel_shift
vel_vmax -= vel_shift
vmin = vel_vmin
vmax = vel_vmax
elif keyyarr[j] == 'dispersion':
if inst_corr:
im = np.sqrt(im ** 2 - obs.instrument.lsf.dispersion.to(
u.km / u.s).value ** 2)
disp_vmin = max(0, np.sqrt(disp_vmin**2 - obs.instrument.lsf.dispersion.to(u.km / u.s).value ** 2))
disp_vmax = np.sqrt(disp_vmax**2 - obs.instrument.lsf.dispersion.to(u.km / u.s).value ** 2)
vmin = disp_vmin
vmax = disp_vmax
imax = ax.imshow(im, cmap=cmap, interpolation=int_mode,
vmin=vmin, vmax=vmax, origin=origin)
ax.set_ylabel(keyytitlearr[j])
for pos in ['top', 'bottom', 'left', 'right']:
ax.spines[pos].set_visible(False)
ax.set_xticks([])
ax.set_yticks([])
if show_contours:
ax = plot_contours_2D_multitype(im, ax=ax, mapname=keyyarr[j], plottype=k,
vmin=vmin, vmax=vmax, kwargs=plot_kwargs)
ax = plot_major_minor_axes_2D(ax, obs, model, im, obs.model_data.mask)
if show_ruler:
pixscale = obs.instrument.pixscale.value
if (len_ruler is None):
len_arcsec = 1.
elif (len_ruler is not None):
if ruler_unit.lower() == 'arcsec':
len_arcsec = len_ruler
elif ruler_unit.lower() == 'kpc':
if dscale is not None:
len_arcsec = len_ruler * dscale
else:
logger.warning("No 'dscale' provided! Using 1arcsec ruler")
len_arcsec = 1.
ruler_unit = 'arcsec'
ax = plot_ruler_arcsec_2D(ax, pixscale, len_arcsec=len_arcsec,
ruler_loc=ruler_loc, color=color_annotate,
ruler_unit=ruler_unit, dscale=dscale)
cbar = ax.cax.colorbar(imax)
cbar.ax.tick_params(labelsize=8)
f.suptitle(obs.name)
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
#############################################################
def plot_model_comparison_2D(obs1=None, obs2=None,
model1=None, model2=None,
show_models=True,
label_gal1='Gal1',
label_gal2='Gal2',
label_residuals='Residuals: Gal2-Gal1',
symmetric_residuals=True,
max_residual=100.,
fileout=None,
vcrop = False,
vcrop_value = 800.,
inst_corr=True,
show_contours=True,
apply_mask=True,
**kwargs):
# Set contour defaults, if not specifed:
for key in _kwargs_contour_defaults.keys():
if key not in kwargs.keys():
kwargs[key] = _kwargs_contour_defaults[key]
######################################
# Setup plot:
ncols = 3
if show_models:
nrows = 3
else:
nrows = 1
padx = pady = 0.25
xextra = 0.25 #0.15
yextra = 0.25
scale = 2.5
f = plt.figure()
f.set_size_inches((ncols+(ncols-1)*padx+xextra)*scale, (nrows+(nrows-1)*pady+yextra)*scale)
padx = 0.2
pady = 0.1
gs02 = gridspec.GridSpec(nrows, ncols, wspace=padx, hspace=pady)
grid_2D = []
for jj in range(nrows):
for mm in range(ncols):
grid_2D.append(plt.subplot(gs02[jj,mm]))
if inst_corr:
inst_corr_sigma = obs1.instrument.lsf.dispersion.to(u.km/u.s).value
inst_corr_sigma2 = obs2.instrument.lsf.dispersion.to(u.km/u.s).value
# Check values are equivalent:
if inst_corr_sigma != inst_corr_sigma2:
raise ValueError
else:
inst_corr_sigma = 0.
# ----------------------------------------------------------------------
# 2D plotting
pixscale = obs1.instrument.pixscale.value
if show_models:
keyyarr = ['gal1', 'gal2', 'residual']
keyytitlearr = [label_gal1, label_gal2, label_residuals]
else:
keyyarr = ['residual']
keyytitlearr = [label_residuals]
keyxarr = ['flux', 'velocity', 'dispersion']
keyxtitlearr = ['Flux', r'$V$', r'$\sigma$']
int_mode = "nearest"
origin = 'lower'
cmap = copy.copy(cmap_spectral_r)
bad_color = 'white'
color_annotate = 'black'
# bad_color = 'black'
# color_annotate = 'white'
cmap.set_bad(color=bad_color)
cmap_resid = copy.copy(cmap_rdbu_r)
cmap_resid.set_bad(color=bad_color)
cmap_resid.set_over(color='magenta')
cmap_resid.set_under(color='blueviolet')
# -----------------------
if show_models:
vel_vmin = disp_vmin = flux_vmin = 999.
vel_vmax = disp_vmax = flux_vmax = -999.
for obs in [obs1, obs2]:
vel_vmin = np.min([vel_vmin, obs.model_data.data['velocity'][obs.model_data.mask].min()])
vel_vmax = np.max([vel_vmin, obs.model_data.data['velocity'][obs.model_data.mask].max()])
if inst_corr:
im = obs.model_data.data['dispersion'].copy()
im = np.sqrt(im ** 2 - inst_corr_sigma ** 2)
msk = obs.model_data.mask.copy()
msk[~np.isfinite(im)] = False
disp_vmin = np.min([disp_vmin, im[msk].min()])
disp_vmax = np.max([disp_vmax, im[msk].max()])
else:
disp_vmin = np.min([disp_vmin, obs.model_data.data['dispersion'][obs.model_data.mask].min()])
disp_vmax = np.max([disp_vmax, obs.model_data.data['dispersion'][obs.model_data.mask].max()])
flux_vmin = np.min([flux_vmin, obs.model_data.data['flux'][obs.model_data.mask].min()])
flux_vmax = np.max([flux_vmax, obs.model_data.data['flux'][obs.model_data.mask].max()])
# Apply vel shift from model:
vel_shift = model1.geometries[obs1.name].vel_shift.value
vel_vmin -= vel_shift
vel_vmax -= vel_shift
# Check for not too crazy:
if vcrop:
if vel_vmin < -vcrop_value:
vel_vmin = -vcrop_value
if vel_vmax > vcrop_value:
vel_vmax = vcrop_value
if disp_vmin < 0:
disp_vmin = 0
if disp_vmax > vcrop_value:
disp_vmax = vcrop_value
alpha_unmasked = 1.
alpha_masked = 0.5
alpha_bkgd = 1.
alpha_aper = 0.8
for mm in range(len(keyyarr)):
for j in range(len(keyxarr)):
kk = mm*len(keyyarr) + j
k = keyyarr[mm]
ax = grid_2D[kk]
xt = keyxtitlearr[j]
yt = keyytitlearr[mm]
# -----------------------------------
if (k == 'gal1') | (k == 'gal2'):
if (k == 'gal1'):
obs = obs1
model = model1
elif (k == 'gal2'):
obs = obs2
model = model2
if keyxarr[j] == 'velocity':
im = obs.model_data.data['velocity'].copy()
im -= model.geometries[obs.name].vel_shift.value
vmin = vel_vmin
vmax = vel_vmax
elif keyxarr[j] == 'dispersion':
im = obs.model_data.data['dispersion'].copy()
im = np.sqrt(im ** 2 - inst_corr_sigma ** 2)
vmin = disp_vmin
vmax = disp_vmax
elif keyxarr[j] == 'flux':
im = obs.model_data.data['flux'].copy()
vmin = flux_vmin
vmax = flux_vmax
mask = obs.model_data.mask
cmaptmp = cmap
gal = None
elif k == 'residual':
if keyxarr[j] == 'velocity':
im = obs2.model_data.data['velocity'].copy() - obs1.model_data.data['velocity'].copy()
im -= model2.geometries[obs2.name].vel_shift.value - model1.geometries[obs1.name].vel_shift.value
if symmetric_residuals:
vmin = -max_residual
vmax = max_residual
elif keyxarr[j] == 'dispersion':
im_model1 = obs1.model_data.data['dispersion'].copy()
im_model1 = np.sqrt(im_model1 ** 2 - inst_corr_sigma ** 2)
im_model2 = obs2.model_data.data['dispersion'].copy()
im_model2 = np.sqrt(im_model2 ** 2 - inst_corr_sigma ** 2)
im = im_model2 - im_model1
if symmetric_residuals:
vmin = -max_residual
vmax = max_residual
elif keyxarr[j] == 'flux':
im = obs2.model_data.data['flux'].copy() - obs1.model_data.data['flux'].copy()
if symmetric_residuals:
if show_models:
fabsmax = np.max(np.abs([flux_vmin, flux_vmax]))
else:
fabsmax = np.max(np.abs(im[np.isfinite(im)]))
vmin = -fabsmax
vmax = fabsmax
if not symmetric_residuals:
vmin = im[np.isfinite(im)].min()
vmax = im[np.isfinite(im)].max()
mask = obs1.model_data.mask
cmaptmp = cmap_resid
else:
raise ValueError("key not supported.")
if apply_mask:
im[~mask] = np.NaN
imax = ax.imshow(im, cmap=cmaptmp, interpolation=int_mode,
vmin=vmin, vmax=vmax, origin=origin)
if show_contours:
ax = plot_contours_2D_multitype(im, ax=ax, mapname=keyxarr[j], plottype=k,
vmin=vmin, vmax=vmax, kwargs=kwargs)
####################################
# Show a 1arcsec line:
ax = plot_ruler_arcsec_2D(ax, pixscale, len_arcsec=1.,
ruler_loc='lowerright', color=color_annotate)
####################################
ax = plot_major_minor_axes_2D(ax, obs1, model1, im, obs1.model_data.mask)
if j == 0:
ax.set_ylabel(yt)
for pos in ['top', 'bottom', 'left', 'right']:
ax.spines[pos].set_visible(False)
ax.set_xticks([])
ax.set_yticks([])
if mm == 0:
ax.set_title(xt)
#########
cax, kw = colorbar.make_axes_gridspec(ax, pad=0.01,
fraction=4.75/101., aspect=20.)
cbar = plt.colorbar(imax, cax=cax, **kw)
cbar.ax.tick_params(labelsize=8)
################
#############################################################
# Save to file:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
#############################################################
def plot_rotcurve_components(gal,
overwrite=False,
overwrite_curve_files=False,
outpath = None,
plotfile = None,
fname_model_matchdata = None,
fname_model_finer = None,
fname_intrinsic = None,
plot_type='pdf',
**plot_kwargs):
if (plotfile is None) & (outpath is None):
raise ValueError
if fname_intrinsic is None:
fname_intrinsic = '{}{}_vcirc_tot_bary_dm.dat'.format(outpath, gal.name)
# Check if the rot curves are done:
if overwrite_curve_files:
curve_files_exist_int = False
else:
curve_files_exist_int = os.path.isfile(fname_intrinsic)
if not curve_files_exist_int:
# *_vcirc_tot_bary_dm.dat
create_vel_profile_files_intrinsic(gal=gal, outpath=outpath,
fname_intrinsic=fname_intrinsic,
overwrite=overwrite_curve_files)
for obs_name in gal.observations:
obs = gal.observations[obs_name]
if obs.instrument.ndim == 1:
plot_single_obs_rotcurve_components(obs, gal.model,
gal_name=gal.name, overwrite=overwrite,
dscale=gal.dscale, z=gal.z,
overwrite_curve_files=overwrite_curve_files,
outpath = outpath, plotfile = plotfile,
fname_model_matchdata = fname_model_matchdata,
fname_model_finer = fname_model_finer,
fname_intrinsic = fname_intrinsic,
plot_type=plot_type, **plot_kwargs)
return None
def plot_single_obs_rotcurve_components(obs, model,
dscale=None, z=None,
gal_name = None,
overwrite=False, overwrite_curve_files=False,
outpath = None,
plotfile = None,
fname_model_matchdata = None,
fname_model_finer = None,
fname_intrinsic = None,
plot_type='pdf',
**plot_kwargs):
if (plotfile is None) & (outpath is None):
raise ValueError
fbase = '{}{}_{}'.format(outpath, gal_name, obs.name)
if plotfile is None:
plotfile = '{}_rot_components.{}'.format(fbase, plot_type)
if fname_model_matchdata is None:
fname_model_matchdata = '{}_out-1dplots.txt'.format(fbase)
if fname_model_finer is None:
fname_model_finer = '{}_out-1dplots_finer_sampling.txt'.format(fbase)
# check if the file exists:
if overwrite:
file_exists = False
else:
file_exists = os.path.isfile(plotfile)
# Check if the rot curves are done:
if overwrite_curve_files:
# curve_files_exist_int = False
curve_files_exist_obs1d = False
file_exists = False
else:
# curve_files_exist_int = os.path.isfile(fname_intrinsic)
curve_files_exist_obs1d = os.path.isfile(fname_model_finer)
if not curve_files_exist_obs1d:
# *_out-1dplots_finer_sampling.txt, *_out-1dplots.txt
create_vel_profile_files_obs1d(obs=obs, model=model, dscale=dscale,
gal_name=gal_name, outpath=outpath,
fname_finer=fname_model_finer,
fname_model_matchdata=fname_model_matchdata,
overwrite=overwrite_curve_files)
if not file_exists:
# ---------------------------------------------------------------------------
# Read in stuff:
model_obs = read_bestfit_1d_obs_file(filename=fname_model_finer)
model_int = read_model_intrinsic_profile(filename=fname_intrinsic)
deg2rad = np.pi/180.
sini = np.sin(model.geometries[obs.name].inc.value*deg2rad)
vel_asymm_drift_sq = model.kinematic_options.get_asymm_drift_profile(model_int.rarr,
model, tracer=obs.tracer)
vsq = model_int.data['vcirc_tot'] ** 2 - vel_asymm_drift_sq
vsq[vsq<0] = 0.
model_int.data['vrot'] = np.sqrt(vsq)
model_int.data['vrot_sini'] = model_int.data['vrot']*sini
sini_l = np.sin(np.max([model.geometries[obs.name].inc.value - 5., 0.])*deg2rad)
sini_u = np.sin(np.min([model.geometries[obs.name].inc.value + 5., 90.])*deg2rad)
model_int.data['vcirc_tot_linc'] = np.sqrt((model_int.data['vrot_sini']/sini_l)**2 + vel_asymm_drift_sq )
model_int.data['vcirc_tot_uinc'] = np.sqrt((model_int.data['vrot_sini']/sini_u)**2 + vel_asymm_drift_sq )
######################################
# Setup plot:
f = plt.figure()
scale = 3.5
ncols = 2
nrows = 1
wspace = 0.2
hspace = 0.2
f.set_size_inches(1.1*ncols*scale, nrows*scale)
gs = gridspec.GridSpec(nrows, ncols, wspace=wspace, hspace=hspace)
axes = []
for i in range(nrows):
for j in range(ncols):
# Comparison:
axes.append(plt.subplot(gs[i,j]))
keyxtitle = r'Radius [arcsec]'
keyxtitle_alt = r'Radius [kpc]'
keyytitle = r'Velocity [km/s]'
keyytitle_fdm = r'DM Fraction'
errbar_lw = 0.5
errbar_cap = 1.5
lw = 1.5
fontsize_ticks = 9.
fontsize_label = 10.
fontsize_ann = 8.
fontsize_title = 10.
fontsize_leg= 7.5
color_arr = ['mediumblue', 'mediumturquoise', 'orange', 'red', 'blueviolet', 'dimgrey']
# ++++++++++++++++++++++++++++++++++++
ax = axes[0]
xlim = [-0.05, np.max([np.max(np.abs(obs.data.rarr)) + 0.5, 2.0])]
xlim2 = np.array(xlim) / dscale
ylim = [0., np.max(model_int.data['vcirc_tot'])*1.15]
ax.plot( model_obs.rarr, model_obs.data['velocity'],
c='red', lw=lw, zorder=3., label=r'$V_{\mathrm{rot}} \sin(i)$ observed')
ax.axhline(y=model.dispersions[obs.tracer].sigma0.value, ls='--', color='blueviolet',
zorder=-20., label=r'Intrinsic $\sigma_0$')
msk = obs.data.mask
if hasattr(obs.data, 'mask_velocity'):
if obs.data.mask_velocity is not None:
msk = obs.data.mask_velocity
# Masked points
ax.errorbar( np.abs(obs.data.rarr[~msk]), np.abs(obs.data.data['velocity'][~msk]),
xerr=None, yerr = obs.data.error['velocity'][~msk],
marker=None, ls='None', ecolor='lightgrey', zorder=4.,
alpha=0.75, lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None )
ax.scatter( np.abs(obs.data.rarr[~msk]), np.abs(obs.data.data['velocity'][~msk]),
edgecolor='lightgrey', facecolor='whitesmoke', marker='s', s=25, lw=1,
zorder=5., label=None)
# Unmasked points
ax.errorbar( np.abs(obs.data.rarr[msk]), np.abs(obs.data.data['velocity'][msk]),
xerr=None, yerr = obs.data.error['velocity'][msk],
marker=None, ls='None', ecolor='dimgrey', zorder=4.,
alpha=0.75, lw = errbar_lw,capthick= errbar_lw,capsize=errbar_cap,label=None )
ax.scatter( np.abs(obs.data.rarr[msk]), np.abs(obs.data.data['velocity'][msk]),
edgecolor='dimgrey', facecolor='white', marker='s', s=25, lw=1, zorder=5., label='Data')
ax2 = ax.twiny()
ax2.plot(model_int.rarr, model_int.data['vrot_sini'],
c='orange', lw=lw, label=r'$V_{\mathrm{rot}} \sin(i)$ intrinsic', zorder=1.)
#
ax2.plot(model_int.rarr, model_int.data['vrot'],
c='mediumturquoise', lw=lw, label=r'$V_{\mathrm{rot}}$ intrinsic', zorder=0.)
####
ax2.fill_between(model_int.rarr, model_int.data['vcirc_tot_linc'], model_int.data['vcirc_tot_uinc'],
color='mediumblue', alpha=0.1, lw=0, label=None, zorder=-1.)
####
ax2.plot(model_int.rarr, model_int.data['vcirc_tot'],
c='mediumblue', lw=lw, label=r'$V_{\mathrm{c}}$ intrinsic', zorder=2.)
ax.annotate(r'{}: {} $z={:0.1f}$'.format(gal_name, obs.name, z), (0.5, 0.96),
xycoords='axes fraction', ha='center', va='top', fontsize=fontsize_title)
ax.set_xlabel(keyxtitle, fontsize=fontsize_label)
ax.set_ylabel(keyytitle, fontsize=fontsize_label)
ax2.set_xlabel(keyxtitle_alt, fontsize=fontsize_label)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax2.set_xlim(xlim2)
ax.tick_params(labelsize=fontsize_ticks)
ax2.tick_params(labelsize=fontsize_ticks)
ax.xaxis.set_major_locator(MultipleLocator(0.5))
ax.xaxis.set_minor_locator(MultipleLocator(0.1))
ax.yaxis.set_major_locator(MultipleLocator(50.))
ax.yaxis.set_minor_locator(MultipleLocator(10.))
ax2.xaxis.set_major_locator(MultipleLocator(5.))
ax2.xaxis.set_minor_locator(MultipleLocator(1.))
ax.set_zorder(ax2.get_zorder()+1)
ax.patch.set_visible(False)
handles, lbls = ax.get_legend_handles_labels()
handles2, lbls2 = ax2.get_legend_handles_labels()
handles_arr = handles2[::-1]
labels_arr = lbls2[::-1]
handles_arr.extend(handles)
labels_arr.extend(lbls)
frameon = True
borderpad = 0.25
markerscale = 1.
labelspacing= 0.25
handletextpad= 0.05
handlelength = 0.
legend = ax.legend(handles_arr, labels_arr,
labelspacing=labelspacing, borderpad=borderpad,
markerscale=markerscale,
handletextpad=handletextpad,
handlelength=handlelength,
loc='lower right',
frameon=frameon, numpoints=1,
scatterpoints=1,
fontsize=fontsize_leg)
for ind,text in enumerate(legend.get_texts()):
legend.legend_handles[ind].set_visible(False)
#text.set_color(color_arr[ind])
try:
text.set_color(legend.legend_handles[ind]._color)
except:
text.set_color(legend.legend_handles[ind]._original_edgecolor)
# ++++++++++++++++++++++++++++++++++++
ax = axes[1]
color_arr = ['mediumblue', 'limegreen', 'purple']
xlim2 = [xlim2[0], np.max(model_int.rarr)+0.1]
xlim = np.array(xlim2) * dscale
ax2 = ax.twiny()
ax2.plot(model_int.rarr, model_int.data['vcirc_tot'],
c='mediumblue', lw=lw, label=r'$V_{\mathrm{c}}$ intrinsic', zorder=2.)
ax2.plot(model_int.rarr, model_int.data['vcirc_bar'],
c='limegreen', lw=lw, label=r'$V_{\mathrm{bar}}$ intrinsic', zorder=1.)
ax2.plot(model_int.rarr, model_int.data['vcirc_dm'],
c='purple', lw=lw, label=r'$V_{\mathrm{DM}}$ intrinsic', zorder=0.)
###
if 'disk+bulge' in model.components.keys():
ax2.axvline(x=model.components['disk+bulge'].r_eff_disk.value, ls='--', color='dimgrey', zorder=-10.)
ax2.annotate(r'$R_{\mathrm{eff}}$',
(model.components['disk+bulge'].r_eff_disk.value + 0.05*(xlim2[1]-xlim2[0]), 0.025*(ylim[1]-ylim[0])), # 0.05
xycoords='data', ha='left', va='bottom', color='dimgrey', fontsize=fontsize_ann)
###
ax3 = ax2.twinx()
fdm = model_int.data['vcirc_dm']**2/model_int.data['vcirc_tot']**2
ax3.plot(model_int.rarr, fdm, ls='-', lw=1,
color='grey', alpha=0.8)
#
ax2.set_zorder(ax3.get_zorder()+1+ax.get_zorder())
ax2.patch.set_visible(False)
ax.patch.set_visible(False)
ax.set_xlabel(keyxtitle, fontsize=fontsize_label)
ax.set_ylabel(keyytitle, fontsize=fontsize_label)
ax2.set_xlabel(keyxtitle_alt, fontsize=fontsize_label)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax2.set_xlim(xlim2)
ax.tick_params(labelsize=fontsize_ticks)
ax2.tick_params(labelsize=fontsize_ticks)
ax.yaxis.set_major_locator(MultipleLocator(50.))
ax.yaxis.set_minor_locator(MultipleLocator(10.))
ax.xaxis.set_major_locator(MultipleLocator(1.))
ax.xaxis.set_minor_locator(MultipleLocator(0.2))
ax2.xaxis.set_major_locator(MultipleLocator(10.))
ax2.xaxis.set_minor_locator(MultipleLocator(2.))
ax3.set_ylabel(keyytitle_fdm, fontsize=fontsize_label-2, color='grey')
ax3.tick_params(axis='y', direction='in', color='black', labelsize=fontsize_ticks-2, colors='grey')
ax3.set_ylim([0.,1.])
ax.yaxis.set_ticks_position('left')
ax2.yaxis.set_ticks_position('right')
ax3.yaxis.set_major_locator(MultipleLocator(0.2))
ax3.yaxis.set_minor_locator(MultipleLocator(0.05))
handles2, lbls2 = ax2.get_legend_handles_labels()
handles_arr = handles2
labels_arr = lbls2
#
legend = ax2.legend(handles_arr, labels_arr,
labelspacing=labelspacing, borderpad=borderpad,
markerscale=markerscale,
handletextpad=handletextpad,
handlelength=handlelength,
loc='lower right',
frameon=frameon, numpoints=1,
scatterpoints=1,
fontsize=fontsize_leg)
#
for ind,text in enumerate(legend.get_texts()):
legend.legend_handles[ind].set_visible(False)
text.set_color(color_arr[ind])
#############################################################
# Save to file:
plt.savefig(plotfile, bbox_inches='tight', dpi=300)
plt.close()
return None
def make_clean_bayesian_plot_names(bayesianResults, short=False):
names = []
for key in bayesianResults.free_param_names.keys():
for i in range(len(bayesianResults.free_param_names[key])):
param = bayesianResults.free_param_names[key][i]
key_nice = " ".join(key.split("_"))
param_nice = " ".join(param.split("_"))
if short:
names.append(param_nice)
else:
names.append(key_nice+': '+param_nice)
return names
#############################################################
def extract_1D_2D_obs_from_cube(obs, model, slit_width=None, slit_pa=None,
aper_dist=None, inst_corr=True, fill_mask=False):
obs_extract = OrderedDict()
obs_extract['extract_2D'] = extract_2D_from_cube(obs.data.data, obs,
errcube=obs.data.error,
modcube=obs.model_data.data,
inst_corr=inst_corr)
obs_extract['extract_1D'] = extract_1D_from_cube(obs.data.data, obs, model,
errcube=obs.data.error,
modcube=obs.model_data.data,
slit_width=slit_width, slit_pa=slit_pa,
aper_dist=aper_dist, inst_corr=inst_corr,
fill_mask=fill_mask)
# Add geometries for these obs to model:
geom1d = copy.deepcopy(model.geometries[obs.name])
geom1d.name = "{}_1D".format(model.geometries[obs.name].name)
geom1d.obs_name = obs_extract['extract_1D'].name
geom2d = copy.deepcopy(model.geometries[obs.name])
geom2d.name = "{}_2D".format(model.geometries[obs.name].name)
geom2d.obs_name = obs_extract['extract_2D'].name
model.add_component(geom1d)
model.add_component(geom2d)
return obs_extract, model
def aper_centers_arcsec_from_cube(datacube, obs, model, mask=None,
slit_width=None, slit_pa=None,
aper_dist=None, fill_mask=False):
if slit_width is None:
try:
slit_width = obs.instrument.beam.major.to(u.arcsec).value
except:
slit_width = obs.instrument.beam.major_fwhm.to(u.arcsec).value
if slit_pa is None:
slit_pa = model.geometries[obs.name].pa.value
if mask is None:
mask = obs.data.mask.copy()
pixscale = obs.instrument.pixscale.value
rpix = slit_width/pixscale/2.
#############################
if aper_dist is None:
# # EVERY PIXEL
aper_dist_pix = 1. #pixscale #rstep
else:
aper_dist_pix = aper_dist/pixscale
# Aper centers: pick roughly number fitting into size:
nx = datacube.shape[2]
ny = datacube.shape[1]
try:
center_pixel = (obs.mod_options.xcenter + model.geometries[obs.name].xshift.value,
obs.mod_options.ycenter + model.geometries[obs.name].yshift.value)
except:
center_pixel = (int(nx / 2.) + model.geometries[obs.name].xshift,
int(ny / 2.) + model.geometries[obs.name].yshift)
cPA = np.cos(slit_pa * np.pi/180.)
sPA = np.sin(slit_pa * np.pi/180.)
###################
if fill_mask:
diag_step = False
if np.abs(cPA) >= np.abs(sPA):
nmax = ny / np.abs(cPA)
if diag_step & (aper_dist_pix == 1.):
aper_dist_pix *= 1. / np.abs(cPA)
else:
nmax = nx / np.abs(sPA)
if diag_step & (aper_dist_pix == 1.):
aper_dist_pix *= 1. / np.abs(sPA)
rMA_arr = None
else:
# Just use unmasked range:
maskflat = np.sum(mask, axis=0)
maskflat[maskflat>0] = 1
mask2D = np.array(maskflat, dtype=bool)
rstep_A = 0.25
rMA_tmp = 0
rMA_arr = []
# PA is to Blue; rMA_arr is [Blue (neg), Red (pos)]
# but for PA definition blue will be pos step; invert at the end
for fac in [1.,-1.]:
ended_MA = False
while not ended_MA:
rMA_tmp += fac * rstep_A
xtmp = rMA_tmp * -sPA + center_pixel[0]
ytmp = rMA_tmp * cPA + center_pixel[1]
if (xtmp < 0) | (xtmp > mask2D.shape[1]-1) | (ytmp < 0) | (ytmp > mask2D.shape[0]-1):
rMA_arr.append(-1.*(rMA_tmp - fac*rstep_A)) # switch sign: pos / blue for calc becomes neg
rMA_tmp = 0
ended_MA = True
elif not mask2D[int(np.round(ytmp)), int(np.round(xtmp))]:
rMA_arr.append(-1.*rMA_tmp) # switch sign: pos / blue for calc becomes neg
rMA_tmp = 0
ended_MA = True
nmax = None
rMA_arr = np.array(rMA_arr)
if rMA_arr is not None:
aper_centers_pix = np.arange(np.sign(rMA_arr[0])*np.floor(np.abs(rMA_arr[0])),
np.sign(rMA_arr[1])*np.floor(np.abs(rMA_arr[1]))+1., 1.)
else:
nap = int(np.floor(nmax/aper_dist_pix)) # If aper_dist_pix = 1, than nmax apertures.
# Otherwise, fewer and more spaced out
# Make odd
if nap % 2 == 0.:
if fill_mask:
nap -= 1
else:
nap += 1
aper_centers_pix = np.linspace(0.,nap-1, num=nap) - int(nap / 2.)
######################################
aper_centers_arcsec = aper_centers_pix*aper_dist_pix*pixscale
return aper_centers_arcsec
###################################################################
###################################################################
###################################################################
def extract_1D_from_cube(datacube, obs, model, errcube=None, modcube=None, mask=None,
slit_width=None, slit_pa=None, aper_dist=None,
inst_corr=True, fill_mask=False):
# ############################################
# Set up the observation
obs1d = Observation(name="{}_1d_extract".format(obs.name),
tracer=obs.tracer)
for key in ['flux', 'velocity', 'dispersion']:
obs1d.fit_options.__dict__['fit_{}'.format(key)] = True
inst1d = copy.deepcopy(obs.instrument)
inst1d.ndim = 1
if slit_width is None:
try:
slit_width = obs.instrument.beam.major.to(u.arcsec).value
except:
slit_width = obs.instrument.beam.major_fwhm.to(u.arcsec).value
if slit_pa is None:
slit_pa = model.geometries[obs.name].pa.value
if mask is None:
mask = obs.data.mask.copy()
#pixscale = obs.instrument.pixscale.value
inst1d.slit_width = slit_width
inst1d.slit_pa = slit_pa
obs1d.mod_options.xcenter = obs.mod_options.xcenter
obs1d.mod_options.ycenter = obs.mod_options.ycenter
data1d, inst1d = extract_1D_from_cube_general(datacube, obs, model, inst1d,
errcube=errcube, mask=mask,
slit_width=slit_width, slit_pa=slit_pa,
aper_dist=aper_dist,
inst_corr=inst_corr, fill_mask=fill_mask)
obs1d.instrument = inst1d
obs1d.data = data1d
if modcube is not None:
mod1d, _ = extract_1D_from_cube_general(modcube, obs, model, copy.deepcopy(inst1d),
errcube=None, mask=mask,
slit_width=slit_width, slit_pa=slit_pa,
aper_dist=aper_dist,
inst_corr=inst_corr, fill_mask=fill_mask)
obs1d.model_data = mod1d
return obs1d
def extract_1D_from_cube_general(datacube, obs, model, inst1d, errcube=None, mask=None,
slit_width=None, slit_pa=None, aper_dist=None,
inst_corr=True, fill_mask=False):
# ############################################
pixscale = inst1d.pixscale.value
slit_width = inst1d.slit_width
slit_pa = inst1d.slit_pa
rpix = slit_width/pixscale/2.
# Aper centers: pick roughly number fitting into size:
nx = datacube.shape[2]
ny = datacube.shape[1]
try:
center_pixel = (obs.mod_options.xcenter + model.geometries[obs.name].xshift.value,
obs.mod_options.ycenter + model.geometries[obs.name].yshift.value)
except:
center_pixel = (int(nx / 2.) + model.geometries[obs.name].xshift,
int(ny / 2.) + model.geometries[obs.name].yshift)
aper_centers_arcsec = aper_centers_arcsec_from_cube(datacube, obs, model,
mask=mask, slit_width=slit_width, slit_pa=slit_pa,
aper_dist=aper_dist, fill_mask=fill_mask)
#######
vel_arr = datacube.spectral_axis.to(u.km/u.s).value
apertures = CircApertures(rarr=aper_centers_arcsec, slit_PA=slit_pa, rpix=rpix,
nx=nx, ny=ny, center_pixel=center_pixel, pixscale=pixscale,
moment=obs.instrument.moment)
data_scaled = datacube.unmasked_data[:].value
if errcube is not None:
ecube = errcube.unmasked_data[:].value * mask
else:
ecube = None
aper_centers, flux1d, vel1d, disp1d = apertures.extract_1d_kinematics(spec_arr=vel_arr,
cube=data_scaled, mask=mask, err=ecube,
center_pixel = center_pixel, pixscale=pixscale)
if not fill_mask:
# Remove points where the fit was bad
ind = np.isfinite(vel1d) & np.isfinite(disp1d)
data1d = Data1D(r=aper_centers[ind], velocity=vel1d[ind],
vel_disp=disp1d[ind], flux=flux1d[ind],
inst_corr=inst_corr)
apertures_redo = CircApertures(rarr=aper_centers_arcsec[ind], slit_PA=slit_pa, rpix=rpix,
nx=nx, ny=ny, center_pixel=center_pixel, pixscale=pixscale,
moment=obs.instrument.moment)
inst1d.apertures = apertures_redo
else:
data1d = Data1D(r=aper_centers, velocity=vel1d,vel_disp=disp1d, flux=flux1d,
inst_corr=inst_corr)
inst1d.apertures = apertures
data1d.profile1d_type = 'circ_ap_cube'
if fill_mask:
# Unmask any fully masked bits, and fill with the other mask:
mask2d = np.sum(mask, axis=0)
whzero = np.where(mask2d == 0)
maskspec = np.sum(np.sum(mask, axis=2), axis=1)
maskspec[maskspec>0] = 1
mask_filled = np.tile(maskspec.reshape((maskspec.shape[0],1,1)), (1, data_scaled.shape[1], data_scaled.shape[2]))
mask[:, whzero[0], whzero[1]] = mask_filled[:, whzero[0], whzero[1]]
if errcube is not None:
ecube = errcube.unmasked_data[:].value * mask
else:
ecube = None
aper_centers, flux1d, vel1d, disp1d = apertures.extract_1d_kinematics(spec_arr=vel_arr,
cube=data_scaled, mask=mask, err=ecube,
center_pixel = center_pixel, pixscale=pixscale)
data1d.filled_mask_data = Data1D(r=aper_centers, velocity=vel1d,
vel_disp=disp1d, flux=flux1d,
inst_corr=inst_corr)
return data1d, inst1d
def extract_2D_from_cube(cube, obs, errcube=None, modcube=None, inst_corr=True):
# Set up the observation
obs2d = Observation(name="{}_2d_extract".format(obs.name),
tracer=obs.tracer)
for key in ['flux', 'velocity', 'dispersion']:
obs2d.fit_options.__dict__['fit_{}'.format(key)] = True
inst2d = copy.deepcopy(obs.instrument)
inst2d.ndim = 2
# cube must be SpectralCube instance!
data2d = extract_2D_from_cube_general(cube, err=errcube, mask=obs.data.mask,
inst_corr=inst_corr,
directly_correct_LSF=False,
LSF_disp_kms=0,
moment=obs.instrument.moment,
pixscale=obs.instrument.pixscale.value,
gauss_extract_with_c=obs.mod_options.gauss_extract_with_c)
obs2d.instrument = inst2d
obs2d.data = data2d
if modcube is not None:
mod2d = extract_2D_from_cube_general(modcube, err=None, mask=obs.data.mask,
inst_corr=inst_corr,
directly_correct_LSF=False,
LSF_disp_kms=0,
moment=obs.instrument.moment,
pixscale=obs.instrument.pixscale.value,
gauss_extract_with_c=obs.mod_options.gauss_extract_with_c)
obs2d.model_data = mod2d
return obs2d
def extract_2D_from_cube_general(cube, err=None, mask=None,
inst_corr=True, pixscale=None, moment=False,
directly_correct_LSF=False, LSF_disp_kms=0,
gauss_extract_with_c=True):
# cube must be SpectralCube instance!
orig_mask = copy.deepcopy(mask)
# mask = BooleanArrayMask(mask= np.array(mask, dtype=bool),
# wcs=cube.wcs)
if orig_mask is None:
mask_start = np.ones(cube.shape)
else:
mask_start = copy.deepcopy(mask)
mask = BooleanArrayMask(mask= np.array(mask_start, dtype=bool),
wcs=cube.wcs)
datacube = SpectralCube(data=cube.unmasked_data[:].value, mask=mask, wcs=cube.wcs)
if err is not None:
err_cube = SpectralCube(data=err.unmasked_data[:].value, mask=mask, wcs=err.wcs)
if moment:
extrac_type = 'moment'
else:
extrac_type = 'gauss'
if extrac_type == 'moment':
flux = datacube.moment0().to(u.km/u.s).value
vel = datacube.moment1().to(u.km/u.s).value
disp = datacube.linewidth_sigma().to(u.km/u.s).value
elif extrac_type == 'gauss':
mom0 = datacube.moment0().to(u.km/u.s).value
mom1 = datacube.moment1().to(u.km/u.s).value
mom2 = datacube.linewidth_sigma().to(u.km/u.s).value
# Clean up NaNs in moms for initial guesses:
mom0[~np.isfinite(mom0)] = 0.0
mom1[~np.isfinite(mom1)] = 0.0
mom2[~np.isfinite(mom2)] = 20.0
flux = np.zeros(mom0.shape)
vel = np.zeros(mom0.shape)
disp = np.zeros(mom0.shape)
# ++++++++++
my_least_chi_squares_1d_fitter = None
if (gauss_extract_with_c) & (_loaded_LeastChiSquares1D):
if gauss_extract_with_c is not None and \
gauss_extract_with_c is not False:
# we will use the C++ LeastChiSquares1D to run the 1d spectral fitting
# but note that if a spectrum has data all too close to zero, it will fail.
# try to prevent this by excluding too low data
# CLEAN DATA
data_cleaned = datacube.unmasked_data[:,:,:].value
if err is not None:
dataerr = err_cube.unmasked_data[:,:,:].value
dataerr[dataerr==0.] = 99.
dataerr[dataerr==99.] = dataerr.min()
else:
dataerr = None
# data_cleaned = copy.deepcopy(datacube.unmasked_data[:,:,:].value)
# this_fitting_mask = 'auto'
this_fitting_mask = None
if mask_start is not None:
this_fitting_mask = copy.copy(mask_start)
# # Only do verbose if logging level is DEBUG or lower
# if logger.level <= logging.DEBUG:
# this_fitting_verbose = True
# else:
# this_fitting_verbose = False
## Force non-verbose, because multiprocessing pool
## resets logging to logger.level = 0....
this_fitting_verbose = False
# do the least chisquares fitting
my_least_chi_squares_1d_fitter = LeastChiSquares1D(\
x = datacube.spectral_axis.to(u.km/u.s).value,
data = data_cleaned,
dataerr = dataerr,
datamask = this_fitting_mask,
initparams = np.array([mom0 / np.sqrt(2 * np.pi) / np.abs(mom2), mom1, mom2]),
nthread = 4,
verbose = this_fitting_verbose)
if my_least_chi_squares_1d_fitter is not None:
logger.debug('my_least_chi_squares_1d_fitter '+str(datetime.datetime.now())) #<DZLIU><DEBUG>#
my_least_chi_squares_1d_fitter.runFitting()
flux = my_least_chi_squares_1d_fitter.outparams[0,:,:] * np.sqrt(2 * np.pi) * my_least_chi_squares_1d_fitter.outparams[2,:,:]
vel = my_least_chi_squares_1d_fitter.outparams[1,:,:]
disp = my_least_chi_squares_1d_fitter.outparams[2,:,:]
flux[np.isnan(flux)] = 0.0 #<DZLIU><DEBUG># 20210809 fixing this bug
logger.debug('my_least_chi_squares_1d_fitter '+str(datetime.datetime.now())) #<DZLIU><DEBUG>#
if np.sum(np.isfinite(vel)) == 0:
alt_fit = True
else:
alt_fit = False
else:
alt_fit = True
if alt_fit:
# HAS ERROR:
if err is not None:
# print("Doing alt fit: astropy model!")
wgt_cube = np.array(datacube.mask._mask, dtype=np.int64) # use mask as weights to get "weighted" solution
try:
ecube = err_cube.filled_data[:].value
ecube[ecube==99.] = ecube.min()
wgt_cube = 1./(ecube)
except:
pass
for i in range(mom0.shape[0]):
for j in range(mom0.shape[1]):
mod = apy_mod.models.Gaussian1D(amplitude=mom0[i,j] / np.sqrt(2 * np.pi * mom2[i,j]**2),
mean=mom1[i,j],
stddev=np.abs(mom2[i,j]))
mod.amplitude.bounds = (0, None)
mod.stddev.bounds = (0, None)
mod.mean.bounds = (datacube.spectral_axis.to(u.km/u.s).value.min(), datacube.spectral_axis.to(u.km/u.s).value.max())
fitter = apy_mod.fitting.LevMarLSQFitter()
wgts = wgt_cube[:,i,j]
if (np.max(np.abs(wgts)) == 0):
wgts = None
########################
# Masked fit:
spec_arr = datacube.spectral_axis.to(u.km/u.s).value
flux_arr = datacube.filled_data[:,i,j].value
if np.isfinite(flux_arr).sum() >= len(mod._parameters):
if wgts is not None:
whgood = (np.isfinite(flux_arr) & np.isfinite(wgts))
else:
whgood = (np.isfinite(flux_arr))
spec_arr = spec_arr[whgood]
if wgts is not None:
wgts = wgts[whgood]
flux_arr = flux_arr[whgood]
########################
try:
best_fit = fitter(mod, spec_arr, flux_arr, weights=wgts)
flux[i,j] = np.sqrt( 2. * np.pi) * best_fit.stddev.value * best_fit.amplitude
vel[i,j] = best_fit.mean.value
disp[i,j] = best_fit.stddev.value
except:
flux[i, j] = vel[i,j] = disp[i,j] = np.nan
else:
flux[i, j] = vel[i,j] = disp[i,j] = np.nan
# NO ERROR:
else:
for i in range(mom0.shape[0]):
for j in range(mom0.shape[1]):
if i==0 and j==0:
logger.debug('gaus_fit_sp_opt_leastsq '+str(mom0.shape[0])+'x'+str(mom0.shape[1])+' '+str(datetime.datetime.now())) #<DZLIU><DEBUG>#
best_fit = gaus_fit_sp_opt_leastsq(datacube.spectral_axis.to(u.km/u.s).value,
datacube.unmasked_data[:,i,j].value,
mom0[i,j], mom1[i,j], mom2[i,j])
flux[i,j] = best_fit[0] * np.sqrt(2 * np.pi) * best_fit[2]
vel[i,j] = best_fit[1]
disp[i,j] = best_fit[2]
if i==mom0.shape[0]-1 and j==mom0.shape[1]-1:
logger.debug('gaus_fit_sp_opt_leastsq '+str(mom0.shape[0])+'x'+str(mom0.shape[1])+' '+str(datetime.datetime.now())) #<DZLIU><DEBUG>#
# ----------
###########################
# Flatten mask: only mask fully masked spaxels:
msk3d_coll = np.sum(mask_start, axis=0)
whmsk = np.where(msk3d_coll == 0)
mask = np.ones((mask.shape[1], mask.shape[2]))
mask[whmsk] = 0
## Check about instrument correction:
if inst_corr & directly_correct_LSF:
raise ValueError("Should not pre-correct for LSF (`directly_correct_LSF=True`) if setting `inst_corr=True` !!")
if directly_correct_LSF & (LSF_disp_kms > 0.):
disp = np.sqrt(disp**2 - LSF_disp_kms**2)
disp[~np.isfinite(disp)] = 0 # Set the dispersion to zero when its below
# below the instrumental dispersion
# Artificially mask the bad stuff:
mask[~np.isfinite(flux)] = 0
mask[~np.isfinite(vel)] = 0
mask[~np.isfinite(disp)] = 0
flux[~np.isfinite(flux)] = 0
vel[~np.isfinite(vel)] = 0
disp[~np.isfinite(disp)] = 0
if np.sum(mask) == 0:
raise ValueError
data2d = Data2D(pixscale=pixscale, velocity=vel, vel_disp=disp, mask=mask,
flux=flux, vel_err=None, vel_disp_err=None, flux_err=None,
inst_corr=inst_corr)
return data2d
def extract_2D_from_cube_file(fname=None, fname_err=None, fname_mask=None,
moment=False, inst_corr=False,
directly_correct_LSF=False, LSF_disp_kms=0,
xcenter=None, ycenter=None,
gauss_extract_with_c=True):
# LOAD STUFF: load cube into SpectralCube instance!
cube_raw = fits.getdata(fname)
hdr = fits.getheader(fname)
# Should be square:
if np.abs(hdr['CDELT1']) != np.abs(hdr['CDELT2']):
raise ValueError
if hdr['CUNIT1'].strip().upper() in ['DEGREE', 'DEG']:
pixscale = np.abs(hdr['CDELT1']) * 3600. # convert from deg CDELT1 to arcsec
elif hdr['CUNIT1'].strip().upper() in ['ARCSEC']:
pixscale = np.abs(hdr['CDELT1'])
# Define WCS:
w = WCS(hdr)
spec_unit = u.km/u.s
cube = SpectralCube(data=cube_raw, wcs=w).with_spectral_unit(spec_unit)
if fname_err is not None:
err_raw = fits.getdata(fname_err)
err = SpectralCube(data=err_raw, wcs=w).with_spectral_unit(spec_unit)
else:
err = None
if fname_mask is not None:
mask = fits.getdata(fname_mask)
else:
mask = None
data2d = extract_2D_from_cube_general(cube, err=err, mask=mask,
moment=moment, inst_corr=inst_corr,
directly_correct_LSF=directly_correct_LSF,
LSF_disp_kms=LSF_disp_kms, pixscale=pixscale,
gauss_extract_with_c=gauss_extract_with_c)
data2d.xcenter = xcenter
data2d.ycenter = ycenter
return data2d
def plot_axes_flux_vel_disp(flux, vel, disp, axes=None,
mask=None, pixscale=None,
center_pixel_kin=None, PA_deg=None,
plottype='data',
label='',
show_contours=True,
apply_mask=True,
vrange_dict=None,
show_xlabels=True,
cmap=None,
**kwargs):
ims = [flux, vel, disp]
keyxarr = ['flux', 'velocity', 'dispersion']
keyxtitlearr = ['Flux', r'$V$', r'$\sigma$']
int_mode = "nearest"
origin = 'lower'
bad_color = 'white'
color_annotate = 'black'
if cmap is None:
cmap = copy.copy(cmap_spectral_r)
cmap.set_bad(color=bad_color)
for i, im in enumerate(ims):
ax = axes[i]
xt = keyxtitlearr[i]
yt = label
if apply_mask:
im[~mask] = np.NaN
if vrange_dict is not None:
vmin = vrange_dict[keyxarr[i]][0]
vmax = vrange_dict[keyxarr[i]][1]
else:
vmin = None
vmax = None
imax = ax.imshow(im, cmap=cmap, interpolation=int_mode,
vmin=vmin, vmax=vmax, origin=origin)
####################################
if show_contours:
ax = plot_contours_2D_multitype(im, ax=ax, mapname=keyxarr[i],
plottype=plottype,
vmin=vmin, vmax=vmax,
kwargs=kwargs)
# Show a 1arcsec line:
if pixscale is not None:
ax = plot_ruler_arcsec_2D(ax, pixscale, len_arcsec=1.,
ruler_loc='lowerright',
color=color_annotate)
if PA_deg is not None:
ax = plot_major_minor_axes_2D_general(ax, im, mask,
center_pixel_kin=center_pixel_kin, PA_deg=PA_deg)
####################################
if i == 0:
ax.set_ylabel(yt)
for pos in ['top', 'bottom', 'left', 'right']:
ax.spines[pos].set_visible(False)
ax.set_xticks([])
ax.set_yticks([])
if show_xlabels:
ax.set_title(xt)
#########
cax, kw = colorbar.make_axes_gridspec(ax, pad=0.01,
fraction=4.75/101.,
aspect=20.)
cbar = plt.colorbar(imax, cax=cax, **kw)
cbar.ax.tick_params(labelsize=8)
axes[i] = ax
return axes
def plot_2D_from_cube_files(fileout=None,
fname_cube=None, fname_cube2=None,
fname_err=None, fname_err2=None,
fname_mask=None, fname_mask2=None,
show_residual=True,
directly_correct_LSF=False, LSF_disp_kms=0,
directly_correct_LSF2=False, LSF_disp_kms2=0,
label='Gal', label2='Gal2',
moment=False, xcenter=None, ycenter=None,
PA_deg=None,
show_contours=True,
symmetric_residuals=True,
max_residual=100.,
match_ranges=True,
**kwargs):
# Load data
data1 = extract_2D_from_cube_file(fname=fname_cube, fname_err=fname_err,
fname_mask=fname_mask,
inst_corr=False,
directly_correct_LSF=directly_correct_LSF,
LSF_disp_kms=LSF_disp_kms, moment=moment,
xcenter=xcenter, ycenter=ycenter)
if fname_cube2 is not None:
data2 = extract_2D_from_cube_file(fname=fname_cube2, fname_err=fname_err2,
fname_mask=fname_mask2,
inst_corr=False,
directly_correct_LSF=directly_correct_LSF2,
LSF_disp_kms=LSF_disp_kms2, moment=moment,
xcenter=xcenter, ycenter=ycenter)
# Setup plotting
# Set contour defaults, if not specifed:
for key in _kwargs_contour_defaults.keys():
if key not in kwargs.keys():
kwargs[key] = _kwargs_contour_defaults[key]
cmap = copy.copy(cmap_spectral_r)
bad_color = 'white'
cmap.set_bad(color=bad_color)
cmap_resid = copy.copy(cmap_rdbu_r)
cmap_resid.set_bad(color=bad_color)
cmap_resid.set_over(color='magenta')
cmap_resid.set_under(color='blueviolet')
# Setup plot structures:
if xcenter is None:
xcenter = (data1.data['flux'].shape[1]-1.)/ 2.
if ycenter is None:
ycenter = (data1.data['flux'].shape[0]-1.)/ 2.
vrange_dict = {'flux': (None,None),
'velocity': (None,None),
'dispersion': (None,None)}
cube1_dict = {'plottype': 'data',
'cmap': cmap,
'flux': data1.data['flux'],
'vel': data1.data['velocity'],
'disp': data1.data['dispersion'],
'mask': data1.mask,
'center_pixel_kin': (xcenter, ycenter),
'PA_deg': PA_deg,
'pixscale': data1.pixscale,
'label': label,
'vrange_dict': vrange_dict}
cubes_list = [cube1_dict]
if fname_cube2 is not None:
if match_ranges:
mins = []
maxs = []
keys = ['flux', 'velocity', 'dispersion']
for key in keys:
min = np.min([data1.data[key][np.isfinite(data1.data[key])].min(),
data2.data[key][np.isfinite(data2.data[key])].min()])
max = np.max([data1.data[key][np.isfinite(data1.data[key])].max(),
data2.data[key][np.isfinite(data2.data[key])].max()])
mins.append(min)
maxs.append(max)
vrange_dict = {'flux': (mins[0],maxs[0]),
'velocity': (mins[1],maxs[1]),
'dispersion': (mins[2],maxs[2])}
cube1_dict['vrange_dict'] = vrange_dict
cubes_list = [cube1_dict]
cube2_dict = {'plottype': 'data',
'cmap': cmap,
'flux': data2.data['flux'],
'vel': data2.data['velocity'],
'disp': data2.data['dispersion'],
'mask': data2.mask,
'center_pixel_kin': (xcenter, ycenter),
'PA_deg': PA_deg,
'pixscale': data2.pixscale,
'label': label2,
'vrange_dict': vrange_dict}
cubes_list.append(cube2_dict)
if show_residual:
if symmetric_residuals:
vmin = -max_residual
vmax = max_residual
im = cube2_dict['flux']-cube1_dict['flux']
fabsmax = np.max(np.abs(im[np.isfinite(im)]))
fmin = -fabsmax
fmax = fabsmax
else:
vmin = None
vmax = None
fmin = None
fmax = None
vrange_dict_resid = {'flux': (fmin,fmax),
'velocity': (vmin,vmax),
'dispersion': (vmin,vmax)}
resid_dict = {'plottype': 'residual',
'cmap': cmap_resid,
'flux': cube2_dict['flux']-cube1_dict['flux'],
'vel': cube2_dict['vel']-cube1_dict['vel'],
'disp': cube2_dict['disp']-cube1_dict['disp'],
'mask': cube1_dict['mask'],
'center_pixel_kin': (xcenter, ycenter),
'PA_deg': PA_deg,
'pixscale': cube1_dict['pixscale'],
'label': 'Residuals: {}-{}'.format(label2,label),
'vrange_dict': vrange_dict_resid}
cubes_list.append(resid_dict)
######################################
# Setup axes:
ncols = 3
nrows = len(cubes_list)
padx = pady = 0.25
xextra = 0.25 #0.15
yextra = 0.25
scale = 2.5
f = plt.figure()
f.set_size_inches((ncols+(ncols-1)*padx+xextra)*scale,
(nrows+(nrows-1)*pady+yextra)*scale)
padx = 0.2
pady = 0.1
gs02 = gridspec.GridSpec(nrows, ncols, wspace=padx, hspace=pady)
axes = []
for jj in range(nrows):
for mm in range(ncols):
axes.append(plt.subplot(gs02[jj,mm]))
for i, cdict in enumerate(cubes_list):
if i == 0:
show_xlabels = True
else:
show_xlabels = False
axes[ncols*i:(ncols)*(i+1)] = plot_axes_flux_vel_disp(cdict['flux'], cdict['vel'], cdict['disp'],
axes=axes[ncols*i:(ncols)*(i+1)],
mask=cdict['mask'],
pixscale=cdict['pixscale'],
center_pixel_kin=cdict['center_pixel_kin'],
PA_deg=cdict['PA_deg'],
plottype=cdict['plottype'],
label=cdict['label'],
vrange_dict=cdict['vrange_dict'],
cmap=cdict['cmap'],
show_xlabels=show_xlabels,
show_contours=show_contours, **kwargs)
#############################################################
# Save plot to file or display directly:
if fileout is not None:
plt.savefig(fileout, bbox_inches='tight', dpi=300)
plt.close()
else:
plt.draw()
plt.show()
return None
###################################################################
###################################################################
###################################################################
#############################################################
#############################################################
#############################################################
# UTILITY FUNCTIONS
#############################################################
def show_1d_apers_plot(ax, obs, data1d, data2d, model=None, obsorig=None, alpha_aper=0.8, remove_shift=True):
aper_centers = data1d.rarr
slit_width = data1d.slit_width
slit_pa = data1d.slit_pa
rstep = obs.instrument.pixscale.value
try:
rstep1d = obs.instrument1d.pixscale.value
except:
rstep1d = rstep
rpix = slit_width/rstep/2.
aper_centers_pix = aper_centers/rstep
if aper_centers[0] <= 0:
# Starts from neg -> pos:
pa = slit_pa
else:
# Goes "backwards": pos -> neg [but the aper shapes are ALSO stored this way]
pa = slit_pa + 180
print(" ndim={}: xshift={}, yshift={}, vsys2d={}".format(obsorig.data.ndim,
model.geometries[obs.name].xshift.value,
model.geometries[obs.name].yshift.value,
model.geometries[obs.name].vel_shift.value))
nx = data2d.data['velocity'].shape[1]
ny = data2d.data['velocity'].shape[0]
try:
center_pixel_kin = (obs.mod_options.xcenter + model.geometries[obs.name].xshift.value*rstep/rstep1d,
obs.mod_options.ycenter + model.geometries[obs.name].yshift.value*rstep/rstep1d)
except:
center_pixel_kin = (int(nx / 2.) + model.geometries[obs.name].xshift.value*rstep/rstep1d,
int(ny / 2.) + model.geometries[obs.name].yshift.value*rstep/rstep1d)
if not remove_shift:
try:
center_pixel = (obs.mod_options.xcenter, obs.mod_options.ycenter)
except:
center_pixel = None
else:
center_pixel = center_pixel_kin
if center_pixel is None:
center_pixel = (int(nx / 2.) + model.geometries[obs.name].xshift.value*rstep/rstep1d,
int(ny / 2.) + model.geometries[obs.name].yshift.value*rstep/rstep1d)
# +++++++++++++++++
pyoff = 0.
ax.scatter(center_pixel[0], center_pixel[1], color='magenta', marker='+')
ax.scatter(center_pixel_kin[0], center_pixel_kin[1], color='cyan', marker='+')
ax.scatter(int(nx / 2), int(ny / 2), color='lime', marker='+')
# +++++++++++++++++
# First determine the centers of all the apertures that fit within the cube
xaps, yaps = calc_pix_position(aper_centers_pix, pa, center_pixel[0], center_pixel[1])
cmstar = cmap_plasma
cNorm = mplcolors.Normalize(vmin=0, vmax=len(xaps)-1)
cmapscale = mpl.cm.ScalarMappable(norm=cNorm, cmap=cmstar)
for mm, (rap, xap, yap) in enumerate(zip(aper_centers, xaps, yaps)):
circle = plt.Circle((xap+pyoff, yap+pyoff), rpix, color=cmapscale.to_rgba(mm, alpha=alpha_aper), fill=False)
ax.add_artist(circle)
if (mm == 0):
ax.scatter(xap+pyoff, yap+pyoff, color=cmapscale.to_rgba(mm), marker='.')
return ax
def plot_major_minor_axes_2D(ax, obs, model, im, mask, finer_step=True,
lw_major = 3., lw_minor = 2.25, fac2 = 0.66, fac_len_minor_marker = 1./20.,
color_kin_axes = 'black', color_kin_axes2 = 'white'):
####################################
# Show MAJOR AXIS line, center:
try:
center_pixel_kin = (obs.mod_options.xcenter + model.geometries[obs.name].xshift.value,
obs.mod_options.ycenter + model.geometries[obs.name].yshift.value)
except:
center_pixel_kin = ((im.shape[1]-1.)/ 2. + model.geometries[obs.name].xshift.value,
(im.shape[0]-1.)/ 2. + model.geometries[obs.name].yshift.value)
return plot_major_minor_axes_2D_general(ax, im, mask,
center_pixel_kin=center_pixel_kin,
PA_deg=model.geometries[obs.name].pa.value,
finer_step=finer_step,
lw_major=lw_major, lw_minor=lw_minor,
fac2=fac2, fac_len_minor_marker=fac_len_minor_marker,
color_kin_axes=color_kin_axes,
color_kin_axes2=color_kin_axes2)
def plot_major_minor_axes_2D_general(ax, im, mask,
center_pixel_kin=None, PA_deg=None,
finer_step=True,
lw_major = 3., lw_minor = 2.25, fac2 = 0.66, fac_len_minor_marker = 1./20.,
color_kin_axes = 'black', color_kin_axes2 = 'white'):
####################################
# Show MAJOR AXIS line, center:
if center_pixel_kin is None:
center_pixel_kin = ((im.shape[1]-1.)/ 2., (im.shape[0]-1.)/ 2.)
# Start going to neg, pos of center, at PA, and check if mask True/not
# in steps of pix, then rounding. if False: stop, and set 1 less as the end.
cPA = np.cos(PA_deg * np.pi/180.)
sPA = np.sin(PA_deg * np.pi/180.)
A_xs = []
A_ys = []
if finer_step:
rstep_A = 0.25
else:
rstep_A = 1.
rMA_tmp = 0
rMA_arr = []
# PA is to Blue; rMA_arr is [Blue (neg), Red (pos)]
# but for PA definition blue will be pos step; invert at the end
for fac in [1.,-1.]:
ended_MA = False
while not ended_MA:
rMA_tmp += fac * rstep_A
xtmp = rMA_tmp * -sPA + center_pixel_kin[0]
ytmp = rMA_tmp * cPA + center_pixel_kin[1]
if (xtmp < 0) | (xtmp > mask.shape[1]-1) | (ytmp < 0) | (ytmp > mask.shape[0]-1):
A_xs.append((rMA_tmp - fac*rstep_A) * -sPA + center_pixel_kin[0])
A_ys.append((rMA_tmp - fac*rstep_A) * cPA + center_pixel_kin[1])
rMA_arr.append(-1.*(rMA_tmp - fac*rstep_A)) # switch sign: pos / blue for calc becomes neg
rMA_tmp = 0
ended_MA = True
elif not mask[int(np.round(ytmp)), int(np.round(xtmp))]:
A_xs.append((rMA_tmp) * -sPA + center_pixel_kin[0])
A_ys.append((rMA_tmp) * cPA + center_pixel_kin[1])
rMA_arr.append(-1.*rMA_tmp) # switch sign: pos / blue for calc becomes neg
rMA_tmp = 0
ended_MA = True
B_xs = []
B_ys = []
len_minor_marker = fac_len_minor_marker * im.shape[0] # CONSTANT SIZE REL TO AXIS
for fac in [-1.,1.]:
B_xs.append(fac*len_minor_marker*0.5 * cPA + center_pixel_kin[0])
B_ys.append(fac*len_minor_marker*0.5 * sPA + center_pixel_kin[1])
ax.plot(A_xs, A_ys, color=color_kin_axes, lw=lw_major, ls='-')
ax.plot(B_xs, B_ys,color=color_kin_axes, lw=lw_minor, ls='-')
ax.plot(A_xs, A_ys, color=color_kin_axes2, lw=lw_major*fac2, ls='-')
ax.plot(B_xs, B_ys,color=color_kin_axes2, lw=lw_minor*fac2, ls='-')
return ax
def plot_ruler_arcsec_2D(ax, pixscale, len_arcsec=0.5,
ruler_unit='arcsec', dscale=None,
ruler_loc='lowerright', color='black', ybase_offset=0.02,
delx=0.075, dely=0.075, dely_text=0.06):
####################################
# Show a ruler line:
xlim = ax.get_xlim()
ylim = ax.get_ylim()
#len_line_angular = 0.5/(pixscale)
len_line_angular = len_arcsec/(pixscale)
if ruler_unit.lower() == 'arcsec':
if len_arcsec % 1. == 0.:
string = r'{}"'.format(int(len_arcsec))
else:
intpart = str(len_arcsec).split('.')[0]
decpart = str(len_arcsec).split('.')[1]
string = r'{}."{}'.format(intpart, decpart)
elif ruler_unit.lower() == 'kpc':
if len_arcsec/dscale % 1. == 0.:
string = r'{:0.0f} kpc'.format(int(len_arcsec/dscale))
else:
string = r'{:0.1f} kpc'.format(int(len_arcsec/dscale))
#ybase_offset = 0.02 #0.035 #0.065
if 'left' in ruler_loc:
x_base = xlim[0] + (xlim[1]-xlim[0])*delx
sign_x = 1.
ha = 'left'
elif 'right' in ruler_loc:
x_base = xlim[1] - (xlim[1]-xlim[0])*delx
sign_x = -1.
ha = 'right'
if 'upper' in ruler_loc:
y_base = ylim[1] - (ylim[1]-ylim[0])*(ybase_offset+dely)
y_text = y_base - (ylim[1]-ylim[0])*(dely_text)
va = 'center'
elif 'lower' in ruler_loc:
y_base = ylim[0] + (ylim[1]-ylim[0])*(ybase_offset+dely)
y_text = y_base + (ylim[1]-ylim[0])*(dely_text)
va = 'center'
ax.plot([x_base+sign_x*len_line_angular, x_base], [y_base, y_base],
c=color, ls='-',lw=2, solid_capstyle='butt')
ax.annotate(string, xy=(x_base, y_text), xycoords="data", xytext=(0,0),
color=color, textcoords="offset points", ha=ha, va=va, fontsize=9)
return ax
def plot_contours_2D_multitype(im, ax=None, mapname='velocity', plottype='data',
vmin=None, vmax=None, kwargs=None):
# mapname: 'flux', 'velocity', 'dispersion'
# plottype: 'data', 'model', 'gal1', 'gal2', etc -- and then 'residual' handled separately
if mapname == 'velocity':
delta_cont = kwargs['delta_cont_v']
delta_cont_minor = kwargs['delta_cont_v_minor']
delta_cont_minor_resid = kwargs['delta_cont_v_minor_resid']
elif mapname == 'dispersion':
delta_cont = kwargs['delta_cont_disp']
delta_cont_minor = kwargs['delta_cont_disp_minor']
delta_cont_minor_resid = kwargs['delta_cont_disp_minor_resid']
elif mapname == 'flux':
delta_cont = kwargs['delta_cont_flux']
delta_cont_minor = kwargs['delta_cont_flux_minor']
delta_cont_minor_resid = kwargs['delta_cont_flux_minor_resid']
else:
raise ValueError
if vmin is None:
vmin = im[np.isfinite(im)].min()
if vmax is None:
vmax = im[np.isfinite(im)].max()
if (plottype != 'residual'):
#####
if delta_cont_minor is not None:
lo_minor = vmin - (vmin%delta_cont_minor) +delta_cont_minor
hi_minor = vmax -(vmax%delta_cont_minor)
contour_levels_tmp2_minor = np.arange(lo_minor, hi_minor+delta_cont_minor,
delta_cont_minor)
ax.contour(im,levels=contour_levels_tmp2_minor,
colors=kwargs['colors_cont_minor'],
alpha=kwargs['alpha_cont_minor'],
linestyles=kwargs['ls_cont_minor'],
linewidths=kwargs['lw_cont_minor'])
#####
if delta_cont is not None:
lo = vmin - (vmin%delta_cont) +delta_cont
hi = vmax -(vmax%delta_cont)
contour_levels_tmp2 = np.arange(lo, hi+delta_cont, delta_cont)
ax.contour(im,levels=contour_levels_tmp2,
colors=kwargs['colors_cont'], alpha=kwargs['alpha_cont'],
linestyles=kwargs['ls_cont'], linewidths=kwargs['lw_cont'])
elif plottype == 'residual':
# Check that the residual isn't all basically 0:
if delta_cont_minor is not None:
compval = np.min([delta_cont_minor / 10., 0.01])
else:
compval = 0.01
if ((im[np.isfinite(im)].max() - im[np.isfinite(im)].min()) >= compval):
#####
if delta_cont_minor_resid is not None:
# Minor minor contours:
lo_mminor = vmin - (vmin%delta_cont_minor_resid) +delta_cont_minor_resid
hi_mminor = vmax -(vmax%delta_cont_minor_resid)
contour_levels_tmp2_mminor = np.arange(lo_mminor,
hi_mminor+delta_cont_minor_resid,
delta_cont_minor_resid)
ax.contour(im,levels=contour_levels_tmp2_mminor,
colors=kwargs['colors_cont_minor_resid'],
alpha=kwargs['alpha_cont_minor_resid'],
linestyles=kwargs['ls_cont_minor_resid'],
linewidths=kwargs['lw_cont_minor_resid'])
if delta_cont_minor is not None:
lo_minor = vmin - (vmin%delta_cont_minor) +delta_cont_minor
hi_minor = vmax -(vmax%delta_cont_minor)
contour_levels_tmp2_minor = np.arange(lo_minor, hi_minor+delta_cont_minor,
delta_cont_minor)
ax.contour(im,levels=contour_levels_tmp2_minor,
colors=kwargs['colors_cont_minor'],
alpha=kwargs['alpha_cont_minor'],
linestyles=kwargs['ls_cont_minor'],
linewidths=kwargs['lw_cont_minor'])
#####
if delta_cont is not None:
lo = vmin - (vmin%delta_cont) +delta_cont
hi = vmax -(vmax%delta_cont)
contour_levels_tmp2 = np.arange(lo, hi+delta_cont, delta_cont)
ax.contour(im,levels=contour_levels_tmp2,
colors=kwargs['colors_cont'], alpha=kwargs['alpha_cont'],
linestyles=kwargs['ls_cont'], linewidths=kwargs['lw_cont'])
return ax
############
# def _count_nObs_ndim(gal, ndim):
# nObsnD = 0
# for obs_name in gal.observations:
# obs = gal.observations[obs_name]#
# if obs.instrument.ndim == ndim:#
# nObsnD += 1
# return nObsnD