Example Dysmalpy built-in masking routine
This notebook expands on dysmalpy_example_fitting_wrapper_3D.ipynb to demonstrate Dysmalpy’s built-in masking function.
It can mask out skylines and exclude low SNR regions.
Manual inspection of the masking quality is strongly advised.
For more advanced masking, users can provide their own 3D mask, which is beyond the scope of this tutorial.
1) Setup steps
Import modules
import sys,os
from os.path import dirname
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from dysmalpy.fitting_wrappers import dysmalpy_fit_single
from dysmalpy.fitting_wrappers import utils_io
from dysmalpy import fitting, plotting
import numpy as np
import dysmalpy
from astropy.io import fits
from matplotlib.patches import Rectangle
Setup notebook
# Setup plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
mpl.rcParams['figure.dpi']= 300
mpl.rc("savefig", dpi=300)
from IPython.core.display import Image
Set data, output paths
# Check the location of the repository
dysmalpy.__path__
# Data directory (datadir = /YOUR/DATA/PATH/)
filepath = os.path.abspath(fitting.__file__)
datadir = os.sep.join(os.path.dirname(filepath).split(os.sep)[:-1]+["tests", "test_data", ""])
# Load the parameters file from the examples directory
param_path = os.sep.join(os.path.dirname(filepath).split(os.sep)[:-1]+["examples", "examples_param_files", ""])
param_filename = param_path+'fitting_3D_mpfit.params'
# Where to save output files (output = /YOUR/OUTPUTS/PATH). Below the paths of some of the maintainers of the code
outdir = '/Users/lilianlee/Documents/dysmalpy_testing_output_temp/'
outdir = '/Users/jespejo/Dropbox/Postdoc/Data/dysmalpy_test_examples/masking_tutorial/'
outdir_mpfit = outdir + 'MPFIT/'
Settings in parameter file:
Note there are many commented out options / parameters. These given an more complete overview of the settings & parameters that can be specified with the fitting wrapper parameter files.
with open(param_filename, 'r') as f:
print(f.read())
# Example parameters file for fitting a single object with 1D data
# Note: DO NOT CHANGE THE NAMES IN THE 1ST COLUMN AND KEEP THE COMMAS!!
# See README for a description of each parameter and its available options.
# ******************************* OBJECT INFO **********************************
galID, GS4_43501 # Name of your object
z, 1.613 # Redshift
# ****************************** DATA INFO *************************************
datadir, None # Optional: Full path to data directory.
fdata_cube, gs4-43501_h250_21h30.fits.gz # Full path to vel map. Alternatively, just the filename if 'datadir' is set.
fdata_err, noise_gs4-43501_h250_21h30.fits.gz # Full path to vel. err map. Alternatively, just the filename if 'datadir' is set.
fdata_mask, GS4_43501_mask.fits # Full path to mask, if a separate mask is to be used.
# # -- strongly recommended to create a mask, or else use 'auto_gen_3D_mask, True' below
# Alternatively, just the filename if 'datadir' is set.
### FOR 3D, MUST SET CUBE, ETC by ** HAND **.
spec_orig_type, wave
spec_line_rest, 6564.
spec_line_rest_unit, angstrom
spec_vel_trim, -500. 500.
# l r b t
spatial_crop_trim, 37 68 34 65
xcenter, 51.75 #None # x/ycenter in UNCROPPED PIXELS. If None, assumed to be center of cube.
ycenter, 48.75 #None
data_inst_corr, False
auto_gen_3D_mask, True
auto_gen_mask_apply_skymask_first, True
auto_gen_mask_snr_thresh_pixel, None
auto_gen_mask_sig_segmap_thresh, 3.25
auto_gen_mask_npix_segmap_min, 5
auto_gen_mask_sky_var_thresh, 2. # H/J BAND # 3. # K/Y
auto_gen_mask_snr_int_flux_thresh, 3.
# ***************************** OUTPUT *****************************************
outdir, GS4_43501_3D_out/ # Full path for output directory
# ***************************** OBSERVATION SETUP ******************************
# Instrument Setup
# ------------------
pixscale, 0.125 # Pixel scale in arcsec/pixel
fov_npix, 37 # Number of pixels on a side of model cube
spec_type, velocity # DON'T CHANGE!
spec_start, -1000. # Starting value for spectral axis // generally don't change
spec_step, 10. # Step size for spectral axis in km/s // generally don't change
nspec, 201 # Number of spectral steps // generally don't change
# LSF Setup
# ---------
use_lsf, True # True/False if using an LSF
sig_inst_res, 51.0 # Instrumental dispersion in km/s
# PSF Setup
# ---------
psf_type, Gaussian # Gaussian, Moffat, or DoubleGaussian
psf_fwhm, 0.55 # PSF FWHM in arcsecs
psf_beta, -99. # Beta parameter for a Moffat PSF
# ## ELLIPTICAL PSF:
# psf_type, Gaussian # Gaussian, Moffat, or DoubleGaussian
# psf_fwhm_major, 0.55 # PSF major axis FWHM in arcsecs
# psf_fwhm_minor, 0.25 # PSF minor axis FWHM in arcsecs
# psf_PA, 0. # PA of PSF major axis, in deg E of N. (0=N, 90=E)
# psf_beta, -99. # Beta parameter for a Moffat PSF
# # DoubleGaussian: settings instead of psf_fwhm
# psf_type, DoubleGaussian
# psf_fwhm1, 0.16 # FWHM of PSF component 1, in arcsecs. SINFONI AO: 0.16
# psf_fwhm2, 0.48 # FWHM of PSF component 1, in arcsecs. SINFONI AO: 0.48
# psf_scale1, 0.368 # Flux scaling (*not* peak height) of component 1. SINFONI AO: 0.368
# psf_scale2, 0.632 # Flux scaling (*not* peak height) of component 2. SINFONI AO: 0.632
# **************************** SETUP MODEL *************************************
# Model Settings
# -------------
# List of components to use: SEPARATE WITH SPACES
## MUST always keep: geometry zheight_gaus
## RECOMMENDED: always keep: disk+bulge const_disp_prof
components_list, disk+bulge const_disp_prof geometry zheight_gaus halo
# possible options:
# disk+bulge, sersic, blackhole
# const_disp_prof, geometry, zheight_gaus, halo,
# radial_flow, uniform_planar_radial_flow, uniform_bar_flow, uniform_wedge_flow,
# unresolved_outflow, biconical_outflow,
# CAUTION: azimuthal_planar_radial_flow, variable_bar_flow, spiral_flow
# List of components that emit light. SEPARATE WITH SPACES
## Current options: disk+bulge / bulge / disk [corresponding to the mass disk+bulge component],
## also: light_sersic, light_gaussian_ring
light_components_list, disk
# NOTE: if a separate light profile (eg light_sersic) is used,
# this MUST be changed to e.g., 'light_components_list, light_sersic'
adiabatic_contract, False # Apply adiabatic contraction?
pressure_support, True # Apply assymmetric drift correction?
noord_flat, True # Apply Noordermeer flattenning?
oversample, 1 # Spatial oversample factor
oversize, 1 # Oversize factor
zcalc_truncate, True # Truncate in zgal direction when calculating or not
n_wholepix_z_min, 3 # Minimum number of whole pixels in zgal dir, if zcalc_truncate=True
# ********************************************************************************
# DISK + BULGE
# ------------
# Initial Values
total_mass, 11.0 # Total mass of disk and bulge log(Msun)
bt, 0.3 # Bulge-to-Total Ratio
r_eff_disk, 5.0 # Effective radius of disk in kpc
n_disk, 1.0 # Sersic index for disk
invq_disk, 5.0 # disk scale length to zheight ratio for disk
n_bulge, 4.0 # Sersic index for bulge
invq_bulge, 1.0 # disk scale length to zheight ratio for bulge
r_eff_bulge, 1.0 # Effective radius of bulge in kpc
# Fixed? True if its a fixed parameter, False otherwise
total_mass_fixed, False
r_eff_disk_fixed, False
bt_fixed, True
n_disk_fixed, True
r_eff_bulge_fixed, True
n_bulge_fixed, True
# Parameter bounds. Lower and upper bounds
total_mass_bounds, 10.0 13.0
bt_bounds, 0.0 1.0
r_eff_disk_bounds, 0.1 30.0
n_disk_bounds, 1.0 8.0
r_eff_bulge_bounds, 1.0 5.0
n_bulge_bounds, 1.0 8.0
# # ********************************************************************************
# # BLACK HOLE
# # ------------
#
# # Initial Values
# BH_mass, 11. # log(Msun)
#
# # Fixed? True if its a fixed parameter, False otherwise
# BH_mass_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# BH_mass_bounds, 6. 18.
# # ********************************************************************************
# # Separate light profile: (Truncated) Sersic profile
# # ------------
# # Initial values
# L_tot_sersic, 1. # arbitrary units
# lr_eff, 4. # kpc
# lsersic_n, 1. # Sersic index of light profile
# lsersic_rinner, 0. # [kpc] Inner truncation radius of sersic profile. 0 = no truncation
# lsersic_router, inf # [kpc] Outer truncation radius of sersic profile. inf = no truncation
#
# # Fixed? True if its a fixed parameter, False otherwise
# L_tot_sersic_fixed, True
# lr_eff_fixed, False
# lsersic_n_fixed, True
# lsersic_rinner_fixed, True
# lsersic_router_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# L_tot_sersic_bounds, 0. 2. # arbitrary units
# lr_eff_bounds, 0.5 15. # kpc
# lsersic_n_bounds, 0.5 8.
# lsersic_rinner_bounds, 0. 5. # kpc
# lsersic_router_bounds, 4. 20. # kpc
# # ********************************************************************************
# # Separate light profile: Gaussian ring
# # ------------
# # Initial values
# L_tot_gaus_ring, 1. # arbitrary units
# R_peak_gaus_ring, 6. # kpc
# FWHM_gaus_ring, 1. # kpc
#
# # Fixed? True if its a fixed parameter, False otherwise
# L_tot_gaus_ring_fixed, True
# R_peak_gaus_ring_fixed, True
# FWHM_gaus_ring_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# L_tot_gaus_ring_bounds, 0. 2. # arbitrary units
# R_peak_gaus_ring_bounds, 0. 15. # kpc
# FWHM_gaus_ring_bounds, 0.1 10. # kpc
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# DARK MATTER HALO
# ----------------
# Halo type: options: NFW / twopowerhalo / burkert / einasto / dekelzhao
halo_profile_type, NFW
# ** NOTE **: Uncomment the section below corresponding to the selected halo type.
# ********************************************************************************
# NFW halo
# Initial Values
mvirial, 11.5 # Halo virial mass in log(Msun)
halo_conc, 5.0 # Halo concentration parameter
fdm, 0.5 # Dark matter fraction at r_eff_disk
# Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
mvirial_fixed, False
halo_conc_fixed, True
fdm_fixed, False
# Parameter bounds. Lower and upper bounds
mvirial_bounds, 10.0 13.0
halo_conc_bounds, 1.0 20.0
fdm_bounds, 0.0 1.0
# Tie the parameters?
fdm_tied, True # for NFW, fdm_tied=True determines fDM from Mvirial (+baryons)
mvirial_tied, False # for NFW, mvirial_tied=True determines Mvirial from fDM (+baryons)
# ********************************************************************************
# # ********************************************************************************
# # Two-power halo
#
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# alpha, 1. # TPH: inner slope. NFW has alpha=1
# beta, 3. # TPH: outer slope. NFW has beta=3
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# alpha_fixed, False
# beta_fixed, True
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# alpha_bounds, 0.0 3.0
# beta_bounds, 1.0 4.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# alpha_tied, False # for TPH, alpha_tied=True determines alpha from free fDM + other parameters.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Burkert halo
#
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# rB, 10. # Burkert: Halo core radius, in kpc
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# rB_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# rB_bounds, 1.0 20.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# rB_tied, False # for Burkert, rB_tied=True determines rB from free fDM + other parameters.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Einasto halo
# # Initial Values
# mvirial, 11.5 # Halo virial mass in log(Msun)
# halo_conc, 5.0 # Halo concentration parameter
# fdm, 0.5 # Dark matter fraction at r_eff_disk
# alphaEinasto, 1. # Einasto: Halo profile index
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# halo_conc_fixed, True
# fdm_fixed, False
# alphaEinasto_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0
# halo_conc_bounds, 1.0 20.0
# fdm_bounds, 0.0 1.0
# alphaEinasto_bounds, 0.0 2.0
#
# # Tie the parameters?
# fdm_tied, True # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
# mvirial_tied, True # for non-NFW, mvirial_tied=True determines Mvirial from SMHM+fgas + baryon total_mass
# alphaEinasto_tied, False # for Einasto, alphaEinasto_tied=True determines alphaEinasto from free fDM + other params.
#
# ### OTHER SETTINGS:
# mhalo_relation, Moster18 ## SMHM relation to use for tying Mvir to Mbar. options: Moster18 / Behroozi13
#
# fgas, 0.5 # Gas fraction for SMHM inference of Mvir if 'mvirial_tied=True'
# lmstar, -99. # Currently code uses fgas to infer lmstar
# # from fitting baryon total_mass for SMHM relation
# # ********************************************************************************
# # ********************************************************************************
# # Dekel-Zhao halo
# # Initial Values
# mvirial, 12.0 # Halo virial mass in log(Msun)
# s1, 1.5 # Inner logarithmic slope (at resolution r1=0.01*Rvir)
# c2, 25.0 # Concentration parameter (defined relative to c, a)
# fdm, 0.5 # Dark matter fraction at r_eff_disk
#
# # Fixed? True if its a fixed parameter, False otherwise. Also set False if it will be tied (below)
# mvirial_fixed, False
# s1_fixed, False
# c2_fixed, False
# fdm_fixed, False
#
# # Parameter bounds. Lower and upper bounds
# mvirial_bounds, 10.0 13.0 # log(Msun)
# s1_bounds, 0.0 2.0
# c2_bounds, 0.0 40.0
# fdm_bounds, 0.0 1.0
#
# # Tie the parameters?
# mvirial_tied, True # mvirial_tied=True determines Mvirial from fDM, s1, c2.
# s1_tied, True # Tie the s1 to M*/Mvir using best-fit Freundlich+20 (Eqs 45, 47, 48, Table 1)
# c2_tied, True # Tie the c2 to M*/Mvir using best-fit Freundlich+20 (Eqs 47, 49, Table 1)
# fdm_tied, False # for non-NFW, fdm_tied=True determines fDM from other halo params (+baryons)
#
# ### OTHER SETTINGS:
# lmstar, 10.5 # Used to infer s1, c2 if s1_tied or c2_tied = True
#
# # ********************************************************************************
# ********************************************************************************
# INTRINSIC DISPERSION PROFILE
# ------------------
# Initial Values
sigma0, 39.0 # Constant intrinsic dispersion value
# Fixed? True if its a fixed parameter, False otherwise
sigma0_fixed, False
# Parameter bounds. Lower and upper bounds
sigma0_bounds, 5.0 300.0
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# HIGHER ORDER COMPONENTS: INFLOW, OUTFLOW
# ----------------------------------------
# # ********************************************************************************
# # UNIFORM SPHERICAL RADIAL FLOW -- in rhat direction in spherical coordinates
# # radial_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# # ********************************************************************************
# # UNIFORM PLANAR RADIAL FLOW -- in Rhat direction in cylindrical coordinates
# # (eg, radial in galaxy midplane)
# # uniform_planar_radial_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# # ********************************************************************************
# # UNIFORM BAR FLOW -- in xhat direction along bar in cartesian coordinates,
# # with bar at an angle relative to galaxy major axis (blue)
# # uniform_bar_flow
# # -------------------
#
# vbar, -90. # Bar flow [km/s]. Positive: Outflow. Negative: Inflow.
# phi, 90. # Azimuthal angle of bar [degrees], counter-clockwise from blue major axis.
# # Default is 90 (eg, along galaxy minor axis)
# bar_width, 2 # Width of the bar perpendicular to bar direction.
# # Bar velocity only is nonzero between -bar_width/2 < ygal < bar_width/2.
# # ********************************************************************************
# # UNIFORM WEDGE FLOW -- in planar radial flow in cylindrical coordinates, restricted to pos, neg wedges
# # uniform_wedge_flow
# # -------------------
#
# vr, -90. # Radial flow [km/s]. Positive: Outflow. Negative: Inflow.
# theta, 60. # Opening angle of wedge [deg]. (the full angular span)
# phi, 90. # Angle offset relative to the galaxy angle, so the wedge center is at phi.
# # Default: 90 deg, so centered along minor axis
# # ********************************************************************************
# # UNRESOLVED OUTFLOW -- at galaxy center (ie, AGN unresolved outflow)
# # unresolved_outflow
# # -------------------
#
# vcenter, 0. # Central velocity of the Gaussian in km/s
# fwhm, 1000. # FWHM of the Gaussian in km/s
# amplitude, 1.e12 # Amplitude of the Gaussian, for flux in ~M/L=1 luminosity units
# # with the dimming applied ... roughly ....
# # ********************************************************************************
# # BICONICAL OUTFLOW
# # biconical_outflow
# # -------------------
#
# n, 0.5 # Power law index
# vmax, 500. # Maximum velocity of the outflow in km/s
# rturn, 5. # Turn-over radius in kpc of the velocty profile
# thetain, 30. # Half inner opening angle in degrees. Measured from the bicone axis
# dtheta, 20. # Difference between inner and outer opening angle in degrees
# rend, 10. # Maximum radius of the outflow in kpc
# norm_flux, 8. # Log flux amplitude of the outflow at r = 0.
# # Need to check dimming/flux conventions
# tau_flux, 1. # Exponential decay rate of the flux
# biconical_profile_type, both # Type of velocity profile:
# # 'both', 'increase', 'decrease', 'constant'
# biconical_outflow_dispersion, 80. # Dispersion (stddev of gaussian) of biconical outflow, km/s
# # ********************************************************************************
# # VARIABLE BAR FLOW -- in xhat direction along bar in cartesian coordinates,
# # with bar at an angle relative to galaxy major axis (blue)
# # CAUTION!!!
# # variable_bar_flow
# # -------------------
#
# vbar_func_bar_flow, -90.*np.exp(-R/5.) # Bar flow FUNCTION [km/s]. Positive: Outflow. Negative: Inflow.
# phi, 90. # Azimuthal angle of bar [degrees], counter-clockwise from blue major axis.
# # Default is 90 (eg, along galaxy minor axis)
# bar_width, 2 # Width of the bar perpendicular to bar direction.
# # Bar velocity only is nonzero between -bar_width/2 < ygal < bar_width/2.
# # ********************************************************************************
# # AZIMUTHAL PLANAR RADIAL FLOW -- in Rhat direction in cylindrical coordinates
# # (eg, radial in galaxy midplane), with an added azimuthal term
# # CAUTION!!!
# # azimuthal_planar_radial_flow
# # -------------------
#
# vr_func_azimuthal_planar_flow, -90.*np.exp(-R/5.) # Radial flow [km/s].
# # Positive: Outflow. Negative: Inflow.
# m, 2 # Number of modes in the azimuthal pattern. m=0 gives a purely radial profile.
# phi0, 0. # Angle offset relative to the galaxy angle [deg],
# # so the azimuthal variation goes as cos(m(phi_gal - phi0))
# # ********************************************************************************
# # SPIRAL DENSIY WAVE FLOW -- as in Davies et al. 2009, ApJ, 702, 114
# # Here assuming CONSTANT velocity -- try to match real Vrot...
# # CAUTION!!! NO SPACES IN FUNCTION DEFINITONS!
# # spiral_flow
# # -------------------
#
# Vrot_func_spiral_flow, 150.+0.*R # Unperturbed rotation velocity of the galaxy
# dVrot_dR_func_spiral_flow, 0.*R # Derivative of Vrot(R) -- ideally evaluated analytically, otherwise very slow.
# rho0_func_spiral_flow, 1.e11*np.exp(-R/5.) # Unperturbed midplane density profile of the galaxy
# f_func_spiral_flow, (np.sqrt(m**2-2.)*Vrot(R)/cs)*np.log(R) # Function describing the spiral shape, m*phi = f(R)
# # with k = df/dR
# k_func_spiral_flow, (np.sqrt(m**2-2.)*Vrot(R)/cs)/R # Function for radial wavenumber
#
# m, 2 # Number of photometric/density spiral arms.
# cs, 10. # Sound speed of medium, in km/s.
# epsilon, 1. # Density contrast of perturbation (unitless).
# Om_p, 0. # Angular speed of the driving force, Omega_p
# phi0, 0. # Angle offset of the arm winding, in degrees. Default: 0.
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# ********************************************************************************
# ZHEIGHT PROFILE
# ---------------
# Initial Values
sigmaz, 0.9 # Gaussian width of the galaxy in z, in kpc
# Fixed? True if its a fixed parameter, False otherwise
sigmaz_fixed, False
# Parameter bounds. Lower and upper bounds
sigmaz_bounds, 0.1 1.0
# Tie the zheight to the effective radius of the disk?
# If set to True, make sure sigmaz_fixed is False
zheight_tied, True
# GEOMETRY
# --------
# Initial Values
inc, 62. # Inclination of galaxy, 0=face-on, 90=edge-on
pa, 142.
xshift, 0. # pixels
yshift, 0. # pixels
vel_shift, 0. # km/s
# Fixed? True if its a fixed parameter, False otherwise
inc_fixed, False
pa_fixed, False
xshift_fixed, False
yshift_fixed, False
vel_shift_fixed, False
# Parameter bounds. Lower and upper bounds
inc_bounds, 42. 82.
pa_bounds, 132. 152.
xshift_bounds, -1.5 1.5 # pixels
yshift_bounds, -1.5 1.5 # pixels
vel_shift_bounds, -100. 100. # km/s
# **************************** Fitting Settings ********************************
fit_method, mpfit # mcmc, nested, or mpfit
do_plotting, True # Produce all output plots?
# *** Note ***: fitflux and fitdispersion are not supported for 3D,
# because these are implicitly included in the cube!
# MPFIT Settings
#----------------
maxiter, 200 # Maximum number of iterations before mpfit quits
Case 1: Single emission IFU data contaminated by skylines
Visualise the masking regions
The skylines to be masked can be visualised by a position-velocity (PV) diagram extracted along the kinematic major axis at a PA = \(142^{\circ}\).
The PV diagram here is plotted with velocity (wavelength) on the horizontal axis and position on the vertical axis.
The S-shaped emission within the red rectangle that arises from the rotation of the galaxy is our region of interest.
On either side of this S-shaped emission, there are vertical stripes along the velocity axis. These vertical stripes represent the unwanted skylines that we aim to mask out.
We will also mask the regions where the integrated S/N is too low for robust kinematic extraction, which in general corresponds to the outskirt of the galaxies, so the regions on top and below of the “S-shaped” emission will also be masked.
pv_maj = fits.open(f"{datadir}/../test_data_masking/gs4-43501_h250_21h30_pvmaj.fits")
fig, ax = plt.subplots(1,1,tight_layout=True,figsize = (7,1))
ax.imshow(pv_maj[0].data,vmin=-3.383e-18, vmax=3.825e-18, origin="lower")
ha_pv_region = Rectangle((670,0), 25, 90, edgecolor="r", facecolor="none", lw=1)
ax.add_patch(ha_pv_region)
plt.axis('off')
(-0.5, 1374.5, -0.5, 90.5)
Setup masking parameters
Setting up masking parameters is done differently here. Instead of relying on the parameters defined in the parameter file, we now specify the values as arguments when calling utils_io.generate_3D_mask. In our first demo, we will use the default parameters originally specified in the parameter file.
First, you need to download the test fits files from:
https://www.mpe.mpg.de/resources/IR/DYSMALPY/dysmalpy_optional_extra_files/gs4-43501_h250_21h30.fits.gz
and
https://www.mpe.mpg.de/resources/IR/DYSMALPY/dysmalpy_optional_extra_files/noise_gs4-43501_h250_21h30.fits.gz
You can either copy (or move) these files into your datadir directory, or you can create the envirnoment variables in the terminal as:
export FDATA_CUBE=/YOUR/DATA/PATH/gs4-43501_h250_21h30.fits.gz
export FDATA_ERR=/YOUR/DATA/PATH/noise_gs4-43501_h250_21h30.fits.gz
# You can also run the bash commands from within the notebook as (e.g. a directoy of one of the maintainers):
%env FDATA_CUBE=/Users/jespejo/Dropbox/Postdoc/dysmalpy/dysmalpy/dysmalpy/tests/test_data/gs4-43501_h250_21h30.fits.gz
%env FDATA_ERR=/Users/jespejo/Dropbox/Postdoc/dysmalpy/dysmalpy/dysmalpy/tests/test_data/noise_gs4-43501_h250_21h30.fits.gz
env: FDATA_CUBE=/Users/jespejo/Dropbox/Postdoc/dysmalpy/dysmalpy/dysmalpy/tests/test_data/gs4-43501_h250_21h30.fits.gz
env: FDATA_ERR=/Users/jespejo/Dropbox/Postdoc/dysmalpy/dysmalpy/dysmalpy/tests/test_data/noise_gs4-43501_h250_21h30.fits.gz
params = utils_io.read_fitting_params(fname=param_filename)
# check if the files are inside the datadir, if so load them from there:
if os.path.exists(datadir+"gs4-43501_h250_21h30.fits.gz") and os.path.exists(datadir+"noise_gs4-43501_h250_21h30.fits.gz"):
gal = utils_io.load_galaxy(params=params, datadir=datadir, skip_mask=True, skip_automask=True)
# if not, load them from the environment variables
else:
params["fdata_cube"] = os.getenv('FDATA_CUBE')
params["fdata_err"] = os.getenv('FDATA_ERR')
gal = utils_io.load_galaxy(params=params, datadir=None, skip_mask=True, skip_automask=True)
mask_v1, mask_dict_v1 = utils_io.generate_3D_mask(obs=gal.observations['OBS'],
params=params,
sig_segmap_thresh = 3.25,
npix_segmap_min = 5,
snr_int_flux_thresh = 3,
snr_thresh_pixel = None,
sky_var_thresh = 2,
apply_skymask_first = True
)
INFO:DysmalPy:dysmalpy.Galaxy:
********************************************************************
*** INFO ***
instrument.fov[0,1]=(37,37) is being reset to match 3D cube (31, 31)
********************************************************************
INFO:DysmalPy:dysmalpy.Galaxy:
********************************************************************
*** INFO ***
instrument spectral settings are being reset
(spec_type=velocity, spec_start=-1000.00 km / s, spec_step=10.00 km / s, nspec=201)
to match 3D cube
(spec_type=velocity, spec_start=-484.72 km / s, spec_step=34.08 km / s, nspec=29)
********************************************************************
++++++++++++++++++++++++++++++++++++++++++++++++++++++
Creating 3D auto mask with the following settings:
sig_segmap_thresh = 3.25
npix_segmap_min = 5
snr_int_flux_thresh = 3
snr_thresh_pixel = None
sky_var_thresh = 2
apply_skymask_first = True
++++++++++++++++++++++++++++++++++++++++++++++++++++++
Masking segmap: used exclude_percentile=15.0
Examine mask info
Quick look at the integrated spatial flux maps and segmaps
plotting.plot_3D_data_automask_info(gal.observations['OBS'], mask_dict_v1)
Save the mask
utils_io.save_3D_mask(obs=gal.observations['OBS'], mask=mask_v1,
filename=datadir+params['fdata_mask'],
save_uncropped_size=True, overwrite=True)
Mask inspection (optional)
To examine the masked regions, you can multiply the 3D mask with the data cube which is already done here.
By default, the mask does not contain spectral unit information in the header. However, this doesn’t impact its functionality as a mask, as it matches with the cube’s dimensions.
We once again use PV diagram to inspect if the cubes are sufficiently masked.
pv_maj_masked = fits.open(f"{datadir}/../test_data_masking/gs4-43501_h250_21h30_pvmaj_masked.fits")
fig, ax = plt.subplots(1,1,tight_layout=True,figsize = (7,1))
plt.imshow(pv_maj_masked[0].data,vmin=-3.383e-18, vmax=3.825e-18,origin="lower")
ha_pv_region = Rectangle((670,0), 25, 90, edgecolor="r", facecolor="none", lw=1)
ax.add_patch(ha_pv_region)
plt.axis('off')
(-0.5, 1374.5, -0.5, 90.5)
Case 2: Data cube free from skylines contamination
In this case only spatial masking is necessasry. To demonstrate it, we use a mock data cube that simulates a low SNR IFU data. The spectral axis in this case represents velocity.
param_filename = datadir+ "../test_data_masking/lowsn_mock.params"
datadir = os.sep.join(os.path.dirname(filepath).split(os.sep)[:-1]+["tests", "test_data_masking", ""])
params = utils_io.read_fitting_params(fname=param_filename)
gal = utils_io.load_galaxy(params=params,
datadir=datadir,
skip_mask=True,
skip_automask=True,
skip_auto_truncate_crop=True
)
INFO:DysmalPy:dysmalpy.Galaxy:
********************************************************************
*** INFO ***
instrument spectral settings are being reset
(spec_type=velocity, spec_start=-500.00 km / s, spec_step=10.00 km / s, nspec=101)
to match 3D cube
(spec_type=velocity, spec_start=-1000.00 km / s, spec_step=10.00 km / s, nspec=201)
********************************************************************
PV diagram
As previously, the horizontal axis represents the velocity, while the vertical axis represents position.
In contrast to Case 1, this cube is clearly free from skylines contamination, so no spectral masking is necessary.
pv_maj = fits.open(datadir+"../test_data_masking/lowsn_mock_pvmaj.fits")
fig, ax = plt.subplots(1,1,tight_layout=True,figsize = (7,1))
ax.imshow(pv_maj[0].data,vmin=-0.2, vmax=0.45, origin="lower")
mock_pv_region = Rectangle((80,0), 40, 50, edgecolor="r", facecolor="none", lw=1)
ax.add_patch(mock_pv_region)
plt.axis('off')
(-0.5, 200.5, -0.5, 50.5)
mask_v3, mask_dict_v3 = utils_io.generate_3D_mask(obs=gal.observations['OBS'],
params=params,
sig_segmap_thresh = 1.5, #lower than before as the cube is more noisy
npix_segmap_min = 5,
snr_int_flux_thresh = 3,
snr_thresh_pixel = None,
sky_var_thresh = None, #None because this cube doesn't contain skylines
apply_skymask_first = True
)
++++++++++++++++++++++++++++++++++++++++++++++++++++++
Creating 3D auto mask with the following settings:
sig_segmap_thresh = 1.5
npix_segmap_min = 5
snr_int_flux_thresh = 3
snr_thresh_pixel = None
sky_var_thresh = 2
apply_skymask_first = True
++++++++++++++++++++++++++++++++++++++++++++++++++++++
Examine mask info
Quick look at the integrated spatial flux maps and segmaps
plotting.plot_3D_data_automask_info(gal.observations['OBS'], mask_dict_v3)
Save the mask
utils_io.save_3D_mask(obs=gal.observations['OBS'], mask=mask_v3,
filename=datadir+params['fdata_mask'],
save_uncropped_size=True, overwrite=True)
Mask inspection (optional)
pv_maj = fits.open(datadir+"../test_data_masking/lowsn_mock_pvmaj_masked.fits")
fig, ax = plt.subplots(1,1,tight_layout=True,figsize = (7,1))
ax.imshow(pv_maj[0].data,vmin=-0.2, vmax=0.45, origin="lower")
mock_pv_region = Rectangle((80,0), 40, 50, edgecolor="r", facecolor="none", lw=1)
ax.add_patch(mock_pv_region)
plt.axis('off')
(-0.5, 200.5, -0.5, 50.5)
