Models calibration

In the following sections, we will see how to use SpectraPy to obtain the final 2D extracted spectra. As training example we will use a MODS1R MOS frame. In some cases one longslit LUCI frame is also used to highlight the differences between MOS and LS cases.

Optical Model calibration

In order to define the distortions map of the FOV, we must initialize the Optical Model. This step requires: the instrument configuration file and the mask description file. We will start by loading an arc lamp frame for MODS1R

>>> # The MODS instrument configuration file
>>> mods1r = "conf/instruments/mods_G670L.icf"
>>> # The mask description file
>>> mods_mask = "examples/data/mods1r/ID532016.mdf"
>>> mods_arc = "examples/data/mods1r/mods1r.20180121.0073.fits.bz2"

We must initialize the ModelsCalibration class: this is the class used to calibrate all the models.

>>> from spectrapy.modelscalib.calib import ModelsCalibration
>>> calib = ModelsCalibration(mods1r, mask=mods_mask)

The last call opens a DS9 instance used by SpectraPy to display images and regions. All the calibration processes consist in moving and adjusting regions on the frame in order to compute proper models. The idea is that the ModelsCalibration instance shows us the current models solutions plotting regions on the DS9 frame. We can adjust these regions and refit the model using the new regions positions.

First of all, we must create a new Optical Model from scratch, because no OPTModel is yet available for this instrument configuration. In this case, we choose to describe the OPTModel with a polynomial of order 2 in both (\(x\) and \(y\)) directions.

>>> opt = calib.new_opt_model(2, 2)

OPTModel assumes the mask oriented with the \(X\) axis horizontal (from left to right) and \(Y\) axis vertical (from bottom to up). In the next figure we show on the left a dispersed frame, on the right the mask as described by the Optical Model.


As we can see, in this case the mask and the image don’t match, slits are not in the expected position, we can check this discrepancy looking at the tilted slit or at the reference square slit.

Due to instrument particularities and optical reflections, this initial assumption (with the axes oriented like in the picture above) can be not true for all the instruments. To solve this problem, the OPTModel can be flipped on both directions by the methods: flipx and flipy. In this MODS1R case, we flipped the model vertically.

>>> # For optical reasons we must (at least in MODS) flip the model along the Y axis
>>> opt.flipy()

The OPTModel locates the mask slits on the FOV. Since we are working with dispersed images, slits are not visible on the frame. For this reason we will use an arc frame to tune the model: we will use the lambda reference position on the frame as virtual slit.


As lambda reference position we suggest to choose a bright isolated line in the mid region of the dispersion range.

>>> calib.load_image(mods_arc)
>>> calib.plot_opt_model()

The arc image and the expected slits as green boxes will be displayed in the DS9 viewer.


The first time, usually slits will not be in the proper position, we have to tune the model moving 1 these slits on the expected reference line position 2, using the standard DS9 regions commands


Once we have moved all the boxes over the correct lines, we must recompute the model solution, refitting it. The library reads from DS9 the current position of the lists, and adjust the model according with these positions.


In case we don’t want to use some slit and discard it from the fits, we can just remove this region from DS9.

>>> calib.fit_opt_model()

We can visually check the results, plotting again of the recomputed model on the frame in use.

>>> calib.plot_opt_model()

If slits regions, still remain in the proper position, the model is good and we can save it.

>>> opt.writeto("examples/tmp/MODS1R.opt", overwrite=True)

In case the model does not match the proper slit positions, we can try to increase the polynomial order of the model and repeat the previous operations.


In case regions are frozen and you are not able to move them, select the Region option in the DS9 Edit menu.


In case you want follow the exercise, without losing time moving regions, you can use already prepared region file in examples/data/mods1r/regions/opt.reg. You must delete all current regions and replace them with regions contained into the file.

The longslit case

In this section, we will calibrate the optical model for a longlist case. For this example we will use LUCI1 frames acquired with a slit 0.75” width and 60” long, using the low resolution grism. The main difference is that the number of slits in the FOV is just 1.

