Computing surface brightness maps¶
If an N-body model snapshot.hdf5 contains all necessary informations (see Setup the N-body model),
the Mockimgs modules provides a script to directly compute surface brightness maps as a .fits file:
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const --rsp_val 1.5 -o output.fits
See Scripts for the full documentation.
It is possible to save several maps obtained from different line of sight using the --nlos option:
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const --rsp_val 1.5 --nlos 4 -o output.fits
Sightlines¶
There are different ways to define one or multiple lines-of-sight.
by default the observer is along the z axis looking at the origin.
one specific line-of-sight is defined with the option:
--los 0 0 1where the three number gives the position of the observer looking at the origin.
using:
--nlos 9 --random_loswill use 9 randomly defined lines-of-sight. A random seed may be defined with the option
--random_seedif only
--nlosis provided, different lines-of-sight are generated by splitting an hemisphere equally. Note that the effective number of line of sightsnloswill be different than the integer provided, as the procedure guarantee that the number of division in phi is twice the one in theta and keep only one polar ligne of sight.if a yaml parameter file is provided with the
-poption, additional parameter may be provided. With the following example:LineOfSights: los : [0,1,0] # line of sight (position of the observer looking at [0,0,0]) grid: n_phi : 3 # number of divisions in azimuthal angle n_theta : 3 # number of divisions in elevation angle d_phi : 10 # azimuthal variation in degree [-d_phi,-d_phi] d_theta : 10 # elevation variation in degree [-d_theta,-d_theta]9 sightlines will be used, centered on [0,1,0] and slightly rotated by 10 degree (covering a 3x3 grid in angle space).
List of instruments¶
The list of available implemented instruments is obtained with the --instruments_list command:
mockimgs_sb_compute_images --instruments_list
Multiple outputs¶
If multiples lines-of-sight are used, a similar number of fits files will
be generated and named lines-of-sight output_LOS0i.fits, with
i from 0 to the number of -1.
Concatenated outputs¶
Instead of saving muliple fits files, it is possible to dump maps in a .pkl file:
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const --rsp_val 1.5 --nlos 4 -o output.pkl
Warning
In .pkl files stored images are not surface brightness images but fluxes.
To convert into a surface brightness image, the following procedure can be used, where images is a list of images contained in the file:
# split into a matrix and parameter
FluxMap,params = images[i]
# apply some filters
FluxMap = libMockimgs.FilterFluxMap(FluxMap,None)
# compute the surface brightness map from the luminosity/flux map
SBmap = libMockimgs.FluxToSurfaceBrightness(FluxMap,params["object_distance"],params["ccd_pixel_area"])
The .pkl file can be displayed with:
mockimgs_sb_show_images output.pkl
The list of pre-defined instruments is given here List of Instruments.
Defining an instrument¶
It is possible to provide mockimgs_sb_compute_images with an instrument which is not pre-defined.
To do so, just define the properties of an instrument in a text file (e.g. instrument.py) with the following content like:
instrument = instrument.Instrument(
name = "myInstrument",
telescope = telescope.Telescope(name="iSIM-170",focal=1500*u.mm),
ccd = ccd.CCD(name="arrakihs_vis",shape=[4096,4096],size=[41*u.mm,41*u.mm]),
filter_type = filters.Filter("GAEA_J"),
)
as explained in The Mockimgs module.
To create the suface brightness map, provide the name of your file to the option --instrument:
mockimgs_sb_compute_images snapshot.hdf5 --instrument instrument.py --distance 35 --rsp_mode const --rsp_val 1.5 -o output.fits
Note that is it convenient to check the properties of the instrument using the --info option:
mockimgs_sb_compute_images --instrument instrument.py --info
Surface brigthness smoothing¶
To get more realistic surface brigthness images, the flux of each particle is spread among several pixels at projection time using a 2D gaussian of fwmh of \(r_h\). There are several parameters we can use to influence this smoothing.
mode option |
additional options |
description |
|---|---|---|
|
|
the softening radius of each particles (stored in |
|
|
the softening radius of each particles is fixed to value provided by |
|
|
the softening radius of each particles is set to \(s_{\max}/(\pi/2) \arctan(r_h/s)\), where \(r_h=\rm{np.rsp}\), \(s_{\max}\) is given by |
|
|
the softening radius of each particles is set to \(\log(r_h/s+1)\), where \(r_h=\rm{np.rsp}\) and \(s\) is given by |
Store the IDs of projected particles¶
In certain cases, it can be useful to store the IDs of particles that fall within a given pixel during projection. This can, for example, help trace back the 3D origin of a structure identified on a 2D surface brightness map.
The mockimgs_sb_compute_images script offers this possibility via the option --IDsMap.
If this option is added and supplemented with a filename as for example:
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --IDsMap IDsMap.npy
the file IDsMap.npy (numpy format) is created and stors the IDs of particles according to the pixel they fall-in.
The following code demonstrates how to retrieve the list of particle IDs that fall within the central
pixel of the image, defined by a pixel radius of 10:
import numpy as np
mat = np.load('IDsMap.npy',allow_pickle=True)
# define a grid
nx = mat.shape[0]
ny = mat.shape[1]
x = np.arange(nx)
y = np.arange(ny)
x,y = np.meshgrid(x,y)
# create a mask
rpix = 10
mask = np.sqrt((x-xc)**2 + (y-yc)**2) < rpix
# ravel the data
x = x.ravel()
y = y.ravel()
mask = mask.ravel()
# keep only pixels in the mask
x = np.compress(mask,x)
y = np.compress(mask,y)
# this will become the list of IDs
IDs = np.array([])
for k in range(len(x)):
ix = int(x[k])
iy = int(y[k])
if mat[ix,iy] is None:
continue
IDs = np.append(IDs,mat[ix,iy])
This second example shows how to create a new matrix that contains the number of particles that fall in each pixel:
import numpy as np
mat = np.load('IDsMap.npy',allow_pickle=True)
matn = np.zeros(mat.shape)
for ix in range(mat.shape[0]):
for iy in range(mat.shape[1]):
if mat[ix,iy] is None:
matn[ix,iy] = 0
else:
matn[ix,iy] = len(mat[ix,iy])
matn = matn.astype(float)
Extracting particles from a snapshot using a 2D mask¶
The list of IDs attached to the pixels of a surface brightness image can be used to extract from a 3D model, any feature detected in 2D in this image.
Lets start with a snapshot snapshot.hdf5 from which we created a surface brightness images
snapshot.fits.gz using mockimgs_sb_compute_images as well as an ID file
snapshot.npy using the --IDsMap option of the mockimgs_sb_compute_images.
Imaging now that we created masks on 2D features identified on the snapshot.fits.gz image.
Those mask serves to tag features. They are saved on a new fits file masks.fits.gz with pixels
belonging to a features having a value that corresponds to the ID of the feature.
To extract the particles from snapshot.hdf5 that fall on a mask, and are responsible to creating
the feature, we can use the following command:
mockimgs_sb_extract_3Dfeatures_from_2Dmask snapshot.npy --masks-file masks.fits.gz -n 1 --nb-file snapshot.hdf5 -o snapshot-LSB1.hdf5
The output file snapshot-LSB1.hdf5 will contains only the extracted particles.
Detail of surface brightness computation¶
Analytical formulae¶
We first define the following quantities:
\(L\) is the luminosity in \(\rm{erg/s}\)
\(d\) is the distance in unit of \(10\,\rm{pc}\)
\(A\) is the pixel area in \(arcsec^2\)
\(M_0\) is the magnitude zero point.
The absolute magnitude \(M\) is related to the luminosity by:
\[M = M_0 -2.5 \log_{10}(L)\]
The apparent magnitude \(m\) is then:
\[m = M + 5 \log_{10}(d)\]
The surface brightness magnitude is defined as
\[SB = m + 2.5 \log_{10}(A)\]
or equivalently
\[SB = M_0 -2.5 \log_{10}(\frac{L}{d^2 A}) = M_0 -2.5 \log_{10}(L) + 5 \log_{10}(d) + 2.5 \log_{10}(A)\]
Implementation¶
The goal is to create a surface brightness map from particles for which we know the magnitude.
The first step is to convert those magnitudes into an extensive quantity akin to the flux at \(10\,\rm{pc}\).
As the flux is proportional to the luminosity, we can use the previous relation, moving the proportionality
constant to an additive constant that we can mergre with the zero point. Moreover, we ignore the
zero point. This is possible a long as we carry the exact inverse transformation once moving from flux to magnitudes.
In instrument.ComputeFLuMap this is coded the following way:
# convert to flux (ignore the zero point)
F = 10**(-M/2.5)
Where M is the magnitude for a given particle.
Once this extensive quantity is defined, we can project particles on a grid and
sum their contribution and get a map akin to a flux at \(10\,\rm{pc}\):
# store in the mass field
nb.mass = F.astype(np.float32)
# create the map
FluxMap = self.ccd.MapUsingGauss(nb.pos,nb.mass,nb.rsp)
The surface brightness is then computed in SurfaceBrightnessMap
(lib.FluxToSurfaceBrightness). The three following steps are performed following the above derivation:
- Transform the extensive quantity akin to a flux (the zero point is ignored) into an absolute magnitude
- \[M = -2.5 \log_{10}(F)\]
- Compute the apparent magnitude, i.e., add the contribution to the distance \(d\) in unit of \(10\,\rm{pc}\)
- \[m = M - 5 \log_{10}(d)\]
- Compute the surface brightness magnitude
- \[SB = m + 2.5 \log_{10}(A)\]
where \(A\) is the pixel area in \(arcsec^2\) The corresponding code is:
# compute the absolute magnitude in each pixel (as before we ignore the zero point)
c = FluxMap==0
FluxMap = np.where(FluxMap==0,1e-40, FluxMap)
MagMap = np.where(c,0,-2.5*np.log10(FluxMap))
# compute the apparent magnitude in each pixel
# Mpc -> 10pc
d = object_distance.to(10*u.pc).value
magMap = MagMap + 5*np.log10(d)
# compute the surface brightness in each pixel
SBMap = np.where(MagMap==0, 100, magMap + 2.5*np.log10(ccd_pixel_area.to(u.arcsec**2).value) )
Scripts¶
mockimgs_sb_compute_images¶
usage: mockimgs_sb_compute_images [-h] [-o OUTPUTFILENAME] [-v] [--unitsfile FILE]
[--distance FLOAT] [--instrument STR] [--magFromField STR]
[--fov FLOAT] [--filter STR] [--rsp_mode STRING]
[--rsp_val FLOAT] [--rsp_max FLOAT] [--rsp_sca FLOAT]
[--rsp_fac FLOAT] [--info] [-p PARAMETER_FILE]
[--remove_identical_coords] [--no-remove_identical_coords]
[--mapping_method MAPPING_METHOD] [--nlos INT]
[--random_los] [--los LOS LOS LOS] [--random_seed IRAND]
[--instruments_list] [--filters_list] [--psf FILE]
[--output-flux] [--IDsMap IDSMAP]
[FILE ...]
compute surface brightness images
positional arguments:
FILE a file name
options:
-h, --help show this help message and exit
-o OUTPUTFILENAME Name of the output file
-v, --verbose verbose mode
--unitsfile FILE pNbody unitsparameters file
--distance FLOAT distance of the object in Mpc
--instrument STR instrument name
--magFromField STR take the magnitude from the pNbody field with name magFromField
--fov FLOAT field of view in arcsec
--filter STR filter name
--rsp_mode STRING smoothing mode
--rsp_val FLOAT smoothing value when set to constant
--rsp_max FLOAT max smoothing value
--rsp_sca FLOAT smoothing scale
--rsp_fac FLOAT rsp factor
--info get info (list of keys) and exit.
-p PARAMETER_FILE Name of the parameter file
--remove_identical_coords
Remove particles with identical 3D coordinates
--no-remove_identical_coords
--mapping_method MAPPING_METHOD
Select kernel for mapping particles to CCD [gauss_old, gauss,
spline]
--nlos INT number of line of sights
--random_los use random line of sights
--los LOS LOS LOS coordinates of a line of sight
--random_seed IRAND random seed
--instruments_list give info on available instruments
--filters_list give info on available filters
--psf FILE a psf fits file
--output-flux store flux assuming
--IDsMap IDSMAP compute and store a matrix that contains IDs of particles
Examples:
--------
mockimgs_sb_compute_images --instruments_list
mockimgs_sb_compute_images --filters_list
mockimgs_sb_compute_images --instrument arrakihs_vis_G --info
mockimgs_sb_compute_images --instrument arrakihs_vis_G --filter SB99_ARK_NIR2 --info
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const -o output.fits
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --nlos 9 -o output.fits
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --nlos 9 --random_los -o output.fits
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --los 1 0 0 -o output.fits
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --fov 60 -o output.fits
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 -p params.yml -o output.fits
mockimgs_sb_compute_images snapshot.hdf5 --instrument arrakihs_vis_G --distance 35 --rsp_mode const --rsp_val 1.5 --nlos 4 -o output.pkl
mockimgs_sb_compute_images snapshot.hdf5 --instrument instrument.py --distance 35 --rsp_mode const --rsp_val 1.5 --nlos 4 -o output.pkl
in this last example, the file instrument.py contains something like :
instrument = instrument.Instrument(
name = "arrakihs_vis_SDSSr",
telescope = telescope.Telescope(name="iSIM-170",focal=1477.5*u.mm),
ccd = ccd.CCD(name="CCD273-84" ,shape=[3400,3400],pixel_size=[12*u.micron,12*u.micron]),
filter_type = filters.Filter("BastI_SDSS_r"),
)
mockimgs_sb_extract_3Dfeatures_from_2Dmask¶
usage: mockimgs_sb_extract_3Dfeatures_from_2Dmask [-h] [-o OUTPUTFILENAME]
[--masks-file MASKS_FILE]
[--nb-file NB_FILE] [-n INT]
FILE
From a 2D mask, extract the corresponding 3D features from the snapshot and save a new snapshot file.
positional arguments:
FILE a numpy file containing IDs of particles, output of
mockimgs_sb_compute_images (with option --IDsMap)
options:
-h, --help show this help message and exit
-o OUTPUTFILENAME Name of the output file containing only extracted particles
--masks-file MASKS_FILE
Name of the fits file containing the masks
--nb-file NB_FILE Name of the nbody file from which particles will be selected
-n INT feature ID
Examples:
--------
mockimgs_sb_extract_3Dfeatures_from_2Dmask snapshot.npy --masks-file masks.fits.gz -n 1 --nb-file snapshot.hdf5 -o snapshot-LSB1.hdf5
mockimgs_sb_show_images¶
usage: mockimgs_sb_show_images [-h] [-o OUTPUTFILENAME] [--sbmin FLOAT] [--sbmax FLOAT]
[--sbcontours [FLOAT ...]] [--ext_kpc FLOAT] [-n INT]
[FILE ...]
display surface brightness images
positional arguments:
FILE a file name
options:
-h, --help show this help message and exit
-o OUTPUTFILENAME Name of the output file
--sbmin FLOAT surface brightness minimum
--sbmax FLOAT surface brightness maximum
--sbcontours [FLOAT ...]
surface brightness contours
--ext_kpc FLOAT field of view extension in kpc
-n INT number of images to display
Examples:
--------
mockimgs_sb_show_images output.pkl
mockimgs_sb_show_images output.pkl -o output.png
mockimgs_sb_show_images output.pkl -n 1 -o output.png
mockimgs_sb_show_images output.pkl --ext_kpc 100 --sbmin 25 --sbcontours 28 30.5 -o output.png
mockimgs_sb_display_fits¶
usage: mockimgs_sb_display_fits [-h] [-o OUTPUTFILENAME] [-n INT] [--ax_unit STRING]
[--ax_max FLOAT] [--sbmin FLOAT] [--sbmax FLOAT]
[--add_axes] [--colorbar] [--colormap COLORMAP]
[--colormap_reverse] [--sbcontours [FLOAT ...]]
[--add-title] [--title STR] [--abs-value]
[FILE ...]
display a set of fits images
positional arguments:
FILE a list of files
options:
-h, --help show this help message and exit
-o OUTPUTFILENAME Name of the output file
-n INT number of images to display
--ax_unit STRING axes units (kpc, arcsec, pixels)
--ax_max FLOAT extention of the image in the axes units
--sbmin FLOAT surface brightness minimum
--sbmax FLOAT surface brightness maximum
--add_axes add axes to the figure
--colorbar add colorbar to the figure
--colormap COLORMAP matplotlib colormap name (e.g. mycmap, tab20c, Greys, jet, binary)
--colormap_reverse reverse colormap
--sbcontours [FLOAT ...]
surface brightness contours
--add-title add a title (automatic)
--title STR title to add
--abs-value plot the absolute value of the flux
Examples:
--------
mockimgs_sb_display_fits --add_axes --ax_unit kpc --ax_max 300 --sbmin 25 --sbmax 31.5 --colorbar --colormap Greys --colormap_reverse *.fits
mockimgs_sb_display_fits *.fits -o output.png
mockimgs_sb_display_fits *.fits
mockimgs_sb_display_fits *.fits --sbcontours 28 30.5 --ax_max 150 --colorbar -o output.png
mockimgs_sb_display_fits *.fits --sbcontours 28 30.5 --ax_max 150 --colorbar --add-title -o output.png
mockimgs_sb_display_fits *.fits --sbcontours 28 30.5 --ax_max 150 --colorbar --colormap mycmap -o output.png
mockimgs_sb_display_fits *.fits --sbcontours 28 30.5 --ax_max 150 --colorbar --colormap mycmap --colormap_reverse -o output.png
mockimgs_sb_display_fits *.fits --sbcontours 28 30.5 --ax_max 150 --colorbar --colormap mycmap --colormap_reverse -o output.png
mockimgs_sb_display_fits *.fits --ax_unit pixels --add_axes
mockimgs_sb_display_fits image.fits --add_axes --ax_unit kpc --ax_max 3
mockimgs_sb_display_fits image.fits --add_axes --ax_unit arcsec --ax_max 10
mockimgs_sb_display_fits image.fits --add_axes --ax_unit arcmin --ax_max 60
mockimgs_sb_display_fits image.fits --sbmin 25 --sbmax 32
mockimgs_sb_display_fits image.fits --sbcontours 28 30.5
mockimgs_sb_display_fits image.fits --colorbar