>>> luci_mask = "conf/masks/luci_LS_0.75.mdf"
>>> luci1 = "conf/instruments/luci_G200LoRes_1.93_1.8.icf"
>>> luci_file = "examples/data/luci1LoRes/luci1.20180202.0181.fits.bz2"

In this case we are using a science frame, since in the instrument configuration file we choose an OH sky line as reference position.

Again we create the OPTModel from scratch and display the slit.

>>> from spectrapy.modelscalib.calib import ModelsCalibration
>>> calib = ModelsCalibration(luci1, mask=luci_mask)
>>> calib.load_image(luci_file)
>>> opt = calib.new_opt_model(1, 1)
>>> calib.plot_opt_model(edit=True)

In this case we have just one slit, and the slit mutal positions can not be used to derive the scaling factor of the Optical Model. You can see how the nominal scale factor is not enough accurate since the slit size (green box on the frame) doesn’t fit the size on the dispersed frame.

So we have to adjust both the slit position and the slit dimensions. This is the reason why we set the parameter edit=True. With this flag ON we can move also the slit corners to fit the real slit position. 3


Be careful: slit is described by a polygon. DON’T ADD corners to this polygon, just move the already existing corners, otherwise the fit will fail!

This resizing allows us to calibrate the OPTModel scale. Once done we can refit the models

>>> calib.fit_opt_model()

And save the fit result

>>> opt.writeto("examples/tmp/LUCI1.opt", overwrite=True)

Already prepared region file is available in examples/data/luci1LoRes/regions/opt.reg

Curvature Model calibration

Once the OPTModel has been calibrated, i.e. the reference lambda is properly located for each list, we can handle the spectra curvatures.

The Curvature Model calibration can be performed either using the previous ModelsCalibration class instance or initializing a new class, loading the saved OPTModel.

>>> #The variables already initialized
>>> mods_mask = "examples/data/mods1r/ID532016.mdf"
>>> mods1r = "conf/instruments/mods_G670L.icf"
>>> mods_arc = "examples/data/mods1r/mods1r.20180121.0073.fits.bz2"
>>> from spectrapy.modelscalib.calib import ModelsCalibration

>>> opt = "examples/tmp/MODS1R.opt"
>>> calib = ModelsCalibration(mods1r, mask=mods_mask, opt=opt)

Like in the OPTModel case, we create a Curvature Model from scratch and display it

>>> crv=calib.new_crv_model(1, 2, 2)

The previous line of code creates a curvature model which is locally described by a straight line, i.e. each spectra trace is described by one 1st order polynomial. The coefficients of these lines change along the FOV, i.e. the curvature of each trace is slightly different spectrum by spectrum. In this example we decide to describe this variation by a 2D polynomial of order 2 by 2.

Since we want to follow the trace of the spectra, the best frame to use for the CRVModel is a through slit flat, which shows clearly the trace edges of the spectra.

>>> mods_flat="examples/data/mods1r/mods1r.20180121.0067.fits.bz2"
>>> calib.load_image(mods_flat)
>>> calib.set_trace_limits(800, 800)

The CRVModel, starting from the reference position defined by the OPTModel, follows the geometry of the spectra both in the blue and the red directions.

The last line of code defines (in pixels) the tracing range we are going to calibrate. In this case we decide to calibrate 800 pixels (both directions) along the dispersion direction around the reference line 4.

>>> calib.plot_crv_model(9)

In this tutorial we decided to use 9 points (marked by green crosses) equally spaced along these 1600 pixels.


DS9 will show the slits position (the blue region) and 9 crosses for each slits along the expected position of the spectra traces. Like we did for the OPTModel, we must adjust in DS9 the crosses along the left side the spectra traces. 5

By default only the left edges of the spectra are used to calibrate the CRVModel. The left and right edges of the spectra are defined by SpectraPy on frames with BU dispersion direction (frames where the dispersion direction goes from bottom to up). According with that: in the frames with UB the dispersion direction the left edges for SpectraPy are the right edges on the frames, if the dispersion direction is LR is the upper side is the left edge and in the RL case the left edge is the lower edge on the frame.



Don’t be afraid to delete some crosses if they follow out of frame or they are in a region where the spectrum signal is too faint, the fitting procedure will use only available crosses

>>> calib.fit_crv_model()
>>> calib.plot_crv_model(20)

Once done, we can refit the model and check it again. During check we can increase the number of points (20 in this case) to better visualize the new solution


If we are satisfied by this solution, we can save it pass on to the IDSModel calibration.

>>> crv.writeto("examples/tmp/MODS1R.crv", overwrite=True)

Otherwise we can increase the degree of the polynomial and repeat the previous steps.


The degree of the local polynomial (1 in this example) can be increased according with the number of crosses used for each slit.


The degrees of the global 2D polynomial is strictly related to the number of slit in the FOV and we MUST take into account of the slit number when we decide its degree.


The first parameter of the set_trace_limits defines the number of pixels in the blue direction, the second parameter the pixels in the red


An already prepared region file is available in examples/data/mods1r/regions/crv.reg

The longslit case

The longlist case shows a single slit very extended along the cross dispersion direction. Due to optical distortions, the spectra produced by this kind of slit could have different sizes in the blue and in the red area. Namely the distance between the left and right edge in blue region slightly differs from the one in the red region.


SpectraPy allow us to address this issue fitting both edges of the spectra.

>>> luci_mask = "conf/masks/luci_LS_0.75.mdf"
>>> luci1 = "conf/instruments/luci_G200LoRes_1.93_1.8.icf"
>>> luci_file = "examples/data/luci1LoRes/luci1.20180202.0181.fits.bz2"
>>> from spectrapy.modelscalib.calib import ModelsCalibration

>>> opt = "examples/tmp/LUCI1.opt"
>>> calib = ModelsCalibration(luci1, mask=luci_mask, opt=opt)
>>> crv=calib.new_crv_model(2, 0, 1)
The last call defines a local curvature model of the 2nd order, described by a global 2D model of:
  • order 0 on X axis: we have just a single slit, it can not change in the FOV moving along X axis

  • order 1 on Y axis: we want a model capable of fitting both edges of the spectrum, i.e. this model can change the spectra curvature along the cross dispersion direction (the Y axis)

>>> calib.set_trace_limits(1000, 1000)
>>> calib.load_image(luci_file)
>>> calib.plot_crv_model(7,  pos=(0, 1))

Unlike the MOS case, we want to plot both edges of the spectra, this is achieved by the pos parameter of the plot_crv_model methos.

SpectraPy describes the slit with a Bezier curve parameterized by a real value t with goes from 0 up to 1. So the left edge of the slit is the slit at t=0 and the right edges is the slit at t=1. The parameter pos=(0,1) in the plot_crv_model call, define which region of the slit we want trace (the edges in this case).


pos parameter can be any number between 0 and 1, that means SpectraPy can show you tracing on every point of the slit. This feature can be useful to work out with very problematic data.

In this case we must adjust both edges of the slit and refit the models 6

>>> calib.fit_crv_model()

And check the results

>>> calib.plot_crv_model(100, pos=(0, 0.325, 1))

And finally, as usual, save the model

>>> crv.writeto("examples/tmp/LUCI1.crv", overwrite=True)

An already prepared region file is available in examples/data/luci1LoRes/regions/crv.reg

Inverse Dispersion Solution calibration

Once the spectra have been located on detector by the Optical Model and geometrically described by the Curvature Model, the wavelength calibration of the 2D spectra can be carried out.

The Inverse Dispersion Solution Model is the model which, following the spectrum curvatures, gives us a relation between wavelength (in Angstrom) and detector pixels.

Like for the CRVModel, the IDSModel can be calibrated either using the previous ModelsCalibration instance or initializing a new class which loads the already saved OPTModel and CRVModel instances.

>>> #The variables already initialized
>>> mods_mask = "examples/data/mods1r/ID532016.mdf"
>>> mods1r = "conf/instruments/mods_G670L.icf"
>>> mods_arc = "examples/data/mods1r/mods1r.20180121.0073.fits.bz2"
>>> from spectrapy.modelscalib.calib import ModelsCalibration

>>> opt = "examples/tmp/MODS1R.opt"
>>> crv = "examples/tmp/MODS1R.crv"
>>> calib = ModelsCalibration(mods1r, mask=mods_mask, opt=opt, crv=crv)

The mathematical description of the IDSModel is similar to the CRVModel: for each slit, one mono dimensional polynomial locates the wavelength positions along the curves described by the CRVModel . The FOV variations of the coefficients of this polynomial are described by a 2D polynomial.

We start creating the IDSModel from scratch

>>> ids = calib.new_ids_model(3, 2, 2)

In this way, we created a model described locally by a 3th order polynomial and globally by 2x2 bi-dimensional polynomial.

In this example for the IDSModel calibration, we use the arc frame

>>> calib.load_image(mods_arc)
>>> NeHg_cat="conf/catalogs/NeHg_hr.dat"

In addition to the arc frame, we need a line catalog to know the expected arc lines position. Catalogs are ASCII files containing line positions, their formats is described the Catalogs paragraph.


Each arcs lamps, can be acquired in different frames separately. During IDS model calibration, it could be useful to add together frames of different arc lamps, in order to span a wider wavelength range. In case we provide a list of files to the load_image method, SpectraPy sums them together and displays the stacked image in DS9.

The global mode

Now we can start with the model calibration. First of all we will show the global approach, by this strategy we will calibrate all the MOS slits at the same time.

>>> calib.plot_ids_model(NeHg_cat)

As we can see in DS9 we will see the expected DS9 lines positions. Since lines catalogs are related to a specific instrument configuration, they can span a wavelength range larger than range covered by the data, this is the reason why we see line region positions of of the frame.

Spectrapy allow us to load just a part of the catalog, using only lines in the relevant range. This is done by the wstart and wend parameter. In this case we ca limit the lines to the MODS red regions

>>> calib.plot_ids_model(NeHg_cat, wstart=5000., wend=10000.)

The method plot_ids_model displays the frame and over it, green arrows on the expected line positions. For each list, we must move the arrowhead in the proper position and then refit the model. 7 In this figure we can see how arrow positions don’t match the line positions. The adjustment is not a simple global shift, we are suppose to adjust the model slit by slit, since the distortions are different along the FOV. Moreover each arrow position on the slit, must be moved of a different amount.


During this adjustment, don’t be afraid to delete lines which can not be properly adjusted, e.g. lines which fall out of frame.

>>> calib.fit_ids_model()

After the fit recomputing, the line positions according with the new model will be plotted.


In case the new computed model well fits the line positions, we can save it.

>>> ids.writeto("examples/tmp/MODS1R.ids",  overwrite=True)

Otherwise we can change the orders of the model and repeat the previous steps.


Already prepared region file is available in examples/data/mods1r/regions/ids_global.reg

Spectra slicing

The first time we are dealing with new data, wavelength calibration can be a very tricky exercise. We have to label emission lines in the raw frame, associating to each of them the proper wavelength value. In literature an on the web, we can find very useful references to complete this task .

For the optical range we can suggest atomic data tables (National Institute of Standards and Technology) or for the near infrared the Rousselot et al. article could be a very useful reference for the OH lines positions.

We can find plots of lamp spectra already calibrated, like this on the LBTO page.


The comparison of this plots with 2D raw frames can be not so straightforward. For this reason SpectraPy allow us to obtain a 1D slice of a single slit using the plot_slice method. We each slit we can cut a slice passing the slit ID to the ModelsCalibration.plot_slice method. In the following example we are plotting a slice for the slit 49.

>>> %pylab
>>> calib.plot_slice('49')

In this case this auxiliary plot shows the 1D slice of the slit 49 allowing us to compare the reference plot, with this slice and adjust properly the regions on the raw frame.

The interactive mode

In some cases, adjusting the entire mask could be very tricky. For example if we are dealing with a very crowded mask, this approach could not be the best choice.

For this reason, beside this global approach, the SpectraPy library provides also an interactive mode. This mode allow us to calibrate the model slit by slit.

The interactive mode starts from the slit closest to center of the FOV (where the distortions should be smaller), so we can focus only on this slit and adjust the solution for it; once done we can switch to the next slit, the slit closest to the current one. Moving to the next slit, SpectraPy will apply to the next slit the solution computed on the previous slit. Since the 2 slits are close each other, the new solution is normally reasonable even for this slit. In this way the manual adjustments to be made are quite small.

We will show this methodology by an example. We starts reading the already computed model and creating the IDS model from scratch, like we did in the global approach.

>>> calib = ModelsCalibration(mods1r, mask=mods_mask, opt=opt, crv=crv)
>>> calib.load_image(mods_arc)
>>> NeHg_cat="conf/catalogs/NeHg_hr.dat"
>>> ids = calib.new_ids_model(3, 2, 2)

and start the iteration process

>>> calib.ids_iter(NeHg_cat, wstart=5000.)

SpectraPy will display in DS9 only the central slit.


We must adjust the solution for this slit 8; once done, we can switch to the next slit, by calling the next method


Calling next SpectraPy performs the following actions:

  • it fits the solution for the current slit already adjust

  • it applies this new solution to the next slit

  • it splits DS9 in 2 frames: in upper frame (or left in case dispersion is bottom-up) there is the slit already adjusted as reference, in the lower frame the new slit to adjust


We must work on this second frame and adjust the solution for this new slit. Once done we can iterate the process calling again the next() method.

When we reach the last slit and adjusted it, we can stop the iteration process and fit the overall model. This is done calling the stop_iter method.

>>> calib.stop_iter()

Calling stop_iter we automatically perform the fit of the IDSModel and restore the visualization with the whole single frame


If we are satisfied of the solution, we can save the model.

>>> ids.writeto('/tmp/MODS1R.ids',  overwrite=True)


In case we need to go back and refine the solution of some slit, the prev method is also available .


Every time we go back and forward with the prev and next method, the current slit solution is refitted by default. In case we may want to browse the slits solutions, without refit their solutions, we can call call prev and next methods with the parameter fit=False.


In case we want stop the iteration without refit the model, we can also call the stop_iter method with the parameter fit=False.


The already prepared region file is available in examples/data/mods1r/regions/ids120.reg. One file for each region is available in the same directory.

The longslit case

Unlike the MOS case, in the longlist case the slit can be quite long, mask flexures and instrument distortions can produce spectrum with curved lines like in frame below. That means the solution of the IDSModel is different in the center of the slit with respect of to the edges.

To address this issue, SpectraPy allows to describe the dispersion solutions, defining the IDSModel solution in several parts of the slits .

>>> luci_mask = "conf/masks/luci_LS_0.75.mdf"
>>> luci1 = "conf/instruments/luci_G200LoRes_1.93_1.8.icf"
>>> luci_file = "examples/data/luci1LoRes/luci1.20180202.0181.fits.bz2"
>>> from spectrapy.modelscalib.calib import ModelsCalibration

>>> opt="examples/tmp/LUCI1.opt"
>>> crv="examples/tmp/LUCI1.crv"
>>> calib = ModelsCalibration(luci1, mask=luci_mask, opt=opt, crv=crv)

We decide to generate a local IDSModel described by 3rd order polynomial, which does not change along the X axis (we have just 1 slit), but which changes along the cross dispersion direction since we want to follow the line curvatures

>>> ids = calib.new_ids_model(3, 0, 2)

Now we can start the calibration procedure, splitting the slit in many pieces (7 in this example)

>>> sky_cat = "conf/catalogs/sky_lr.dat"
>>> calib.load_image(luci_file)
>>> calib.plot_ids_model(sky_cat, wstart=15000., nsplit=7)

The same line position is showed 7 times along the slit. Adjusting the line position for each piece, we will instruct the model to follow the lines curvatures 9


Once done, we refit the model and check the solution on the frame. We can also increase the number of slits for a better check.

>>> calib.fit_ids_model()
>>> calib.plot_ids_model(sky_cat, nsplit=20)

And save the model

>>> ids.writeto("examples/tmp/LUCI1.ids", overwrite=True)

An already prepared region file is available in examples/data/luci1/regions/ids.reg