Skip to content

Construction of templates

Stellar templates

Handles the creation of the stellar models.

StellarModel

Bases: TemplateFramework

The StellarModel is responsible for the creation of the stellar template for each sub-Instrument.

It allows the user to apply spectral_conditions to select the observations that will be in use for this process.

This object supports the following user parameters:

Note: Also check the User parameters of the parent classes for further customization options of SBART Note: All disk- and user-facing interactions are handled by the parent class

Source code in src/ASTRA/template_creation/StellarModel.py 
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
class StellarModel(TemplateFramework):
    """The StellarModel is responsible for the creation of the stellar template for each sub-Instrument.

    It allows the user to apply [`spectral_conditions`][ASTRA.utils.spectral_conditions] to select the observations that will be in use for this process.

    This object supports the following user parameters:

    **Note:** Also check the **User parameters** of the parent classes for further customization options of SBART
    **Note:** All disk- and user-facing interactions are handled by the parent class

    """

    _object_type = "SpectralModel"
    _name = "Stellar"

    model_type = "Stellar"

    template_map: ClassVar[dict[STELLAR_CREATION_MODE, StellarTemplate]] = {
        STELLAR_CREATION_MODE.Sum: SumStellar,
        STELLAR_CREATION_MODE.OBSERVATION: OBS_Stellar,
        STELLAR_CREATION_MODE.Median: MedianStellar,
        STELLAR_CREATION_MODE.PHOENIX: PHOENIX,
    }

    _default_params = TemplateFramework._default_params + DefaultValues(
        CREATION_MODE=UserParam(
            STELLAR_CREATION_MODE.Sum,
            constraint=ValueFromIterable(STELLAR_CREATION_MODE),
        ),
        ALIGNEMENT_RV_SOURCE=UserParam(
            "DRS", constraint=ValueFromIterable(["DRS", "SBART"])
        ),
        PREVIOUS_SBART_PATH=UserParam("", constraint=ValueFromDtype((str, Path))),
        USE_MERGED_RVS=UserParam(False, constraint=BooleanValue),
    )

    def __init__(
        self, root_folder_path: UI_PATH, user_configs: Optional[UI_DICT] = None
    ):
        """Instantiate the object.

        Parameters
        ----------
        root_folder_path: Union [str, pathlib.Path]
            Path to the folder inside which SBART will store its outputs
        user_configs: Optional[Dict[str, Any]]
            Dictionary with the keys and values of the user parameters that have been described above

        """
        super().__init__(
            mode="", root_folder_path=root_folder_path, user_configs=user_configs
        )

        self._creation_conditions = Empty_condition()
        self.iteration_number = 0
        self.RV_source = None

    def Generate_Model(
        self,
        dataClass: DataClass | None,
        template_configs: Optional[UI_DICT] = None,
        conditions: Optional[ConditionModel] = None,
        force_computation: bool = False,
        store_templates: bool = True,
        previous_sbart_rv=None,
    ) -> None:
        """Apply the spectral conditions to decide which observations to use.

        Then, returns to the model generation as defined in the parent implementation.

        Parameters
        ----------
        dataClass : [type]
            DataClass with the observations
        template_configs : dict
            Dictionary with the parameters needed to control the creation of the template. Each one has its own
            set of requirements
        conditions: None, Condition
            Either None (to not select the OBS that will be used) or :py:mod:`~ASTRAutils.spectral_conditions` to
                restrict the observations in use.
        force_computation: bool
            If True, recompute the stellar templates, even if they exist on disk. By default False
        store_templates: bool
            If True [default], store the templates to disk
        previous_sbart_rv: SBART.data_objects.RV_outputs.RV_holder
            Object that s-BART loaded from disk, containing the RVs from a previous iteration

        Notes
        -----
        * The conditions that are passed to the StellarModel are **only** used for the creation of the stellar template.
         This will **not** reject observations from the RV extraction

        """
        if conditions is not None:
            logger.info("Applying conditions to creation of stellar template")
            self._creation_conditions = conditions

        self.RV_source = self._internal_configs["ALIGNEMENT_RV_SOURCE"]

        if self._internal_configs["ALIGNEMENT_RV_SOURCE"] == "SBART":
            logger.info(
                "{} using {} RVs as the source for stellar template creation",
                self.name,
                self._internal_configs["ALIGNEMENT_RV_SOURCE"],
            )
            if previous_sbart_rv is None:
                msg = "No previous SBART RVs were provided!"
                logger.critical(msg)
                raise custom_exceptions.InvalidConfiguration(msg)

            self.iteration_number = previous_sbart_rv.iteration_number + 1
            self.RV_source = "SBART"

            if dataClass is not None:
                try:
                    dataClass.load_previous_SBART_results(
                        previous_sbart_rv,
                        use_merged_cube=self._internal_configs["USE_MERGED_RVS"],
                    )
                except Exception as e:
                    logger.opt(exception=True).critical("Failed")
                    raise e
            else:
                logger.warning(
                    "No DataClass provided. The RVs will not be loaded into the observations"
                )

        else:
            logger.info(
                "Using CCF RVs as the basis for the creation of the stellar models"
            )

        self.add_relative_path("Stellar", f"Stellar/Iteration_{self.iteration_number}")

        super().Generate_Model(
            dataClass=dataClass,
            template_configs=template_configs,
            attempt_to_load=not force_computation,
            store_templates=False,
        )

        for subInst, temp in self.templates.items():
            temp.update_RV_source_info(
                iteration_number=self.iteration_number,
                RV_source=self.RV_source,
                merged_source=self._internal_configs["USE_MERGED_RVS"],
            )

        if store_templates:
            self.store_templates_to_disk()

    def initialize_modelling_interfaces(self, new_properties):
        for subInst, temp in self.templates.items():
            temp.initialize_modelling_interface()

    def get_interpol_modes(self) -> set[str]:
        return set(temp.interpol_mode for temp in self.templates.values())

    def _compute_template(
        self, data: DataClass, subInstrument: str, user_configs: dict
    ) -> None:
        chosen_template = self.template_map[self._internal_configs["CREATION_MODE"]]
        key = "ALIGNEMENT_RV_SOURCE"
        if key in user_configs:
            logger.warning(
                f"Key <{key}> from Stellar Model over-riding the one from the template configs"
            )
        user_configs[key] = self._internal_configs[key]
        stellar_template = chosen_template(
            subInst=subInstrument, user_configs=user_configs
        )

        try:
            stellar_template.create_stellar_template(
                dataClass=data, conditions=self._creation_conditions
            )
        except NoDataError:
            logger.info(
                "{} has no available data. The template will be created as an array of zeros",
                subInstrument,
            )

        return stellar_template

    def get_orders_to_skip(self, subInst: str) -> set[int]:
        """Compute orders to skip on a given subInstrument.

        GIven by the orders that the stellar template rejected.

        Args:
            subInst (str): sub-Instrument name

        Raises:
            BadTemplateError: No valid template for the subInstrument

        Returns:
            Set[int]: Set of bad orders.

        """
        if subInst == "all":
            bad_orders: set[int] = set()
            for temp in self.templates.values():
                if not temp.is_valid:
                    logger.critical(
                        "Invalid template <{}> does not have orders to skip", temp
                    )
                    continue
                bad_orders.union(temp.bad_orders)
        else:
            if not self.templates[subInst].is_valid:
                raise BadTemplateError
            bad_orders = self.templates[subInst].bad_orders

        return bad_orders

    def load_templates_from_disk(self) -> None:
        """Load stellar templates from disk.

        Parameters
        ----------
        path : str
            [description]

        """
        super().load_templates_from_disk()

    @property
    def RV_keyword(self) -> str:
        """RV keyword to use when aligning observations."""
        if self._internal_configs["ALIGNEMENT_RV_SOURCE"] == "SBART":
            RV_KW_start = "previous_SBART_RV"
        else:
            RV_KW_start = "DRS_RV"

        return RV_KW_start

model_type = 'Stellar' class-attribute instance-attribute

template_map: dict[STELLAR_CREATION_MODE, StellarTemplate] = {STELLAR_CREATION_MODE.Sum: SumStellar, STELLAR_CREATION_MODE.OBSERVATION: OBS_Stellar, STELLAR_CREATION_MODE.Median: MedianStellar, STELLAR_CREATION_MODE.PHOENIX: PHOENIX} class-attribute

iteration_number = 0 instance-attribute

RV_source = None instance-attribute

RV_keyword: str property

RV keyword to use when aligning observations.

Generate_Model(dataClass: DataClass | None, template_configs: Optional[UI_DICT] = None, conditions: Optional[ConditionModel] = None, force_computation: bool = False, store_templates: bool = True, previous_sbart_rv=None) -> None

Apply the spectral conditions to decide which observations to use.

Then, returns to the model generation as defined in the parent implementation.

Parameters:

Name Type Description Default
dataClass [type]

DataClass with the observations

 required
template_configs dict

Dictionary with the parameters needed to control the creation of the template. Each one has its own set of requirements

 None
conditions Optional[ConditionModel]

Either None (to not select the OBS that will be used) or :py:mod:~ASTRAutils.spectral_conditions to restrict the observations in use.

 None
force_computation bool

If True, recompute the stellar templates, even if they exist on disk. By default False

 False
store_templates bool

If True [default], store the templates to disk

 True
previous_sbart_rv

Object that s-BART loaded from disk, containing the RVs from a previous iteration

 None
Notes
  • The conditions that are passed to the StellarModel are only used for the creation of the stellar template. This will not reject observations from the RV extraction
Source code in src/ASTRA/template_creation/StellarModel.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
def Generate_Model(
    self,
    dataClass: DataClass | None,
    template_configs: Optional[UI_DICT] = None,
    conditions: Optional[ConditionModel] = None,
    force_computation: bool = False,
    store_templates: bool = True,
    previous_sbart_rv=None,
) -> None:
    """Apply the spectral conditions to decide which observations to use.

    Then, returns to the model generation as defined in the parent implementation.

    Parameters
    ----------
    dataClass : [type]
        DataClass with the observations
    template_configs : dict
        Dictionary with the parameters needed to control the creation of the template. Each one has its own
        set of requirements
    conditions: None, Condition
        Either None (to not select the OBS that will be used) or :py:mod:`~ASTRAutils.spectral_conditions` to
            restrict the observations in use.
    force_computation: bool
        If True, recompute the stellar templates, even if they exist on disk. By default False
    store_templates: bool
        If True [default], store the templates to disk
    previous_sbart_rv: SBART.data_objects.RV_outputs.RV_holder
        Object that s-BART loaded from disk, containing the RVs from a previous iteration

    Notes
    -----
    * The conditions that are passed to the StellarModel are **only** used for the creation of the stellar template.
     This will **not** reject observations from the RV extraction

    """
    if conditions is not None:
        logger.info("Applying conditions to creation of stellar template")
        self._creation_conditions = conditions

    self.RV_source = self._internal_configs["ALIGNEMENT_RV_SOURCE"]

    if self._internal_configs["ALIGNEMENT_RV_SOURCE"] == "SBART":
        logger.info(
            "{} using {} RVs as the source for stellar template creation",
            self.name,
            self._internal_configs["ALIGNEMENT_RV_SOURCE"],
        )
        if previous_sbart_rv is None:
            msg = "No previous SBART RVs were provided!"
            logger.critical(msg)
            raise custom_exceptions.InvalidConfiguration(msg)

        self.iteration_number = previous_sbart_rv.iteration_number + 1
        self.RV_source = "SBART"

        if dataClass is not None:
            try:
                dataClass.load_previous_SBART_results(
                    previous_sbart_rv,
                    use_merged_cube=self._internal_configs["USE_MERGED_RVS"],
                )
            except Exception as e:
                logger.opt(exception=True).critical("Failed")
                raise e
        else:
            logger.warning(
                "No DataClass provided. The RVs will not be loaded into the observations"
            )

    else:
        logger.info(
            "Using CCF RVs as the basis for the creation of the stellar models"
        )

    self.add_relative_path("Stellar", f"Stellar/Iteration_{self.iteration_number}")

    super().Generate_Model(
        dataClass=dataClass,
        template_configs=template_configs,
        attempt_to_load=not force_computation,
        store_templates=False,
    )

    for subInst, temp in self.templates.items():
        temp.update_RV_source_info(
            iteration_number=self.iteration_number,
            RV_source=self.RV_source,
            merged_source=self._internal_configs["USE_MERGED_RVS"],
        )

    if store_templates:
        self.store_templates_to_disk()

initialize_modelling_interfaces(new_properties)

Source code in src/ASTRA/template_creation/StellarModel.py 
186
187
188
def initialize_modelling_interfaces(self, new_properties):
    for subInst, temp in self.templates.items():
        temp.initialize_modelling_interface()

get_interpol_modes() -> set[str]

Source code in src/ASTRA/template_creation/StellarModel.py 
190
191
def get_interpol_modes(self) -> set[str]:
    return set(temp.interpol_mode for temp in self.templates.values())

get_orders_to_skip(subInst: str) -> set[int]

Compute orders to skip on a given subInstrument.

GIven by the orders that the stellar template rejected.

Args: subInst (str): sub-Instrument name

Raises: BadTemplateError: No valid template for the subInstrument

Returns: Set[int]: Set of bad orders.

Source code in src/ASTRA/template_creation/StellarModel.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
def get_orders_to_skip(self, subInst: str) -> set[int]:
    """Compute orders to skip on a given subInstrument.

    GIven by the orders that the stellar template rejected.

    Args:
        subInst (str): sub-Instrument name

    Raises:
        BadTemplateError: No valid template for the subInstrument

    Returns:
        Set[int]: Set of bad orders.

    """
    if subInst == "all":
        bad_orders: set[int] = set()
        for temp in self.templates.values():
            if not temp.is_valid:
                logger.critical(
                    "Invalid template <{}> does not have orders to skip", temp
                )
                continue
            bad_orders.union(temp.bad_orders)
    else:
        if not self.templates[subInst].is_valid:
            raise BadTemplateError
        bad_orders = self.templates[subInst].bad_orders

    return bad_orders

load_templates_from_disk() -> None

Load stellar templates from disk.

Parameters:

Name Type Description Default
path str

[description]

 required
Source code in src/ASTRA/template_creation/StellarModel.py 
250
251
252
253
254
255
256
257
258
259
def load_templates_from_disk(self) -> None:
    """Load stellar templates from disk.

    Parameters
    ----------
    path : str
        [description]

    """
    super().load_templates_from_disk()

Different models for the stellar template.

OBS_stellar

OBS_Stellar

Bases: StellarTemplate

Stellar template from the observation with the highest SNR (computed as a sum over all spectral orders).

Source code in src/ASTRA/template_creation/stellar_templates/OBS_stellar.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
class OBS_Stellar(StellarTemplate):
    """Stellar template from the observation with the highest SNR (computed as a sum over all spectral orders)."""

    method_name = "OBSERVATION"
    _default_params = StellarTemplate._default_params + DefaultValues(
        ALIGNEMENT_RV_SOURCE=UserParam(
            "DRS", constraint=ValueFromIterable(["DRS", "SBART"])
        ),
    )

    _default_params.update(
        "MINIMUM_NUMBER_OBS",
        UserParam(1, constraint=NumericValue, mandatory=False),
    )

    def __init__(
        self, subInst: str, user_configs: Union[None, dict] = None, loaded: bool = False
    ):
        super().__init__(subInst=subInst, user_configs=user_configs, loaded=loaded)
        self._selected_frameID = None
        self._found_error = False

    def _run_template_construction(
        self,
        dataClass,
        conditions,
        reference_frame,
    ) -> None:
        """Create the stellar template."""
        # removal may change the first common wavelength; make sure

        logger.info("Searching for frameID with highest sum of orderwise SNRs")
        total_SNR = []
        for frameID in self.frameIDs_to_use:
            total_SNR.append(
                np.nansum(
                    dataClass.get_KW_from_frameID(KW="orderwise_SNRs", frameID=frameID)
                )
            )

        self._selected_frameID = self.frameIDs_to_use[np.argmax(total_SNR)]
        logger.info("Selected frameID={}", self._selected_frameID)
        wavelenghts, spectra, uncertainties, mask = dataClass.get_frame_arrays_by_ID(
            self._selected_frameID
        )

        if self._internal_configs["OVERSAMPLE_TEMPLATE"] > 1:
            raise custom_exceptions.InternalError(
                "OBS stellar does not yet support oversampling"
            )
            # new_wave = np.linspace(
            #     wave_order[0],
            #     wave_order[-1],
            #     wavelenghts.shape[1] * self._internal_configs["OVERSAMPLE_TEMPLATE"],
            # )
            # self.wavelengths[order_index] = new_wave

        else:
            self.wavelengths = wavelenghts
            self.spectra = spectra
            self.uncertainties = uncertainties

        self.spectral_mask = Mask(mask, mask_type="binary")

        instrument_information = dataClass.get_instrument_information()
        epoch_shape = instrument_information["array_size"]
        self.rejection_array = np.zeros((1, epoch_shape[0]))

        self._finish_template_creation()

method_name = 'OBSERVATION' class-attribute instance-attribute

PHOENIX_STELLAR

PHOENIX

Bases: StellarTemplate

Stellar template from a PHOENIX template.

Source code in src/ASTRA/template_creation/stellar_templates/PHOENIX_STELLAR.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
class PHOENIX(StellarTemplate):
    """Stellar template from a PHOENIX template."""

    method_name = "PHOENIX"
    _default_params = StellarTemplate._default_params + DefaultValues(
        PHOENIX_FILE_PATH=UserParam("", constraint=PathValue, mandatory=False),
    )

    def __init__(
        self, subInst: str, user_configs: Union[None, dict] = None, loaded: bool = False
    ) -> None:
        super().__init__(subInst=subInst, user_configs=user_configs, loaded=loaded)
        self._found_error = False

    def _run_template_construction(
        self,
        dataClass,
        conditions,
        reference_frame,
    ) -> None:
        """Createthe stellar template."""
        # removal may change the first common wavelength; make sure

        logger.info("Searching for frameID with highest sum of orderwise SNRs")

        with fits.open(self._internal_configs["PHOENIX_FILE_PATH"]) as hdu:
            data = hdu[1].data

        logger.info("Selected frameID={}", self._selected_frameID)

        self.wavelengths = data["wavelengths"]
        self.spectra = data["flux"] * 10_000
        self.uncertainties = np.zeros(self.spectra.shape)
        mask = np.zeros_like(self.spectra, dtype=bool)
        self.spectral_mask = Mask(mask, mask_type="binary")

        self._finish_template_creation()

    def store_metrics(self): ...

method_name = 'PHOENIX' class-attribute instance-attribute

store_metrics()

Source code in src/ASTRA/template_creation/stellar_templates/PHOENIX_STELLAR.py 
61
def store_metrics(self): ...

Stellar_Template

Base class for stellar template.

StellarTemplate

Bases: BaseTemplate, Spectral_Modelling

Parent class of all StellarTemplates.

Implementing a common set of user configurations, data storage operations, among others.

User parameters:

========================= ================ ================ ================================= ================ Parameter name Mandatory Default Value Valid Values Comment ========================= ================ ================ ================================= ================ NUMBER_WORKERS False 1 Integer >= 0 [2] MEMORY_SAVE_MODE False False boolean [3] MINIMUM_NUMBER_OBS False 3 Integer >= 0 [4] ========================= ================ ================ ================================= ================

[1] - How to propagate the spectral uncertainties [2] - Number of jobs at once [3] - Save RAM by clearing the frame's S2D arrays from memory after using them [4] - Minimum number of valid observations needed to proceed with template creation

.. note:: This class also uses the User parameters defined by the :class:~ASTRAComponents.Modelling.Spectral_Modelling class

.. note:: Also check the User parameters of the parent classes for further customization options of SBART

Subclassed by:
Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
class StellarTemplate(BaseTemplate, Spectral_Modelling):
    """Parent class of all StellarTemplates.

    Implementing a common set of user configurations, data storage
    operations, among others.

    **User parameters:**

    ========================= ================ ================ ================================= ================
    Parameter name               Mandatory      Default Value               Valid Values            Comment
    ========================= ================ ================ ================================= ================
    NUMBER_WORKERS                 False           1            Integer >= 0                       [2]
    MEMORY_SAVE_MODE               False           False           boolean                         [3]
    MINIMUM_NUMBER_OBS             False           3            Integer >= 0                       [4]
    ========================= ================ ================ ================================= ================

    [1] - How to propagate the spectral uncertainties
    [2] - Number of jobs at once
    [3] - Save RAM by clearing the frame's S2D arrays from memory after using them
    [4] - Minimum number of  **valid** observations needed to proceed with template creation

    .. note::
       This class also uses the User parameters defined by the :class:`~ASTRAComponents.Modelling.Spectral_Modelling`
    class

    .. note::
        Also check the **User parameters** of the parent classes for further customization options of SBART

    """

    # Note: the root_path of the stellarTemplate is the templates/Stellar folder
    _name = "Stellar"

    _default_params = BaseTemplate._default_params + DefaultValues(
        CONSTANT_RV_GUESS=UserParam(
            False,
            BooleanValue,
            description="If set to True, the RV for the template construction will be set to zero for all observations. By default is False",
        ),
        NUMBER_WORKERS=UserParam(1, IntegerValue + Positive_Value_Constraint),
        MEMORY_SAVE_MODE=UserParam(
            False, constraint=BooleanValue
        ),  # if True, close the S2D files after using them!
        MINIMUM_NUMBER_OBS=UserParam(
            3, constraint=IntegerValue
        ),  # minimum number of OBS to create stellar template
        OVERSAMPLE_TEMPLATE=UserParam(
            default_value=1,
            constraint=ValueInInterval(interval=[1, 10000], include_edges=True)
            + ValueFromDtype(dtype_list=(int,)),
            description="If different than one, oversample the observations by this amount",
        ),
    )

    template_type = "Stellar"
    method_name = "Base"

    def __init__(
        self, subInst: str, user_configs: Union[None, dict] = None, loaded: bool = False
    ):
        super().__init__(subInst, user_configs, loaded)

        self.rejection_array = None

        self.frameIDs_to_use: list[int] = []

        self._in_shared_mem = False
        self.shm = {}
        self.package_pool = None
        self.output_pool = None
        self.target_name = None
        self._reference_frameID = None
        self._reference_filepath = None
        self.used_fpaths = []
        self._conditions = None

        self._rejection_flags_map = {}

        # Count the current iteration number; Used when improving the stellar template
        self._iter_number: int = 0
        self._RV_source = None
        self._merged_source = None
        self.sourceRVs: Optional[List] = None
        self.header_information = {}
        self.wave_air = None

    #################################
    #           Template creation   #
    #################################

    @custom_exceptions.ensure_invalid_template
    def create_stellar_template(
        self,
        dataClass,
        conditions,
        reference_frame: None | Frame = None,
    ) -> None:
        """Trigger the start of the creatoin of the stellar template.

        Parameters
        ----------
        dataClass : :class:`~ASTRAdata_objects.DataClass.DataClass`
            Data
        conditions: :py:mod:`~ASTRAutils.spectral_conditions`
            Apply conditions over the frames that will be used to create the stellar template
        reference_frame (None | Frame): if not None, ASTRA frame that will be used as the wavelength reference

        Raises
        ------
        StopComputationError
            If the computations is triggered by some event

        """
        logger.info(
            "Starting creation of {} template from {}",
            self.__class__.template_type,
            self._associated_subInst,
        )

        array_size = dataClass.get_instrument_information()["array_size"]
        self.array_size = [
            array_size[0],
            array_size[1] * self._internal_configs["OVERSAMPLE_TEMPLATE"],
        ]

        self._OrderStatus = OrderStatus(array_size[0])
        try:
            self.frameIDs_to_use = dataClass.get_frameIDs_from_subInst(
                self._associated_subInst
            )
        except NoDataError:
            logger.critical(
                "{} has no valid observations. Not computing {} template",
                self._associated_subInst,
                self.__class__.template_type,
            )
            self.add_to_status(
                MISSING_DATA(f"No valid observations from {self._associated_subInst}")
            )
            return

        self._base_checks_for_template_creation()

        self._conditions = conditions

        if conditions is None:
            logger.info(
                "No conditions for the creation of the stellar template. Using all available frames from {}",
                self._associated_subInst,
            )
            IDS_to_use = self.frameIDs_to_use
        else:
            logger.info(
                "Evaluating spectral conditions to select the valid observations from {}",
                self._associated_subInst,
            )
            IDS_to_use = []
            for frameID in self.frameIDs_to_use:
                keep, flags = conditions.evaluate(dataClass.get_frame_by_ID(frameID))
                if keep:
                    IDS_to_use.append(frameID)
                else:
                    logger.info(
                        "FrameID {} rejected from stellar template due to {}",
                        frameID,
                        flags,
                    )
                    self._rejection_flags_map[
                        dataClass.get_filename_from_frameID(frameID, full_path=True)
                    ] = flags

            if len(IDS_to_use) == 0:
                msg = (
                    f"subInstrument {self._associated_subInst} has no available observations after applying the spectral conditions. Stellar "
                    "template creation will halt"
                )

                logger.warning(msg)
                self.add_to_status(MISSING_DATA(msg))

        self._base_checks_for_template_creation()

        self.frameIDs_to_use = IDS_to_use

        if len(self.frameIDs_to_use) < self._internal_configs["MINIMUM_NUMBER_OBS"]:
            logger.critical(
                "Construction of stellar template from {} using less observations ({}) than the limit ({})",
                self._associated_subInst,
                len(self.frameIDs_to_use),
                self._internal_configs["MINIMUM_NUMBER_OBS"],
            )
            self.add_to_status(
                MISSING_DATA(
                    "Can't create stellar template with less than {} observations".format(
                        self._internal_configs["MINIMUM_NUMBER_OBS"],
                    ),
                ),
            )

        self._base_checks_for_template_creation()

        # TODO: ensure that they all observations are consistent!
        first_frame = dataClass.get_frame_by_ID(self.frameIDs_to_use[0])
        self.is_blaze_corrected = first_frame.check_if_data_correction_enabled(
            "is_blaze_corrected"
        )
        self.was_telluric_corrected = first_frame.check_if_data_correction_enabled(
            "was_telluric_corrected"
        )
        self.is_skysub = first_frame.is_skysub
        self.is_BERV_corrected = first_frame.check_if_data_correction_enabled(
            "is_BERV_corrected"
        )
        self.flux_atmos_balance_corrected = (
            first_frame.check_if_data_correction_enabled("flux_atmos_balance_corrected")
        )
        self.flux_dispersion_balance_corrected = (
            first_frame.check_if_data_correction_enabled(
                "flux_dispersion_balance_corrected",
            )
        )
        logger.info(f"Collecting RVs for the stellar tempalte with {self.RV_keyword=}")
        if self._internal_configs["CONSTANT_RV_GUESS"]:
            logger.warning(
                "Setting initial guess for template alignement to zero! Thread carefully in those uncharted areas",
            )
            self.sourceRVs = [0 * kilometer_second for _ in self.frameIDs_to_use]
        else:
            self.sourceRVs = dataClass.collect_RV_information(
                KW=self.RV_keyword,
                frameIDs=self.frameIDs_to_use,
                subInst=self._associated_subInst,
                units=kilometer_second,
                as_value=False,
                include_invalid=False,
            )
        self._run_template_construction(
            dataClass,
            conditions=conditions,
            reference_frame=reference_frame,
        )
        self.post_process_template(dataClass=dataClass)

    def _run_template_construction(
        self,
        dataClass: DataClass,
        conditions,
        reference_frame,
    ) -> None:
        raise NotImplementedError

    def post_process_template(self, dataClass: DataClass) -> None:
        # Computes weighted BJD, SNR and other metrics for the template
        # new_bjd = np.average(dataClass.)
        bjds = [
            dataClass.get_KW_from_frameID("BJD", frameID)
            for frameID in self.frameIDs_to_use
        ]

        errs = [
            dataClass.get_KW_from_frameID("DRS_RV_ERR", frameID)
            for frameID in self.frameIDs_to_use
        ]
        errs = convert_data(errs, new_units=kilometer_second, as_value=True)
        weight = np.array([1 / (err**2) for err in errs])
        weighted_bjds = np.average(bjds, weights=weight)

        self.header_information["STACKED BJD"] = weighted_bjds
        self.header_information["N OBS"] = len(self.frameIDs_to_use)
        for index, flux in enumerate(self.spectra):
            mask = ~self.spectral_mask.get_custom_mask()[index]
            self.header_information[f"ORDER {index} MEDIAN FLUX"] = np.nan_to_num(
                np.nanmedian(flux[mask]), nan=0
            )
        for index, flux in enumerate(self.spectra):
            mask = ~self.spectral_mask.get_custom_mask()[index]

            self.header_information[f"ORDER {index} SNR"] = np.nan_to_num(
                np.nanmedian(flux[mask] / self.uncertainties[index][mask]), nan=0
            )

        frame = dataClass.get_frame_by_ID(self._reference_frameID)
        self.header_information["POET_APERTURE"] = frame.get_KW_value("POET_APERTURE")
        self.header_information["MAX BERV"] = convert_data(
            frame.get_KW_value("MAX_BERV"),
            new_units=kilometer_second,
            as_value=True,
        )
        try:
            self.wave_air = remove_RVshift(
                frame.get_air_wavelength(),
                stellar_RV=convert_data(
                    frame.get_KW_value(self.RV_keyword),
                    new_units=kilometer_second,
                    as_value=True,
                ),
            )
        except Exception as e:
            logger.critical(f"Couldn't get the air wavelengths due to {e}")

    def add_new_frame_to_template(self, frame: Frame):
        """Allow to inject a new observation into a pre-existing model.

        This base function checks for
        a match on the different flux corrections and ensures that the loaded Flag is set to False,
        so that it is possible to update the disk products afterwards.

        Parameters
        ----------
        frame: Frame
            A new frame to inject into the stellar template

        Raises
        ------
        custom_exceptions.InvalidConfiguration:
            If the flux corrections of the Frame do not match those from the stellar template

        """
        logger.info(
            "Adding new frame to pre-existing stellar template. Updating model!"
        )
        self._loaded = False

        keep = True
        for name, val1, val2 in [
            (
                "Blaze",
                self.is_blaze_corrected,
                frame.check_if_data_correction_enabled("is_blaze_corrected"),
            ),
            (
                "Telluric",
                self.was_telluric_corrected,
                frame.check_if_data_correction_enabled("was_telluric_corrected"),
            ),
            (
                "BERV_corrected",
                self.is_BERV_corrected,
                frame.check_if_data_correction_enabled("is_BERV_corrected"),
            ),
            (
                "Flux atmos balance",
                self.flux_atmos_balance_corrected,
                frame.check_if_data_correction_enabled("flux_atmos_balance_corrected"),
            ),
            (
                "flux_dispersion_balance_corrected",
                self.flux_dispersion_balance_corrected,
                frame.check_if_data_correction_enabled(
                    "flux_dispersion_balance_corrected"
                ),
            ),
            ("sub-Instrument", self.sub_instrument, frame.sub_instrument),
        ]:
            if val1 != val2:
                keep = False
                logger.warning(
                    f"Template-frame corrections are different: {name} - template: {val1} - Frame: {val2}"
                )

        if not keep:
            msg = "New frame does not match the corrections from the stellar template"
            logger.critical(msg)
            raise custom_exceptions.InvalidConfiguration(msg)

    def evaluate_bad_orders(self) -> None:
        """Compute orders with too many points masked."""
        entire_mask = self.spectral_mask.get_custom_mask()
        N_pixels_per_order = self.pixels_per_order
        order_cutoff = N_pixels_per_order - 100

        for order in range(self.N_orders):
            order_mask = entire_mask[order]
            bad_pixels = np.sum(order_mask)
            if bad_pixels > order_cutoff:
                msg = f"Stellar template rejecting order {order} due to having more than {order_cutoff}/{N_pixels_per_order} pixels masked"

                logger.warning(msg)

                self._OrderStatus.add_flag_to_order(order, HIGH_CONTAMINATION(msg))

    #################################
    #           Data access        #
    #################################

    def check_if_used_frameID(self, frameID: int) -> bool:
        """Return true if the frameID was used for the construction of the template."""
        return frameID in self.frameIDs_to_use

    #################################
    #           Data storage        #
    #################################
    def trigger_data_storage(self, clobber: bool) -> None:  # noqa: D102
        try:
            super().trigger_data_storage(clobber)
            self._store_json_information()
            if self.disk_save_level != DISK_SAVE_MODE.EXTREME:
                self.store_metrics()
        except custom_exceptions.FailedStorage:
            return

    def _store_json_information(self) -> None:
        logger.info("Storing template flags to disk")
        storage_path = self._internalPaths.get_path_to("RunTimeProds", as_posix=False)

        detailedflags_storage = build_filename(
            storage_path,
            f"DetailedFlags_{self._associated_subInst}",
            fmt="json",
        )
        self._OrderStatus.store_as_json(detailedflags_storage)

        miscinfo = build_filename(
            storage_path,
            f"miscInfo_{self._associated_subInst}",
            fmt="json",
        )

        # Note: when adding things here we must use as a key the exact name of the parameter...
        # Note2: carefull with the datatypes that are stored in here...
        # Note3: The miscInfo file will use the parameters of this dict on a setattr

        with open(miscinfo, mode="w") as file:
            json.dump(self.get_miscInfo_of_template(), file, indent=4)
        logger.debug("Finished storage of template flags to disk")

    def get_miscInfo_of_template(self) -> Dict[str, Any]:
        """Get miscellaneous information from the stellar template, to be stored in output directory.

        Returns:
            Dict[str, Any]: Dictonary with key, value information (in json-compatible format) to be stored

        """
        return {
            "frameIDs_to_use": self.frameIDs_to_use,
            "used_fpaths": [i.as_posix() for i in self.used_fpaths],
            "is_BERV_corrected": self.is_BERV_corrected,
            "is_blaze_corrected": self.is_blaze_corrected,
            "was_telluric_corrected": self.was_telluric_corrected,
            "flux_dispersion_balance_corrected": self.flux_dispersion_balance_corrected,
            "flux_atmos_balance_corrected": self.flux_atmos_balance_corrected,
            "_reference_frameID": self._reference_frameID,
            "_reference_filepath": self._reference_filepath,
            "sourceRVs": convert_data(
                self.sourceRVs, new_units=kilometer_second, as_value=True
            ),
            "RV_keyword": self.RV_keyword,
        }

    def store_template(self, clobber: bool) -> None:
        # TODO: store the frameIDS that were used!
        logger.info("Starting template storage to disk")
        super().store_template(clobber)

        header = fits.Header()

        header["TYPE"] = self.__class__.template_type
        header["subInst"] = self._associated_subInst
        header["VERSION"] = __version__
        for key, config_val in self._internal_configs.items():
            if key in ["SAVE_DISK_SPACE"]:
                continue

            if isinstance(config_val, Enum):
                config_val = config_val.name
            header[f"HIERARCH {key}"] = config_val

        for key, value in self.header_information.items():
            header[f"HIERARCH {key}"] = value

        for key, value in self.get_miscInfo_of_template().items():
            if not isinstance(value, Iterable) and not isinstance(value, str):
                value = str(value)
                true_key = key.upper().replace("_", " ").strip()
                header[f"HIERARCH {true_key}"] = value

        hdu = fits.PrimaryHDU(data=[], header=header)

        hdus_cubes = [hdu]

        hdu_wave = fits.ImageHDU(data=self.wavelengths, header=header, name="WAVE")

        hdu_temp = fits.ImageHDU(data=self.spectra, header=header, name="TEMP")

        mask = self.spectral_mask.get_custom_mask().astype(int)

        hdu_mask = fits.ImageHDU(data=mask, header=header, name="MASK")
        hdu_uncerts = fits.ImageHDU(
            data=self.uncertainties, header=header, name="UNCERTAINTIES"
        )

        for val in [hdu_wave, hdu_temp, hdu_mask, hdu_uncerts]:
            hdus_cubes.append(val)

        if self.wave_air is not None:
            hdu_wave_air = fits.ImageHDU(
                data=self.wave_air, header=header, name="WAVE_AIR"
            )
            hdus_cubes.append(hdu_wave_air)
        hdul = fits.HDUList(hdus_cubes)

        filename = f"{self.storage_name}_{self._associated_subInst}.fits"
        hdul.writeto(self._internalPaths.root_storage_path / filename, overwrite=True)

        filename = f"{self.storage_name}_{self._associated_subInst}_inputs.txt"

        logger.debug("Storing used filepaths to disk")
        with open(
            self._internalPaths.root_storage_path / filename, mode="w"
        ) as to_write:
            if self._conditions is not None:
                self._conditions.write_to_disk(to_write)

            to_write.write(
                self._internal_configs.text_pretty_description(indent_level=0)
            )

            to_write.write(f"\n\nRejected files (N = {len(self._rejection_flags_map)})")
            for filepath, flags in self._rejection_flags_map.items():
                to_write.write(f"\n\t{filepath}:")
                for flag in flags:
                    to_write.write(f"\n\t\t{flag}")

            to_write.write(
                f"\n\nEpochs in use (Total = {len(self.frameIDs_to_use)}): \n"
            )
            for path in self.used_fpaths:
                to_write.write(f"\n{path}")
        logger.info("Finished template storage to disk")

    def store_metrics(self) -> None:
        metrics_path = self._internalPaths.get_path_to("metrics", as_posix=False)
        mask = self.spectral_mask.get_custom_mask().astype(int)

        N_order, N_pixel = self.spectra.shape
        fig, ax = plt.subplots(1, 1)
        custom_cmap = lambda x: "green" if x < 0.2 else "orange" if x < 0.6 else "red"
        for order in range(N_order):
            order_mask = mask[order]
            value = np.sum(order_mask) / N_pixel
            ax.scatter(order, value, color=custom_cmap(value))
        ax.set_xlabel("Order number")
        ax.set_ylabel("% pixels rejected")
        fig.savefig(
            metrics_path
            / f"Template_rejection_percentage_{self._associated_subInst}.png"
        )
        plt.close(fig)

        c = [
            "palegreen",
            "green",
            "darkgreen",
            "orange",
            "chocolate",
            "salmon",
            "red",
            "darkred",
            "purple",
            "black",
        ]
        v = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1]
        l = list(zip(v, c))
        cmap = LinearSegmentedColormap.from_list("rg", l, N=256)

        hori_dire = range(self.rejection_array.shape[1] + 1)
        verti_dire = range(self.rejection_array.shape[0] + 1)

        fig, ax = plt.subplots(figsize=(20, 10), constrained_layout=True)
        data = ax.pcolormesh(
            hori_dire,
            verti_dire,
            self.rejection_array,
            vmin=0,
            vmax=1,
            cmap=cmap,
            shading="auto",
            edgecolors="w",
        )
        fig.colorbar(data)
        ax.set_xlabel("Order")
        ax.set_ylabel("FrameID")
        ax.set_xticks(hori_dire[::10])

        ax.set_yticks(np.asarray(verti_dire)[:-1] + 0.5, labels=self.frameIDs_to_use)
        # ax.set_yticks(verti_dire, minor=False, labels=frameIDs)

        ax.set_ylim([0, np.max(verti_dire)])
        ax.set_xlim([0, np.max(hori_dire)])

        new_array = np.c_[self.frameIDs_to_use, self.rejection_array]

        new_array = np.r_[[[np.nan] + list(hori_dire)[:-1]], new_array]

        np.savetxt(
            metrics_path / f"order_pixel_rejections_{self._associated_subInst}.txt",
            new_array,
            fmt="%.3f",
        )
        fig.savefig(
            metrics_path / f"order_pixel_rejection_{self._associated_subInst}.pdf"
        )
        plt.close(fig)

    def remove_telluric_features(self, Telluric_Template: TelluricModel) -> None:
        """Remove telluric features from the stellar template."""

        subInst_template = Telluric_Template.request_data(self.sub_instrument)

        order_map = {
            i: (np.min(self.wavelengths[i]), np.max(self.wavelengths[i]))
            for i in range(self.N_orders)
        }
        if subInst_template.for_feature_removal:
            for wave_pair in subInst_template.contaminated_regions:

                for order in range(self.N_orders):
                    if check_if_overlap(wave_pair, order_map[order]):
                        indexes = np.where(
                            np.logical_and(
                                self.wavelengths[order] >= wave_pair[0],
                                self.wavelengths[order] <= wave_pair[1],
                            ),
                        )
                        self.spectral_mask.add_indexes_to_mask_order(
                            order, indexes, TELLURIC
                        )
        else:
            raise custom_exceptions.InternalError("Why aren't we removing tellurics?")
        self._applied_telluric_removal = True

    def load_from_file(self, root_path, loading_path: str) -> None:
        super().load_from_file(root_path, loading_path)

        with fits.open(loading_path) as hdulist:
            if hdulist[1].header.get("VERSION", "") != __version__:
                logger.warning(
                    "Loaded template was not created under the current SBART version"
                )
            self._associated_subInst = hdulist["WAVE"].header["subInst"]

            self.spectra = hdulist["TEMP"].data
            self.wavelengths = hdulist["WAVE"].data
            self.uncertainties = hdulist["UNCERTAINTIES"].data
            self.wave_air = hdulist["WAVE_AIR"].data

            mask = hdulist["MASK"].data.astype(bool)
            self.spectral_mask = Mask(initial_mask=mask)
            header = hdulist[0].header

        for kw in ["MAX_BERV", "POET_APERTURE"]:
            self.header_information[kw] = header.get(kw, 0)

        self.array_size = self.spectra.shape
        storage_path = self._internalPaths.get_path_to("RunTimeProds", as_posix=False)

        detailedFlags_storage = build_filename(
            storage_path,
            f"DetailedFlags_{self._associated_subInst}",
            fmt="json",
        )

        try:
            self._OrderStatus = OrderStatus.load_from_json(detailedFlags_storage)
        except FileNotFoundError:
            raise custom_exceptions.TemplateNotExistsError(
                "No data for current ASTRA version"
            )

        miscInfo = build_filename(
            storage_path,
            f"miscInfo_{self._associated_subInst}",
            fmt="json",
        )

        with open(miscInfo) as file:
            json_info = json.load(file)

        json_info["used_fpaths"] = list([Path(i) for i in json_info["used_fpaths"]])

        for key, value in json_info.items():
            if key == "RV_keyword":
                # This one is not supposed to be loaded in the current version
                continue
            setattr(self, key, value)

    # Handle shared memory
    def convert_to_shared_mem(
        self, custom_size: Optional[Iterable[float]] = None
    ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
        """Convert stellar template into shared memory arrays (used for construction of the model)."""
        if self._in_shared_mem:
            # Avoid opening multiple arrays!
            logger.info(f"{self.__class__.name} already in shared memory")
            return None
        logger.info("Putting the stellar template in shared memory")
        self._in_shared_mem = True

        array_of_zeros = (
            np.zeros(self.wavelengths.shape)
            if custom_size is None
            else np.zeros(custom_size)
        )

        buffer_info, shared_uncerts = create_shared_array(array_of_zeros)
        self.shm["template_errors"] = buffer_info
        uncertainties = shared_uncerts

        buffer_info, shared_temp = create_shared_array(array_of_zeros)
        self.shm["template"] = buffer_info
        template = shared_temp

        buffer_info, shared_wave = create_shared_array(np.zeros(self.wavelengths.shape))
        self.shm["template_counts"] = buffer_info
        counts = shared_wave

        buffer_info, shared_wave = create_shared_array(self.wavelengths)
        self.shm["template_wavelength"] = buffer_info
        wavelengths = shared_wave

        return wavelengths, template, uncertainties, counts

    def cleanup_shared_memory(self) -> None:
        """Close shared memory interface (after template construction)."""
        logger.debug(
            "Cleaning the shared memory interfaces from the Stellar template creation"
        )
        self._close_workers()
        self._close_queues()
        self._close_shared_memory_arrays()

    def _close_workers(self) -> None:
        logger.debug("{} closing the workers", self.name)
        for _ in range(self._internal_configs["NUMBER_WORKERS"]):
            self.package_pool.put(np.nan)

    def _close_shared_memory_arrays(self) -> None:
        """Close the shared memory arrays."""
        logger.debug("{} closing the shared memory array", self.name)

        self._in_shared_mem = False
        for mem_block in self.shm.values():
            mem_block[0].close()
            mem_block[0].unlink()

        self.shm = {}

    def _close_queues(self) -> None:
        """Close the communication queues."""
        logger.debug("{} closing the communication queues", self.name)

        self.package_pool.close()
        self.output_pool.close()

        self.package_pool = None
        self.output_pool = None

    @property
    def storage_name(self) -> str:
        """Return the storage name for the stellar templates.

        Will include the name of the method, the RV source for the aligmenet

        """
        base_name = f"{self.__class__.method_name}_{self._RV_source}"

        if self._merged_source is not None:
            base_name += "_"
            base_name += "merged" if self._merged_source else "individual"

        return base_name

    @property
    def iteration_number(self) -> int:
        """Get current iteration number."""
        return self._iter_number

    def update_RV_source_info(
        self, iteration_number: int, RV_source: str, merged_source: bool
    ):
        """Update information from the RV source.

        Args:
            iteration_number (int): Current iteration number
            RV_source (str): Name of the source
            merged_source (bool): merged RVs or not.

        """
        self._iter_number = iteration_number
        self._RV_source = RV_source
        self._merged_source = merged_source

    @property
    def RV_keyword(self) -> str:
        """RV keyword to use when aligning observations."""
        if self._internal_configs["ALIGNEMENT_RV_SOURCE"] == "SBART":
            RV_KW_start = "previous_SBART_RV"
        else:
            RV_KW_start = "DRS_RV"

        return RV_KW_start

template_type = 'Stellar' class-attribute instance-attribute

method_name = 'Base' class-attribute instance-attribute

rejection_array = None instance-attribute

frameIDs_to_use: list[int] = [] instance-attribute

shm = {} instance-attribute

package_pool = None instance-attribute

output_pool = None instance-attribute

target_name = None instance-attribute

used_fpaths = [] instance-attribute

sourceRVs: Optional[List] = None instance-attribute

header_information = {} instance-attribute

wave_air = None instance-attribute

storage_name: str property

Return the storage name for the stellar templates.

Will include the name of the method, the RV source for the aligmenet

iteration_number: int property

Get current iteration number.

RV_keyword: str property

RV keyword to use when aligning observations.

create_stellar_template(dataClass, conditions, reference_frame: None | Frame = None) -> None

Trigger the start of the creatoin of the stellar template.

Parameters:

Name Type Description Default
dataClass :class:`~ASTRAdata_objects.DataClass.DataClass`

Data

 required
conditions

Apply conditions over the frames that will be used to create the stellar template

 required
reference_frame None | Frame
None

Raises:

Type Description
StopComputationError

If the computations is triggered by some event
Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
@custom_exceptions.ensure_invalid_template
def create_stellar_template(
    self,
    dataClass,
    conditions,
    reference_frame: None | Frame = None,
) -> None:
    """Trigger the start of the creatoin of the stellar template.

    Parameters
    ----------
    dataClass : :class:`~ASTRAdata_objects.DataClass.DataClass`
        Data
    conditions: :py:mod:`~ASTRAutils.spectral_conditions`
        Apply conditions over the frames that will be used to create the stellar template
    reference_frame (None | Frame): if not None, ASTRA frame that will be used as the wavelength reference

    Raises
    ------
    StopComputationError
        If the computations is triggered by some event

    """
    logger.info(
        "Starting creation of {} template from {}",
        self.__class__.template_type,
        self._associated_subInst,
    )

    array_size = dataClass.get_instrument_information()["array_size"]
    self.array_size = [
        array_size[0],
        array_size[1] * self._internal_configs["OVERSAMPLE_TEMPLATE"],
    ]

    self._OrderStatus = OrderStatus(array_size[0])
    try:
        self.frameIDs_to_use = dataClass.get_frameIDs_from_subInst(
            self._associated_subInst
        )
    except NoDataError:
        logger.critical(
            "{} has no valid observations. Not computing {} template",
            self._associated_subInst,
            self.__class__.template_type,
        )
        self.add_to_status(
            MISSING_DATA(f"No valid observations from {self._associated_subInst}")
        )
        return

    self._base_checks_for_template_creation()

    self._conditions = conditions

    if conditions is None:
        logger.info(
            "No conditions for the creation of the stellar template. Using all available frames from {}",
            self._associated_subInst,
        )
        IDS_to_use = self.frameIDs_to_use
    else:
        logger.info(
            "Evaluating spectral conditions to select the valid observations from {}",
            self._associated_subInst,
        )
        IDS_to_use = []
        for frameID in self.frameIDs_to_use:
            keep, flags = conditions.evaluate(dataClass.get_frame_by_ID(frameID))
            if keep:
                IDS_to_use.append(frameID)
            else:
                logger.info(
                    "FrameID {} rejected from stellar template due to {}",
                    frameID,
                    flags,
                )
                self._rejection_flags_map[
                    dataClass.get_filename_from_frameID(frameID, full_path=True)
                ] = flags

        if len(IDS_to_use) == 0:
            msg = (
                f"subInstrument {self._associated_subInst} has no available observations after applying the spectral conditions. Stellar "
                "template creation will halt"
            )

            logger.warning(msg)
            self.add_to_status(MISSING_DATA(msg))

    self._base_checks_for_template_creation()

    self.frameIDs_to_use = IDS_to_use

    if len(self.frameIDs_to_use) < self._internal_configs["MINIMUM_NUMBER_OBS"]:
        logger.critical(
            "Construction of stellar template from {} using less observations ({}) than the limit ({})",
            self._associated_subInst,
            len(self.frameIDs_to_use),
            self._internal_configs["MINIMUM_NUMBER_OBS"],
        )
        self.add_to_status(
            MISSING_DATA(
                "Can't create stellar template with less than {} observations".format(
                    self._internal_configs["MINIMUM_NUMBER_OBS"],
                ),
            ),
        )

    self._base_checks_for_template_creation()

    # TODO: ensure that they all observations are consistent!
    first_frame = dataClass.get_frame_by_ID(self.frameIDs_to_use[0])
    self.is_blaze_corrected = first_frame.check_if_data_correction_enabled(
        "is_blaze_corrected"
    )
    self.was_telluric_corrected = first_frame.check_if_data_correction_enabled(
        "was_telluric_corrected"
    )
    self.is_skysub = first_frame.is_skysub
    self.is_BERV_corrected = first_frame.check_if_data_correction_enabled(
        "is_BERV_corrected"
    )
    self.flux_atmos_balance_corrected = (
        first_frame.check_if_data_correction_enabled("flux_atmos_balance_corrected")
    )
    self.flux_dispersion_balance_corrected = (
        first_frame.check_if_data_correction_enabled(
            "flux_dispersion_balance_corrected",
        )
    )
    logger.info(f"Collecting RVs for the stellar tempalte with {self.RV_keyword=}")
    if self._internal_configs["CONSTANT_RV_GUESS"]:
        logger.warning(
            "Setting initial guess for template alignement to zero! Thread carefully in those uncharted areas",
        )
        self.sourceRVs = [0 * kilometer_second for _ in self.frameIDs_to_use]
    else:
        self.sourceRVs = dataClass.collect_RV_information(
            KW=self.RV_keyword,
            frameIDs=self.frameIDs_to_use,
            subInst=self._associated_subInst,
            units=kilometer_second,
            as_value=False,
            include_invalid=False,
        )
    self._run_template_construction(
        dataClass,
        conditions=conditions,
        reference_frame=reference_frame,
    )
    self.post_process_template(dataClass=dataClass)

post_process_template(dataClass: DataClass) -> None

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
def post_process_template(self, dataClass: DataClass) -> None:
    # Computes weighted BJD, SNR and other metrics for the template
    # new_bjd = np.average(dataClass.)
    bjds = [
        dataClass.get_KW_from_frameID("BJD", frameID)
        for frameID in self.frameIDs_to_use
    ]

    errs = [
        dataClass.get_KW_from_frameID("DRS_RV_ERR", frameID)
        for frameID in self.frameIDs_to_use
    ]
    errs = convert_data(errs, new_units=kilometer_second, as_value=True)
    weight = np.array([1 / (err**2) for err in errs])
    weighted_bjds = np.average(bjds, weights=weight)

    self.header_information["STACKED BJD"] = weighted_bjds
    self.header_information["N OBS"] = len(self.frameIDs_to_use)
    for index, flux in enumerate(self.spectra):
        mask = ~self.spectral_mask.get_custom_mask()[index]
        self.header_information[f"ORDER {index} MEDIAN FLUX"] = np.nan_to_num(
            np.nanmedian(flux[mask]), nan=0
        )
    for index, flux in enumerate(self.spectra):
        mask = ~self.spectral_mask.get_custom_mask()[index]

        self.header_information[f"ORDER {index} SNR"] = np.nan_to_num(
            np.nanmedian(flux[mask] / self.uncertainties[index][mask]), nan=0
        )

    frame = dataClass.get_frame_by_ID(self._reference_frameID)
    self.header_information["POET_APERTURE"] = frame.get_KW_value("POET_APERTURE")
    self.header_information["MAX BERV"] = convert_data(
        frame.get_KW_value("MAX_BERV"),
        new_units=kilometer_second,
        as_value=True,
    )
    try:
        self.wave_air = remove_RVshift(
            frame.get_air_wavelength(),
            stellar_RV=convert_data(
                frame.get_KW_value(self.RV_keyword),
                new_units=kilometer_second,
                as_value=True,
            ),
        )
    except Exception as e:
        logger.critical(f"Couldn't get the air wavelengths due to {e}")

add_new_frame_to_template(frame: Frame)

Allow to inject a new observation into a pre-existing model.

This base function checks for a match on the different flux corrections and ensures that the loaded Flag is set to False, so that it is possible to update the disk products afterwards.

Parameters:

Name Type Description Default
frame Frame

A new frame to inject into the stellar template

 required

Raises:

Type Description
custom_exceptions.InvalidConfiguration:

If the flux corrections of the Frame do not match those from the stellar template
Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
def add_new_frame_to_template(self, frame: Frame):
    """Allow to inject a new observation into a pre-existing model.

    This base function checks for
    a match on the different flux corrections and ensures that the loaded Flag is set to False,
    so that it is possible to update the disk products afterwards.

    Parameters
    ----------
    frame: Frame
        A new frame to inject into the stellar template

    Raises
    ------
    custom_exceptions.InvalidConfiguration:
        If the flux corrections of the Frame do not match those from the stellar template

    """
    logger.info(
        "Adding new frame to pre-existing stellar template. Updating model!"
    )
    self._loaded = False

    keep = True
    for name, val1, val2 in [
        (
            "Blaze",
            self.is_blaze_corrected,
            frame.check_if_data_correction_enabled("is_blaze_corrected"),
        ),
        (
            "Telluric",
            self.was_telluric_corrected,
            frame.check_if_data_correction_enabled("was_telluric_corrected"),
        ),
        (
            "BERV_corrected",
            self.is_BERV_corrected,
            frame.check_if_data_correction_enabled("is_BERV_corrected"),
        ),
        (
            "Flux atmos balance",
            self.flux_atmos_balance_corrected,
            frame.check_if_data_correction_enabled("flux_atmos_balance_corrected"),
        ),
        (
            "flux_dispersion_balance_corrected",
            self.flux_dispersion_balance_corrected,
            frame.check_if_data_correction_enabled(
                "flux_dispersion_balance_corrected"
            ),
        ),
        ("sub-Instrument", self.sub_instrument, frame.sub_instrument),
    ]:
        if val1 != val2:
            keep = False
            logger.warning(
                f"Template-frame corrections are different: {name} - template: {val1} - Frame: {val2}"
            )

    if not keep:
        msg = "New frame does not match the corrections from the stellar template"
        logger.critical(msg)
        raise custom_exceptions.InvalidConfiguration(msg)

evaluate_bad_orders() -> None

Compute orders with too many points masked.

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
def evaluate_bad_orders(self) -> None:
    """Compute orders with too many points masked."""
    entire_mask = self.spectral_mask.get_custom_mask()
    N_pixels_per_order = self.pixels_per_order
    order_cutoff = N_pixels_per_order - 100

    for order in range(self.N_orders):
        order_mask = entire_mask[order]
        bad_pixels = np.sum(order_mask)
        if bad_pixels > order_cutoff:
            msg = f"Stellar template rejecting order {order} due to having more than {order_cutoff}/{N_pixels_per_order} pixels masked"

            logger.warning(msg)

            self._OrderStatus.add_flag_to_order(order, HIGH_CONTAMINATION(msg))

check_if_used_frameID(frameID: int) -> bool

Return true if the frameID was used for the construction of the template.

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
426
427
428
def check_if_used_frameID(self, frameID: int) -> bool:
    """Return true if the frameID was used for the construction of the template."""
    return frameID in self.frameIDs_to_use

trigger_data_storage(clobber: bool) -> None

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
433
434
435
436
437
438
439
440
def trigger_data_storage(self, clobber: bool) -> None:  # noqa: D102
    try:
        super().trigger_data_storage(clobber)
        self._store_json_information()
        if self.disk_save_level != DISK_SAVE_MODE.EXTREME:
            self.store_metrics()
    except custom_exceptions.FailedStorage:
        return

get_miscInfo_of_template() -> Dict[str, Any]

Get miscellaneous information from the stellar template, to be stored in output directory.

Returns: Dict[str, Any]: Dictonary with key, value information (in json-compatible format) to be stored

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
def get_miscInfo_of_template(self) -> Dict[str, Any]:
    """Get miscellaneous information from the stellar template, to be stored in output directory.

    Returns:
        Dict[str, Any]: Dictonary with key, value information (in json-compatible format) to be stored

    """
    return {
        "frameIDs_to_use": self.frameIDs_to_use,
        "used_fpaths": [i.as_posix() for i in self.used_fpaths],
        "is_BERV_corrected": self.is_BERV_corrected,
        "is_blaze_corrected": self.is_blaze_corrected,
        "was_telluric_corrected": self.was_telluric_corrected,
        "flux_dispersion_balance_corrected": self.flux_dispersion_balance_corrected,
        "flux_atmos_balance_corrected": self.flux_atmos_balance_corrected,
        "_reference_frameID": self._reference_frameID,
        "_reference_filepath": self._reference_filepath,
        "sourceRVs": convert_data(
            self.sourceRVs, new_units=kilometer_second, as_value=True
        ),
        "RV_keyword": self.RV_keyword,
    }

store_template(clobber: bool) -> None

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
def store_template(self, clobber: bool) -> None:
    # TODO: store the frameIDS that were used!
    logger.info("Starting template storage to disk")
    super().store_template(clobber)

    header = fits.Header()

    header["TYPE"] = self.__class__.template_type
    header["subInst"] = self._associated_subInst
    header["VERSION"] = __version__
    for key, config_val in self._internal_configs.items():
        if key in ["SAVE_DISK_SPACE"]:
            continue

        if isinstance(config_val, Enum):
            config_val = config_val.name
        header[f"HIERARCH {key}"] = config_val

    for key, value in self.header_information.items():
        header[f"HIERARCH {key}"] = value

    for key, value in self.get_miscInfo_of_template().items():
        if not isinstance(value, Iterable) and not isinstance(value, str):
            value = str(value)
            true_key = key.upper().replace("_", " ").strip()
            header[f"HIERARCH {true_key}"] = value

    hdu = fits.PrimaryHDU(data=[], header=header)

    hdus_cubes = [hdu]

    hdu_wave = fits.ImageHDU(data=self.wavelengths, header=header, name="WAVE")

    hdu_temp = fits.ImageHDU(data=self.spectra, header=header, name="TEMP")

    mask = self.spectral_mask.get_custom_mask().astype(int)

    hdu_mask = fits.ImageHDU(data=mask, header=header, name="MASK")
    hdu_uncerts = fits.ImageHDU(
        data=self.uncertainties, header=header, name="UNCERTAINTIES"
    )

    for val in [hdu_wave, hdu_temp, hdu_mask, hdu_uncerts]:
        hdus_cubes.append(val)

    if self.wave_air is not None:
        hdu_wave_air = fits.ImageHDU(
            data=self.wave_air, header=header, name="WAVE_AIR"
        )
        hdus_cubes.append(hdu_wave_air)
    hdul = fits.HDUList(hdus_cubes)

    filename = f"{self.storage_name}_{self._associated_subInst}.fits"
    hdul.writeto(self._internalPaths.root_storage_path / filename, overwrite=True)

    filename = f"{self.storage_name}_{self._associated_subInst}_inputs.txt"

    logger.debug("Storing used filepaths to disk")
    with open(
        self._internalPaths.root_storage_path / filename, mode="w"
    ) as to_write:
        if self._conditions is not None:
            self._conditions.write_to_disk(to_write)

        to_write.write(
            self._internal_configs.text_pretty_description(indent_level=0)
        )

        to_write.write(f"\n\nRejected files (N = {len(self._rejection_flags_map)})")
        for filepath, flags in self._rejection_flags_map.items():
            to_write.write(f"\n\t{filepath}:")
            for flag in flags:
                to_write.write(f"\n\t\t{flag}")

        to_write.write(
            f"\n\nEpochs in use (Total = {len(self.frameIDs_to_use)}): \n"
        )
        for path in self.used_fpaths:
            to_write.write(f"\n{path}")
    logger.info("Finished template storage to disk")

store_metrics() -> None

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
def store_metrics(self) -> None:
    metrics_path = self._internalPaths.get_path_to("metrics", as_posix=False)
    mask = self.spectral_mask.get_custom_mask().astype(int)

    N_order, N_pixel = self.spectra.shape
    fig, ax = plt.subplots(1, 1)
    custom_cmap = lambda x: "green" if x < 0.2 else "orange" if x < 0.6 else "red"
    for order in range(N_order):
        order_mask = mask[order]
        value = np.sum(order_mask) / N_pixel
        ax.scatter(order, value, color=custom_cmap(value))
    ax.set_xlabel("Order number")
    ax.set_ylabel("% pixels rejected")
    fig.savefig(
        metrics_path
        / f"Template_rejection_percentage_{self._associated_subInst}.png"
    )
    plt.close(fig)

    c = [
        "palegreen",
        "green",
        "darkgreen",
        "orange",
        "chocolate",
        "salmon",
        "red",
        "darkred",
        "purple",
        "black",
    ]
    v = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1]
    l = list(zip(v, c))
    cmap = LinearSegmentedColormap.from_list("rg", l, N=256)

    hori_dire = range(self.rejection_array.shape[1] + 1)
    verti_dire = range(self.rejection_array.shape[0] + 1)

    fig, ax = plt.subplots(figsize=(20, 10), constrained_layout=True)
    data = ax.pcolormesh(
        hori_dire,
        verti_dire,
        self.rejection_array,
        vmin=0,
        vmax=1,
        cmap=cmap,
        shading="auto",
        edgecolors="w",
    )
    fig.colorbar(data)
    ax.set_xlabel("Order")
    ax.set_ylabel("FrameID")
    ax.set_xticks(hori_dire[::10])

    ax.set_yticks(np.asarray(verti_dire)[:-1] + 0.5, labels=self.frameIDs_to_use)
    # ax.set_yticks(verti_dire, minor=False, labels=frameIDs)

    ax.set_ylim([0, np.max(verti_dire)])
    ax.set_xlim([0, np.max(hori_dire)])

    new_array = np.c_[self.frameIDs_to_use, self.rejection_array]

    new_array = np.r_[[[np.nan] + list(hori_dire)[:-1]], new_array]

    np.savetxt(
        metrics_path / f"order_pixel_rejections_{self._associated_subInst}.txt",
        new_array,
        fmt="%.3f",
    )
    fig.savefig(
        metrics_path / f"order_pixel_rejection_{self._associated_subInst}.pdf"
    )
    plt.close(fig)

remove_telluric_features(Telluric_Template: TelluricModel) -> None

Remove telluric features from the stellar template.

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
def remove_telluric_features(self, Telluric_Template: TelluricModel) -> None:
    """Remove telluric features from the stellar template."""

    subInst_template = Telluric_Template.request_data(self.sub_instrument)

    order_map = {
        i: (np.min(self.wavelengths[i]), np.max(self.wavelengths[i]))
        for i in range(self.N_orders)
    }
    if subInst_template.for_feature_removal:
        for wave_pair in subInst_template.contaminated_regions:

            for order in range(self.N_orders):
                if check_if_overlap(wave_pair, order_map[order]):
                    indexes = np.where(
                        np.logical_and(
                            self.wavelengths[order] >= wave_pair[0],
                            self.wavelengths[order] <= wave_pair[1],
                        ),
                    )
                    self.spectral_mask.add_indexes_to_mask_order(
                        order, indexes, TELLURIC
                    )
    else:
        raise custom_exceptions.InternalError("Why aren't we removing tellurics?")
    self._applied_telluric_removal = True

load_from_file(root_path, loading_path: str) -> None

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
def load_from_file(self, root_path, loading_path: str) -> None:
    super().load_from_file(root_path, loading_path)

    with fits.open(loading_path) as hdulist:
        if hdulist[1].header.get("VERSION", "") != __version__:
            logger.warning(
                "Loaded template was not created under the current SBART version"
            )
        self._associated_subInst = hdulist["WAVE"].header["subInst"]

        self.spectra = hdulist["TEMP"].data
        self.wavelengths = hdulist["WAVE"].data
        self.uncertainties = hdulist["UNCERTAINTIES"].data
        self.wave_air = hdulist["WAVE_AIR"].data

        mask = hdulist["MASK"].data.astype(bool)
        self.spectral_mask = Mask(initial_mask=mask)
        header = hdulist[0].header

    for kw in ["MAX_BERV", "POET_APERTURE"]:
        self.header_information[kw] = header.get(kw, 0)

    self.array_size = self.spectra.shape
    storage_path = self._internalPaths.get_path_to("RunTimeProds", as_posix=False)

    detailedFlags_storage = build_filename(
        storage_path,
        f"DetailedFlags_{self._associated_subInst}",
        fmt="json",
    )

    try:
        self._OrderStatus = OrderStatus.load_from_json(detailedFlags_storage)
    except FileNotFoundError:
        raise custom_exceptions.TemplateNotExistsError(
            "No data for current ASTRA version"
        )

    miscInfo = build_filename(
        storage_path,
        f"miscInfo_{self._associated_subInst}",
        fmt="json",
    )

    with open(miscInfo) as file:
        json_info = json.load(file)

    json_info["used_fpaths"] = list([Path(i) for i in json_info["used_fpaths"]])

    for key, value in json_info.items():
        if key == "RV_keyword":
            # This one is not supposed to be loaded in the current version
            continue
        setattr(self, key, value)

convert_to_shared_mem(custom_size: Optional[Iterable[float]] = None) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]

Convert stellar template into shared memory arrays (used for construction of the model).

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
def convert_to_shared_mem(
    self, custom_size: Optional[Iterable[float]] = None
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Convert stellar template into shared memory arrays (used for construction of the model)."""
    if self._in_shared_mem:
        # Avoid opening multiple arrays!
        logger.info(f"{self.__class__.name} already in shared memory")
        return None
    logger.info("Putting the stellar template in shared memory")
    self._in_shared_mem = True

    array_of_zeros = (
        np.zeros(self.wavelengths.shape)
        if custom_size is None
        else np.zeros(custom_size)
    )

    buffer_info, shared_uncerts = create_shared_array(array_of_zeros)
    self.shm["template_errors"] = buffer_info
    uncertainties = shared_uncerts

    buffer_info, shared_temp = create_shared_array(array_of_zeros)
    self.shm["template"] = buffer_info
    template = shared_temp

    buffer_info, shared_wave = create_shared_array(np.zeros(self.wavelengths.shape))
    self.shm["template_counts"] = buffer_info
    counts = shared_wave

    buffer_info, shared_wave = create_shared_array(self.wavelengths)
    self.shm["template_wavelength"] = buffer_info
    wavelengths = shared_wave

    return wavelengths, template, uncertainties, counts

cleanup_shared_memory() -> None

Close shared memory interface (after template construction).

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
763
764
765
766
767
768
769
770
def cleanup_shared_memory(self) -> None:
    """Close shared memory interface (after template construction)."""
    logger.debug(
        "Cleaning the shared memory interfaces from the Stellar template creation"
    )
    self._close_workers()
    self._close_queues()
    self._close_shared_memory_arrays()

update_RV_source_info(iteration_number: int, RV_source: str, merged_source: bool)

Update information from the RV source.

Args: iteration_number (int): Current iteration number RV_source (str): Name of the source merged_source (bool): merged RVs or not.

Source code in src/ASTRA/template_creation/stellar_templates/Stellar_Template.py 
818
819
820
821
822
823
824
825
826
827
828
829
830
831
def update_RV_source_info(
    self, iteration_number: int, RV_source: str, merged_source: bool
):
    """Update information from the RV source.

    Args:
        iteration_number (int): Current iteration number
        RV_source (str): Name of the source
        merged_source (bool): merged RVs or not.

    """
    self._iter_number = iteration_number
    self._RV_source = RV_source
    self._merged_source = merged_source

median_stellar

MedianStellar

Bases: StellarTemplate

Stellar template constructed as a median.

User parameters:

This object doesn't introduce unique user parameters.

Note: Check the User parameters of the parent classes for further customization options of this class

Source code in src/ASTRA/template_creation/stellar_templates/median_stellar.py 
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
class MedianStellar(StellarTemplate):
    """Stellar template constructed as a median.

    **User parameters:**

        This object doesn't introduce unique user parameters.

    *Note:* Check the **User parameters** of the parent classes for further customization options of this class

    """

    method_name = "Median"
    _default_params = StellarTemplate._default_params + DefaultValues(
        ALIGNEMENT_RV_SOURCE=UserParam(
            "DRS", constraint=ValueFromIterable(["DRS", "SBART"])
        ),
        FLUX_threshold_for_template=UserParam(
            default_value=1,
            constraint=Positive_Value_Constraint,
            description="Flux threshold for masking the spectral template. Set to one to avoid possible numerical issues with near-zero values",
        ),
    )

    def __init__(
        self, subInst: str, user_configs: Optional[Dict] = None, loaded: bool = False
    ):
        super().__init__(subInst=subInst, user_configs=user_configs, loaded=loaded)

        if not loaded:
            if self._internal_configs["MEMORY_SAVE_MODE"]:
                logger.warning(
                    "Stellar template creation will save RAM usage. This will result in multiple open/close "
                    "operations across the entire SBART pipeline! ",
                )

            self._found_error = False

    def _run_template_construction(
        self,
        dataClass,
        conditions,
        reference_frame,
    ) -> None:
        """Creating the stellar template"""
        # removal may change the first common wavelength; make sure

        instrument_information = dataClass.get_instrument_information()

        epoch_shape = instrument_information["array_size"]
        # Create arrays of zeros in order to open in shared memory and change their values!
        self.spectra = np.zeros(epoch_shape)
        self.rejection_array = np.zeros((len(self.frameIDs_to_use), epoch_shape[0]))

        self.wavelengths = np.zeros(epoch_shape)
        self.uncertainties = np.zeros(epoch_shape)
        self.spectral_mask = None

        try:
            self.package_pool = Queue()
            self.output_pool = Queue()

            self.launch_parallel(dataClass)
            self.evaluate_bad_orders()
            self._finish_template_creation()

        except Exception as e:
            logger.opt(exception=True).critical(
                "Stellar template creation failed due to: {}", e
            )
        finally:
            logger.info("Closing shared memory interfaces of the Stellar template")
            self.cleanup_shared_memory()

    def launch_parallel(self, dataClass):
        inst_info = dataClass.get_instrument_information()
        N_orders = inst_info["array_size"][0]

        epoch_BERVs = dataClass.collect_RV_information(
            KW="BERV",
            subInst=self._associated_subInst,
            frameIDs=self.frameIDs_to_use,
            units=kilometer_second,
            as_value=True,
            include_invalid=False,
        )

        chosen_epochID = self.frameIDs_to_use[np.argmin(epoch_BERVs)]
        self._reference_frameID = chosen_epochID
        self._reference_filepath = dataClass.get_filename_from_frameID(
            self._reference_frameID
        )

        wave_reference, _, _, _ = dataClass.get_frame_arrays_by_ID(chosen_epochID)

        self.wavelengths = remove_RVshift(
            wave_reference,
            stellar_RV=convert_data(
                self.sourceRVs[np.argmin(epoch_BERVs)],
                new_units=kilometer_second,
                as_value=True,
            ),
        )

        logger.info(
            "Using observation from {} as a basis for stellar template construction",
            dataClass.get_frame_by_ID(chosen_epochID),
        )

        logger.info("Using frameIDs: {}", self.frameIDs_to_use)

        kwargs = {
            "valid_epochIDs": self.frameIDs_to_use,
            "chosen_epochID": chosen_epochID,
            "subInst": self._associated_subInst,
            "N_orders": N_orders,
            "dataClass": dataClass,
            "frame_RV_map": {  # construct map between the actual frames and the RVs. Done like this to allow being over-riden easily
                i: j for i, j in zip(self.frameIDs_to_use, self.sourceRVs)
            },
        }

        logger.info("Launching {} workers!", self._internal_configs["NUMBER_WORKERS"])

        # TODO: Avoid error ir we launch this after the template is already in shared memory!
        buffer_sizes = (
            self.wavelengths.shape[0],
            len(self.frameIDs_to_use),
            self.wavelengths.shape[1],
        )
        shr_wave, shr_tmp, shr_uncert, shr_counts = self.convert_to_shared_mem(
            buffer_sizes
        )
        buffers = self.shm

        for _ in range(self._internal_configs["NUMBER_WORKERS"]):
            _ = tqdm.tqdm(
                total=len(self.frameIDs_to_use)
                // self._internal_configs["NUMBER_WORKERS"],
                leave=False,
            )
            p = Process(
                target=self.perform_calculations,
                args=(self.package_pool, self.output_pool, buffers),
                kwargs=kwargs,
            )
            p.start()

        RunTimeRejections = []

        frame_count = 0
        for frameID in self.frameIDs_to_use:
            # to avoid multiple processes opening the arrays at the same time, we open it beforehand
            logger.info("Starting frameID: {}", frameID)
            try:
                _ = dataClass.load_frame_by_ID(frameID)
            except custom_exceptions.FrameError:
                logger.warning("Run Time rejection of frameID = {}", frameID)
                RunTimeRejections.append(frameID)
                continue

            self.used_fpaths.append(
                dataClass.get_filename_from_frameID(frameID, full_path=True)
            )

            total_number_packages = 0
            for order in range(N_orders):
                self.package_pool.put((frameID, order, frame_count))
                total_number_packages += 1
            t = time.time()
            received = 0

            while received != total_number_packages:
                comm_out = self.output_pool.get()
                if not isinstance(comm_out, tuple) and not np.isfinite(comm_out):
                    logger.critical("non finite output")
                    kill_workers(
                        [], self.package_pool, self._internal_configs["NUMBER_WORKERS"]
                    )
                    self._found_error = True
                    raise BadTemplateError("Template creation failed")

                frameID, order, rejection = comm_out
                self.rejection_array[self.frameIDs_to_use.index(frameID), order] = (
                    rejection
                )
                received += 1
            logger.debug(f"Frame took {time.time() - t:0f} seconds")

            if self._internal_configs["MEMORY_SAVE_MODE"]:
                _ = dataClass.close_frame_by_ID(frameID)
            frame_count += 1

        # account for possible Frame Rejections after opening the fits data units
        valid_inds = []
        for frameID in self.frameIDs_to_use:
            if frameID not in RunTimeRejections:
                valid_inds.append(self.frameIDs_to_use.index(frameID))
        self.rejection_array = self.rejection_array[valid_inds]

        for frameID in RunTimeRejections:
            self.frameIDs_to_use.remove(frameID)

        logger.info("Updating template mask")

        self.spectra = np.median(shr_tmp, axis=1)

        new_mask = np.zeros(self.spectra.shape, dtype=bool)

        new_mask[np.where(shr_counts != len(self.frameIDs_to_use))] = True
        new_mask[
            np.where(
                self.spectra < self._internal_configs["FLUX_threshold_for_template"]
            )
        ] = True

        logger.debug("Ensuring increasing wavelenghs in the stellar template")
        # ENsure that we always have increasing wavelengths
        # By construction, the wavelength solution of the template is the "un-masked" wavelengths of the chosen frame
        # When opening the frames we search for non-increasing wavelengths. However, that "mask" is not translated into the
        # actual wavelength solution of the template. The goal of this is to reject those regions (this does not have any
        # impact on the non-affected regions).
        diffs = np.where(np.diff(self.wavelengths, axis=1) < 0)
        if diffs[0].size > 0:
            new_mask[diffs] = True

        self.spectral_mask = Mask(new_mask, mask_type="binary")
        # error propagation for the mean
        self.uncertainties = np.sqrt(np.sum(shr_uncert, axis=1)) / len(
            self.frameIDs_to_use
        )

    def perform_calculations(self, in_queue, out_queue, buffer_info, **kwargs):
        """Compute the stellar template from the input S2D data. Accesses the data from shared memory arrays!"""
        current_subInst = kwargs["subInst"]
        DataClassProxy = kwargs["dataClass"]
        frame_RV_map = kwargs["frame_RV_map"]

        shared_buffers = []
        (
            stellar_template,
            stellar_template_errors,
            stellar_template_wavelengths,
            counts,
            shared_buffers,
        ) = open_buffer(buffer_info, open_type="template", buffers=shared_buffers)

        pixels_in_order = stellar_template[0].size
        try:
            while True:
                data_in = in_queue.get()
                continue_computation = True
                if not isinstance(data_in, (list, tuple)):
                    if not np.isfinite(data_in):
                        return None
                    logger.critical("Wrong data format in the communication queue")
                    raise InvalidConfiguration

                frameID, order, frame_count = data_in

                try:
                    (
                        wavelengths,
                        s2d_data,
                        s2d_uncerts,
                        s2d_mask,
                    ) = DataClassProxy.get_frame_OBS_order(frameID, order)
                except BadOrderError:
                    continue_computation = False

                if continue_computation:
                    current_epochRV = convert_data(
                        frame_RV_map[frameID], new_units=kilometer_second, as_value=True
                    )

                    wavelengths_to_interpolate = np.zeros(
                        stellar_template_wavelengths[order].shape, dtype=bool
                    )

                    # until now the mask has ones in the regions to remove
                    blocks = build_blocks(np.where(~s2d_mask))

                    for block in blocks:
                        start = remove_RVshift(wavelengths[block[0]], current_epochRV)
                        end = remove_RVshift(wavelengths[block[-1]], current_epochRV)
                        interpolation_indexes = np.where(
                            np.logical_and(
                                stellar_template_wavelengths[order] >= start,
                                stellar_template_wavelengths[order] <= end,
                            ),
                        )
                        wavelengths_to_interpolate[interpolation_indexes] = True

                    template_indices = wavelengths_to_interpolate

                    try:
                        interp_ord, interp_err = DataClassProxy.interpolate_frame_order(
                            frameID=frameID,
                            order=order,
                            new_wavelengths=stellar_template_wavelengths[order][
                                template_indices
                            ],
                            shift_RV_by=current_epochRV,
                            RV_shift_mode="remove",
                            include_invalid=False,
                        )

                    except Exception as e:
                        logger.critical("Interpolation failed due to: {}", e)
                        raise e
                    stellar_template[order][frame_count][
                        wavelengths_to_interpolate
                    ] += interp_ord
                    stellar_template_errors[order][frame_count][
                        wavelengths_to_interpolate
                    ] += (interp_err**2)
                    a = counts[order]
                    a[wavelengths_to_interpolate] = a[wavelengths_to_interpolate] + 1
                    counts[order] = a

                    valid_pixels = np.sum(wavelengths_to_interpolate)
                else:
                    valid_pixels = 0
                out_queue.put(
                    (frameID, order, (pixels_in_order - valid_pixels) / pixels_in_order)
                )
        except Exception as e:
            # TODO: fix the procedure for when the workers die

            print(traceback.print_tb(e.__traceback__))
            close_buffers(shared_buffers)

            out_queue.put(np.inf)
            print(f"Template creation dead due to {e}, in order {order}")
            return INTERNAL_ERROR

method_name = 'Median' class-attribute instance-attribute

launch_parallel(dataClass)

Source code in src/ASTRA/template_creation/stellar_templates/median_stellar.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
def launch_parallel(self, dataClass):
    inst_info = dataClass.get_instrument_information()
    N_orders = inst_info["array_size"][0]

    epoch_BERVs = dataClass.collect_RV_information(
        KW="BERV",
        subInst=self._associated_subInst,
        frameIDs=self.frameIDs_to_use,
        units=kilometer_second,
        as_value=True,
        include_invalid=False,
    )

    chosen_epochID = self.frameIDs_to_use[np.argmin(epoch_BERVs)]
    self._reference_frameID = chosen_epochID
    self._reference_filepath = dataClass.get_filename_from_frameID(
        self._reference_frameID
    )

    wave_reference, _, _, _ = dataClass.get_frame_arrays_by_ID(chosen_epochID)

    self.wavelengths = remove_RVshift(
        wave_reference,
        stellar_RV=convert_data(
            self.sourceRVs[np.argmin(epoch_BERVs)],
            new_units=kilometer_second,
            as_value=True,
        ),
    )

    logger.info(
        "Using observation from {} as a basis for stellar template construction",
        dataClass.get_frame_by_ID(chosen_epochID),
    )

    logger.info("Using frameIDs: {}", self.frameIDs_to_use)

    kwargs = {
        "valid_epochIDs": self.frameIDs_to_use,
        "chosen_epochID": chosen_epochID,
        "subInst": self._associated_subInst,
        "N_orders": N_orders,
        "dataClass": dataClass,
        "frame_RV_map": {  # construct map between the actual frames and the RVs. Done like this to allow being over-riden easily
            i: j for i, j in zip(self.frameIDs_to_use, self.sourceRVs)
        },
    }

    logger.info("Launching {} workers!", self._internal_configs["NUMBER_WORKERS"])

    # TODO: Avoid error ir we launch this after the template is already in shared memory!
    buffer_sizes = (
        self.wavelengths.shape[0],
        len(self.frameIDs_to_use),
        self.wavelengths.shape[1],
    )
    shr_wave, shr_tmp, shr_uncert, shr_counts = self.convert_to_shared_mem(
        buffer_sizes
    )
    buffers = self.shm

    for _ in range(self._internal_configs["NUMBER_WORKERS"]):
        _ = tqdm.tqdm(
            total=len(self.frameIDs_to_use)
            // self._internal_configs["NUMBER_WORKERS"],
            leave=False,
        )
        p = Process(
            target=self.perform_calculations,
            args=(self.package_pool, self.output_pool, buffers),
            kwargs=kwargs,
        )
        p.start()

    RunTimeRejections = []

    frame_count = 0
    for frameID in self.frameIDs_to_use:
        # to avoid multiple processes opening the arrays at the same time, we open it beforehand
        logger.info("Starting frameID: {}", frameID)
        try:
            _ = dataClass.load_frame_by_ID(frameID)
        except custom_exceptions.FrameError:
            logger.warning("Run Time rejection of frameID = {}", frameID)
            RunTimeRejections.append(frameID)
            continue

        self.used_fpaths.append(
            dataClass.get_filename_from_frameID(frameID, full_path=True)
        )

        total_number_packages = 0
        for order in range(N_orders):
            self.package_pool.put((frameID, order, frame_count))
            total_number_packages += 1
        t = time.time()
        received = 0

        while received != total_number_packages:
            comm_out = self.output_pool.get()
            if not isinstance(comm_out, tuple) and not np.isfinite(comm_out):
                logger.critical("non finite output")
                kill_workers(
                    [], self.package_pool, self._internal_configs["NUMBER_WORKERS"]
                )
                self._found_error = True
                raise BadTemplateError("Template creation failed")

            frameID, order, rejection = comm_out
            self.rejection_array[self.frameIDs_to_use.index(frameID), order] = (
                rejection
            )
            received += 1
        logger.debug(f"Frame took {time.time() - t:0f} seconds")

        if self._internal_configs["MEMORY_SAVE_MODE"]:
            _ = dataClass.close_frame_by_ID(frameID)
        frame_count += 1

    # account for possible Frame Rejections after opening the fits data units
    valid_inds = []
    for frameID in self.frameIDs_to_use:
        if frameID not in RunTimeRejections:
            valid_inds.append(self.frameIDs_to_use.index(frameID))
    self.rejection_array = self.rejection_array[valid_inds]

    for frameID in RunTimeRejections:
        self.frameIDs_to_use.remove(frameID)

    logger.info("Updating template mask")

    self.spectra = np.median(shr_tmp, axis=1)

    new_mask = np.zeros(self.spectra.shape, dtype=bool)

    new_mask[np.where(shr_counts != len(self.frameIDs_to_use))] = True
    new_mask[
        np.where(
            self.spectra < self._internal_configs["FLUX_threshold_for_template"]
        )
    ] = True

    logger.debug("Ensuring increasing wavelenghs in the stellar template")
    # ENsure that we always have increasing wavelengths
    # By construction, the wavelength solution of the template is the "un-masked" wavelengths of the chosen frame
    # When opening the frames we search for non-increasing wavelengths. However, that "mask" is not translated into the
    # actual wavelength solution of the template. The goal of this is to reject those regions (this does not have any
    # impact on the non-affected regions).
    diffs = np.where(np.diff(self.wavelengths, axis=1) < 0)
    if diffs[0].size > 0:
        new_mask[diffs] = True

    self.spectral_mask = Mask(new_mask, mask_type="binary")
    # error propagation for the mean
    self.uncertainties = np.sqrt(np.sum(shr_uncert, axis=1)) / len(
        self.frameIDs_to_use
    )

perform_calculations(in_queue, out_queue, buffer_info, **kwargs)

Compute the stellar template from the input S2D data. Accesses the data from shared memory arrays!

Source code in src/ASTRA/template_creation/stellar_templates/median_stellar.py 
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
def perform_calculations(self, in_queue, out_queue, buffer_info, **kwargs):
    """Compute the stellar template from the input S2D data. Accesses the data from shared memory arrays!"""
    current_subInst = kwargs["subInst"]
    DataClassProxy = kwargs["dataClass"]
    frame_RV_map = kwargs["frame_RV_map"]

    shared_buffers = []
    (
        stellar_template,
        stellar_template_errors,
        stellar_template_wavelengths,
        counts,
        shared_buffers,
    ) = open_buffer(buffer_info, open_type="template", buffers=shared_buffers)

    pixels_in_order = stellar_template[0].size
    try:
        while True:
            data_in = in_queue.get()
            continue_computation = True
            if not isinstance(data_in, (list, tuple)):
                if not np.isfinite(data_in):
                    return None
                logger.critical("Wrong data format in the communication queue")
                raise InvalidConfiguration

            frameID, order, frame_count = data_in

            try:
                (
                    wavelengths,
                    s2d_data,
                    s2d_uncerts,
                    s2d_mask,
                ) = DataClassProxy.get_frame_OBS_order(frameID, order)
            except BadOrderError:
                continue_computation = False

            if continue_computation:
                current_epochRV = convert_data(
                    frame_RV_map[frameID], new_units=kilometer_second, as_value=True
                )

                wavelengths_to_interpolate = np.zeros(
                    stellar_template_wavelengths[order].shape, dtype=bool
                )

                # until now the mask has ones in the regions to remove
                blocks = build_blocks(np.where(~s2d_mask))

                for block in blocks:
                    start = remove_RVshift(wavelengths[block[0]], current_epochRV)
                    end = remove_RVshift(wavelengths[block[-1]], current_epochRV)
                    interpolation_indexes = np.where(
                        np.logical_and(
                            stellar_template_wavelengths[order] >= start,
                            stellar_template_wavelengths[order] <= end,
                        ),
                    )
                    wavelengths_to_interpolate[interpolation_indexes] = True

                template_indices = wavelengths_to_interpolate

                try:
                    interp_ord, interp_err = DataClassProxy.interpolate_frame_order(
                        frameID=frameID,
                        order=order,
                        new_wavelengths=stellar_template_wavelengths[order][
                            template_indices
                        ],
                        shift_RV_by=current_epochRV,
                        RV_shift_mode="remove",
                        include_invalid=False,
                    )

                except Exception as e:
                    logger.critical("Interpolation failed due to: {}", e)
                    raise e
                stellar_template[order][frame_count][
                    wavelengths_to_interpolate
                ] += interp_ord
                stellar_template_errors[order][frame_count][
                    wavelengths_to_interpolate
                ] += (interp_err**2)
                a = counts[order]
                a[wavelengths_to_interpolate] = a[wavelengths_to_interpolate] + 1
                counts[order] = a

                valid_pixels = np.sum(wavelengths_to_interpolate)
            else:
                valid_pixels = 0
            out_queue.put(
                (frameID, order, (pixels_in_order - valid_pixels) / pixels_in_order)
            )
    except Exception as e:
        # TODO: fix the procedure for when the workers die

        print(traceback.print_tb(e.__traceback__))
        close_buffers(shared_buffers)

        out_queue.put(np.inf)
        print(f"Template creation dead due to {e}, in order {order}")
        return INTERNAL_ERROR

sum_stellar

SumStellar

Bases: StellarTemplate

Typical stellar template, constructed by shifting all observations to a common rest frame.

(i.e. by removing their stellar RVs) and computing a mean of all observations. Uncertainties can be propagated either through an analytical formula or by interpolating the original uncertainties of each observation.

Note: Check the User parameters of the parent classes for further customization options of this class

Source code in src/ASTRA/template_creation/stellar_templates/sum_stellar.py 
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
class SumStellar(StellarTemplate):
    """Typical stellar template, constructed by shifting all observations to a common rest frame.

    (i.e. by removing their stellar RVs) and computing a mean of all observations. Uncertainties can be
    propagated either through an analytical formula or by interpolating the original uncertainties of each
    observation.

    *Note:* Check the **User parameters** of the parent classes for further customization options of this class

    """

    method_name = choices.STELLAR_CREATION_MODE.Sum.value

    _default_params = StellarTemplate._default_params + DefaultValues(
        ALIGNEMENT_RV_SOURCE=UserParam(
            "DRS", constraint=ValueFromIterable(["DRS", "SBART"])
        ),
        FLUX_threshold_for_template=UserParam(
            default_value=1,
            constraint=Positive_Value_Constraint,
            description="Flux threshold for masking the spectral template. Set to one to avoid possible numerical"
            " issues with near-zero values",
        ),
        ENSURE_COMMON_WAVELENGTHS=UserParam(
            default_value=True,
            constraint=BooleanValue,
            description="If True (default) only keep pixels that are informed by all spectra in use. If False, "
            "use everything available. This will mean that the final SNR will not have a consistent SNR across"
            "wavelengths",
        ),
    )

    def __init__(
        self, subInst: str, user_configs: Optional[Dict] = None, loaded: bool = False
    ):
        super().__init__(subInst=subInst, user_configs=user_configs, loaded=loaded)

        if not loaded:
            if self._internal_configs["MEMORY_SAVE_MODE"]:
                logger.warning(
                    "Stellar template creation will save RAM usage. This will result in multiple open/close "
                    "operations across the entire SBART pipeline! ",
                )

            self._found_error = False

    def _run_template_construction(
        self,
        dataClass: DataClass,
        conditions,
        reference_frame,
    ) -> None:
        """Create the stellar template."""
        # removal may change the first common wavelength; make sure

        instrument_information = dataClass.get_instrument_information()

        epoch_shape = instrument_information["array_size"]

        if self._internal_configs["OVERSAMPLE_TEMPLATE"] > 1:
            logger.warning("Oversampling template wavelengths")
            epoch_shape = (
                epoch_shape[0],
                epoch_shape[1] * self._internal_configs["OVERSAMPLE_TEMPLATE"],
            )

        # Create arrays of zeros in order to open in shared memory and change their values!
        self.spectra = np.zeros(epoch_shape)
        self.rejection_array = np.zeros((len(self.frameIDs_to_use), epoch_shape[0]))

        self.wavelengths = np.zeros(epoch_shape)
        self.uncertainties = np.zeros(epoch_shape)
        self.spectral_mask = None

        try:
            self.package_pool = Queue()
            self.output_pool = Queue()

            self.launch_parallel(dataClass, reference_frame=reference_frame)
            self.evaluate_bad_orders()
            self._finish_template_creation()

        except Exception as e:
            logger.opt(exception=True).critical(
                "Stellar template creation failed due to: {}", e
            )
        finally:
            logger.info("Closing shared memory interfaces of the Stellar template")
            self.cleanup_shared_memory()

    def launch_parallel(self, dataClass: DataClass, reference_frame: None | Frame):
        """Launch parallel computation of the stellar template.

        Raises:
            BadTemplateError: The worker dies and computation is halted.

        """
        inst_info = dataClass.get_instrument_information()
        N_orders = inst_info["array_size"][0]

        if reference_frame is not None:
            logger.warning(
                "Using reference frame as the basis for wavelength grid of template"
            )
            chosen_epochID = reference_frame.frameID
            self._reference_filepath = reference_frame.fname
            self._reference_frameID = chosen_epochID

            wave_reference, _, _, _ = reference_frame.get_data_from_full_spectrum()

            self.wavelengths = remove_RVshift(
                wave_reference,
                stellar_RV=convert_data(
                    reference_frame.get_KW_value(self.RV_keyword),
                    new_units=kilometer_second,
                    as_value=True,
                ),
            )

        else:
            epoch_BERVs = dataClass.collect_RV_information(
                KW="BERV",
                subInst=self._associated_subInst,
                frameIDs=self.frameIDs_to_use,
                units=kilometer_second,
                as_value=True,
                include_invalid=False,
            )

            chosen_epochID = self.frameIDs_to_use[np.argmin(epoch_BERVs)]
            self._reference_frameID = chosen_epochID
            self._reference_filepath = dataClass.get_filename_from_frameID(
                self._reference_frameID
            )

            wave_reference, _, _, _ = dataClass.get_frame_arrays_by_ID(chosen_epochID)

            wavelengths = remove_RVshift(
                wave_reference,
                stellar_RV=convert_data(
                    self.sourceRVs[np.argmin(epoch_BERVs)],
                    new_units=kilometer_second,
                    as_value=True,
                ),
            )

            logger.info(
                "Using observation from {} as a basis for stellar template construction",
                dataClass.get_frame_by_ID(chosen_epochID),
            )

        for order_index, wave_order in enumerate(wavelengths):
            if self._internal_configs["OVERSAMPLE_TEMPLATE"] > 1:
                new_wave = np.linspace(
                    wave_order[0],
                    wave_order[-1],
                    wavelengths.shape[1]
                    * self._internal_configs["OVERSAMPLE_TEMPLATE"],
                )
                self.wavelengths[order_index] = new_wave
            else:
                self.wavelengths[order_index] = wavelengths[order_index]

        logger.info("Using frameIDs: {}", self.frameIDs_to_use)
        kwargs = {
            "valid_epochIDs": self.frameIDs_to_use,
            "chosen_epochID": chosen_epochID,
            "subInst": self._associated_subInst,
            "N_orders": N_orders,
            "dataClass": dataClass,
            # construct map between the actual frames and the RVs. Done like this to allow being over-riden easily
            "frame_RV_map": dict(zip(self.frameIDs_to_use, self.sourceRVs)),
        }

        logger.info("Launching {} workers!", self._internal_configs["NUMBER_WORKERS"])

        # TODO: Avoid error ir we launch this after the template is already in shared memory!
        shr_wave, shr_tmp, shr_uncert, shr_counts = self.convert_to_shared_mem()
        buffers = self.shm

        for _ in range(self._internal_configs["NUMBER_WORKERS"]):
            p = Process(
                target=self.perform_calculations,
                args=(self.package_pool, self.output_pool, buffers),
                kwargs=kwargs,
            )
            p.start()

        RunTimeRejections = []

        start_time = time.time()
        for frameID in self.frameIDs_to_use:
            # to avoid multiple processes opening the arrays at the same time, we open it beforehand
            try:
                _ = dataClass.load_frame_by_ID(frameID)
            except custom_exceptions.FrameError:
                logger.warning("Run Time rejection of frameID = {}", frameID)
                RunTimeRejections.append(frameID)
                continue

            self.used_fpaths.append(
                dataClass.get_filename_from_frameID(frameID, full_path=True)
            )

            total_number_packages = 0
            for order in range(N_orders):
                self.package_pool.put((frameID, order))
                total_number_packages += 1
            received = 0

            while received != total_number_packages:
                comm_out = self.output_pool.get()
                if not isinstance(comm_out, tuple) and not np.isfinite(comm_out):
                    logger.critical("non finite output")
                    kill_workers(
                        [], self.package_pool, self._internal_configs["NUMBER_WORKERS"]
                    )
                    self._found_error = True
                    msg = "Template creation failed"
                    raise BadTemplateError(msg)

                frameID, order, rejection = comm_out
                self.rejection_array[self.frameIDs_to_use.index(frameID), order] = (
                    rejection
                )
                received += 1

            if self._internal_configs["MEMORY_SAVE_MODE"]:
                _ = dataClass.close_frame_by_ID(frameID)

        # account for possible Frame Rejections after opening the fits data units
        valid_inds = []
        for frameID in self.frameIDs_to_use:
            if frameID not in RunTimeRejections:
                valid_inds.append(self.frameIDs_to_use.index(frameID))
        self.rejection_array = self.rejection_array[valid_inds]

        for frameID in RunTimeRejections:
            self.frameIDs_to_use.remove(frameID)

        logger.info("Updating template mask")

        self.spectra = shr_tmp[:] / len(self.frameIDs_to_use)

        new_mask = np.zeros(self.spectra.shape, dtype=bool)

        if self._internal_configs["ENSURE_COMMON_WAVELENGTHS"]:
            new_mask[np.where(shr_counts != len(self.frameIDs_to_use))] = True
        new_mask[
            np.where(
                self.spectra < self._internal_configs["FLUX_threshold_for_template"]
            )
        ] = True

        logger.debug("Ensuring increasing wavelenghs in the stellar template")
        # ENsure that we always have increasing wavelengths
        # By construction, the wavelength solution of the template is the "un-masked" wavelengths of the chosen frame
        # When opening the frames we search for non-increasing wavelengths. However, that "mask" is not translated into the
        # actual wavelength solution of the template. The goal of this is to reject those regions (this does not have any
        # impact on the non-affected regions).
        diffs = np.where(np.diff(self.wavelengths, axis=1) < 0)
        if diffs[0].size > 0:
            new_mask[diffs] = True

        self.spectral_mask = Mask(new_mask, mask_type="binary")
        # error propagation for the mean
        self.uncertainties = np.sqrt(shr_uncert[:]) / len(self.frameIDs_to_use)
        logger.info(f"Construction of stellar template took {time.time() - start_time}")

    def add_new_frame_to_template(self, frame: Frame):
        """Add an individual frame to template."""
        super().add_new_frame_to_template(frame)

        if not frame.is_valid:
            logger.critical(
                "Injected frame does not pass the SBART-defined QC flags, template will not be updated"
            )
            return

        self.spectra *= len(self.used_fpaths)
        self.uncertainties *= len(self.used_fpaths)
        self.uncertainties **= 2
        self.used_fpaths.append(frame.file_path)

        for order in range(frame.N_orders):
            if order in self.bad_orders:
                # no need to handle the orders that are marked as "bad" in the original template
                continue

            try:
                (
                    wavelengths,
                    s2d_data,
                    s2d_uncerts,
                    s2d_mask,
                ) = frame.get_data_from_spectral_order(order)
            except BadOrderError:
                self.spectral_mask.add_indexes_to_mask_order(
                    order=order,
                    indexes=np.ones(self.template[order], dtype=bool),
                    mask_type=MISSING_DATA,
                )
                logger.critical(
                    "New frame has an invalid order, completely rejecting spectral order from template"
                )
                continue

            current_epochRV = convert_data(
                frame.get_KW_value("DRS_RV"), new_units=kilometer_second, as_value=True
            )

            wavelengths_to_interpolate = np.zeros(wavelengths.shape, dtype=bool)

            # until now the mask has ones in the regions to remove
            blocks = build_blocks(np.where(~s2d_mask))

            for block in blocks:
                start = remove_RVshift(wavelengths[block[0]], current_epochRV)
                end = remove_RVshift(wavelengths[block[-1]], current_epochRV)
                interpolation_indexes = np.where(
                    np.logical_and(
                        self.wavelengths[order] >= start,
                        self.wavelengths[order] <= end,
                    ),
                )
                wavelengths_to_interpolate[interpolation_indexes] = True

            int_ord, int_err = frame.interpolate_spectrum_to_wavelength(
                order=order,
                new_wavelengths=self.wavelengths[order][wavelengths_to_interpolate],
                shift_RV_by=current_epochRV,
                RV_shift_mode="remove",
                include_invalid=False,
            )

            self.spectra[order][wavelengths_to_interpolate] += int_ord
            self.uncertainties[order][wavelengths_to_interpolate] += int_err**2

            # reject non-interpolated points
            self.spectral_mask.add_indexes_to_mask_order(
                order=order,
                indexes=np.where(wavelengths_to_interpolate == False),
                mask_type=MISSING_DATA,
            )

            # reject low-flux points
            self.spectral_mask.add_indexes_to_mask_order(
                order=order,
                indexes=np.where(
                    self.spectra[order]
                    < self._internal_configs["FLUX_threshold_for_template"]
                ),
                mask_type=MISSING_DATA,
            )

        self.uncertainties = np.sqrt(self.uncertainties) / len(self.used_fpaths)
        self.spectra /= len(self.used_fpaths)

    def perform_calculations(self, in_queue, out_queue, buffer_info, **kwargs):
        """Compute the stellar template from the input S2D data.

        Accesses the data from shared memory arrays!
        """
        DataClassProxy: DataClass = kwargs["dataClass"]
        frame_RV_map = kwargs["frame_RV_map"]

        shared_buffers = []
        (
            stellar_template,
            stellar_template_errors,
            stellar_template_wavelengths,
            counts,
            shared_buffers,
        ) = open_buffer(buffer_info, open_type="template", buffers=shared_buffers)

        pixels_in_order = stellar_template[0].size
        try:
            while True:
                data_in = in_queue.get()
                continue_computation = True
                if not isinstance(data_in, (list, tuple)):
                    if not np.isfinite(data_in):
                        return None
                    logger.critical("Wrong data format in the communication queue")
                    raise InvalidConfiguration

                frameID, order = data_in

                try:
                    (
                        wavelengths,
                        s2d_data,
                        s2d_uncerts,
                        s2d_mask,
                    ) = DataClassProxy.get_frame_OBS_order(frameID, order)
                except BadOrderError:
                    continue_computation = False

                if continue_computation:
                    current_epochRV = convert_data(
                        frame_RV_map[frameID], new_units=kilometer_second, as_value=True
                    )

                    wavelengths_to_interpolate = np.zeros(
                        stellar_template_wavelengths[order].shape, dtype=bool
                    )

                    # until now the mask has ones in the regions to remove
                    blocks = build_blocks(np.where(~s2d_mask))

                    for block in blocks:
                        start = remove_RVshift(wavelengths[block[0]], current_epochRV)
                        end = remove_RVshift(wavelengths[block[-1]], current_epochRV)
                        interpolation_indexes = np.where(
                            np.logical_and(
                                stellar_template_wavelengths[order] >= start,
                                stellar_template_wavelengths[order] <= end,
                            ),
                        )
                        wavelengths_to_interpolate[interpolation_indexes] = True

                    template_indices = wavelengths_to_interpolate

                    try:
                        interp_ord, interp_err = DataClassProxy.interpolate_frame_order(
                            frameID=frameID,
                            order=order,
                            new_wavelengths=stellar_template_wavelengths[order][
                                template_indices
                            ],
                            shift_RV_by=current_epochRV,
                            RV_shift_mode="remove",
                            include_invalid=False,
                        )

                    except Exception as e:
                        logger.critical("Interpolation failed due to: {}", e)
                        raise e

                    stellar_template[order][wavelengths_to_interpolate] += interp_ord
                    stellar_template_errors[order][wavelengths_to_interpolate] += (
                        interp_err**2
                    )
                    a = counts[order]
                    a[wavelengths_to_interpolate] = a[wavelengths_to_interpolate] + 1
                    counts[order] = a

                    valid_pixels = np.sum(wavelengths_to_interpolate)
                else:
                    valid_pixels = 0
                out_queue.put(
                    (frameID, order, (pixels_in_order - valid_pixels) / pixels_in_order)
                )
        except Exception as e:
            # TODO: fix the procedure for when the workers die

            print(traceback.print_tb(e.__traceback__))
            close_buffers(shared_buffers)

            out_queue.put(np.inf)
            print(f"Template creation dead due to {e}, in order {order}")
            return INTERNAL_ERROR

method_name = choices.STELLAR_CREATION_MODE.Sum.value class-attribute instance-attribute

launch_parallel(dataClass: DataClass, reference_frame: None | Frame)

Launch parallel computation of the stellar template.

Raises: BadTemplateError: The worker dies and computation is halted.

Source code in src/ASTRA/template_creation/stellar_templates/sum_stellar.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
def launch_parallel(self, dataClass: DataClass, reference_frame: None | Frame):
    """Launch parallel computation of the stellar template.

    Raises:
        BadTemplateError: The worker dies and computation is halted.

    """
    inst_info = dataClass.get_instrument_information()
    N_orders = inst_info["array_size"][0]

    if reference_frame is not None:
        logger.warning(
            "Using reference frame as the basis for wavelength grid of template"
        )
        chosen_epochID = reference_frame.frameID
        self._reference_filepath = reference_frame.fname
        self._reference_frameID = chosen_epochID

        wave_reference, _, _, _ = reference_frame.get_data_from_full_spectrum()

        self.wavelengths = remove_RVshift(
            wave_reference,
            stellar_RV=convert_data(
                reference_frame.get_KW_value(self.RV_keyword),
                new_units=kilometer_second,
                as_value=True,
            ),
        )

    else:
        epoch_BERVs = dataClass.collect_RV_information(
            KW="BERV",
            subInst=self._associated_subInst,
            frameIDs=self.frameIDs_to_use,
            units=kilometer_second,
            as_value=True,
            include_invalid=False,
        )

        chosen_epochID = self.frameIDs_to_use[np.argmin(epoch_BERVs)]
        self._reference_frameID = chosen_epochID
        self._reference_filepath = dataClass.get_filename_from_frameID(
            self._reference_frameID
        )

        wave_reference, _, _, _ = dataClass.get_frame_arrays_by_ID(chosen_epochID)

        wavelengths = remove_RVshift(
            wave_reference,
            stellar_RV=convert_data(
                self.sourceRVs[np.argmin(epoch_BERVs)],
                new_units=kilometer_second,
                as_value=True,
            ),
        )

        logger.info(
            "Using observation from {} as a basis for stellar template construction",
            dataClass.get_frame_by_ID(chosen_epochID),
        )

    for order_index, wave_order in enumerate(wavelengths):
        if self._internal_configs["OVERSAMPLE_TEMPLATE"] > 1:
            new_wave = np.linspace(
                wave_order[0],
                wave_order[-1],
                wavelengths.shape[1]
                * self._internal_configs["OVERSAMPLE_TEMPLATE"],
            )
            self.wavelengths[order_index] = new_wave
        else:
            self.wavelengths[order_index] = wavelengths[order_index]

    logger.info("Using frameIDs: {}", self.frameIDs_to_use)
    kwargs = {
        "valid_epochIDs": self.frameIDs_to_use,
        "chosen_epochID": chosen_epochID,
        "subInst": self._associated_subInst,
        "N_orders": N_orders,
        "dataClass": dataClass,
        # construct map between the actual frames and the RVs. Done like this to allow being over-riden easily
        "frame_RV_map": dict(zip(self.frameIDs_to_use, self.sourceRVs)),
    }

    logger.info("Launching {} workers!", self._internal_configs["NUMBER_WORKERS"])

    # TODO: Avoid error ir we launch this after the template is already in shared memory!
    shr_wave, shr_tmp, shr_uncert, shr_counts = self.convert_to_shared_mem()
    buffers = self.shm

    for _ in range(self._internal_configs["NUMBER_WORKERS"]):
        p = Process(
            target=self.perform_calculations,
            args=(self.package_pool, self.output_pool, buffers),
            kwargs=kwargs,
        )
        p.start()

    RunTimeRejections = []

    start_time = time.time()
    for frameID in self.frameIDs_to_use:
        # to avoid multiple processes opening the arrays at the same time, we open it beforehand
        try:
            _ = dataClass.load_frame_by_ID(frameID)
        except custom_exceptions.FrameError:
            logger.warning("Run Time rejection of frameID = {}", frameID)
            RunTimeRejections.append(frameID)
            continue

        self.used_fpaths.append(
            dataClass.get_filename_from_frameID(frameID, full_path=True)
        )

        total_number_packages = 0
        for order in range(N_orders):
            self.package_pool.put((frameID, order))
            total_number_packages += 1
        received = 0

        while received != total_number_packages:
            comm_out = self.output_pool.get()
            if not isinstance(comm_out, tuple) and not np.isfinite(comm_out):
                logger.critical("non finite output")
                kill_workers(
                    [], self.package_pool, self._internal_configs["NUMBER_WORKERS"]
                )
                self._found_error = True
                msg = "Template creation failed"
                raise BadTemplateError(msg)

            frameID, order, rejection = comm_out
            self.rejection_array[self.frameIDs_to_use.index(frameID), order] = (
                rejection
            )
            received += 1

        if self._internal_configs["MEMORY_SAVE_MODE"]:
            _ = dataClass.close_frame_by_ID(frameID)

    # account for possible Frame Rejections after opening the fits data units
    valid_inds = []
    for frameID in self.frameIDs_to_use:
        if frameID not in RunTimeRejections:
            valid_inds.append(self.frameIDs_to_use.index(frameID))
    self.rejection_array = self.rejection_array[valid_inds]

    for frameID in RunTimeRejections:
        self.frameIDs_to_use.remove(frameID)

    logger.info("Updating template mask")

    self.spectra = shr_tmp[:] / len(self.frameIDs_to_use)

    new_mask = np.zeros(self.spectra.shape, dtype=bool)

    if self._internal_configs["ENSURE_COMMON_WAVELENGTHS"]:
        new_mask[np.where(shr_counts != len(self.frameIDs_to_use))] = True
    new_mask[
        np.where(
            self.spectra < self._internal_configs["FLUX_threshold_for_template"]
        )
    ] = True

    logger.debug("Ensuring increasing wavelenghs in the stellar template")
    # ENsure that we always have increasing wavelengths
    # By construction, the wavelength solution of the template is the "un-masked" wavelengths of the chosen frame
    # When opening the frames we search for non-increasing wavelengths. However, that "mask" is not translated into the
    # actual wavelength solution of the template. The goal of this is to reject those regions (this does not have any
    # impact on the non-affected regions).
    diffs = np.where(np.diff(self.wavelengths, axis=1) < 0)
    if diffs[0].size > 0:
        new_mask[diffs] = True

    self.spectral_mask = Mask(new_mask, mask_type="binary")
    # error propagation for the mean
    self.uncertainties = np.sqrt(shr_uncert[:]) / len(self.frameIDs_to_use)
    logger.info(f"Construction of stellar template took {time.time() - start_time}")

add_new_frame_to_template(frame: Frame)

Add an individual frame to template.

Source code in src/ASTRA/template_creation/stellar_templates/sum_stellar.py 
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
def add_new_frame_to_template(self, frame: Frame):
    """Add an individual frame to template."""
    super().add_new_frame_to_template(frame)

    if not frame.is_valid:
        logger.critical(
            "Injected frame does not pass the SBART-defined QC flags, template will not be updated"
        )
        return

    self.spectra *= len(self.used_fpaths)
    self.uncertainties *= len(self.used_fpaths)
    self.uncertainties **= 2
    self.used_fpaths.append(frame.file_path)

    for order in range(frame.N_orders):
        if order in self.bad_orders:
            # no need to handle the orders that are marked as "bad" in the original template
            continue

        try:
            (
                wavelengths,
                s2d_data,
                s2d_uncerts,
                s2d_mask,
            ) = frame.get_data_from_spectral_order(order)
        except BadOrderError:
            self.spectral_mask.add_indexes_to_mask_order(
                order=order,
                indexes=np.ones(self.template[order], dtype=bool),
                mask_type=MISSING_DATA,
            )
            logger.critical(
                "New frame has an invalid order, completely rejecting spectral order from template"
            )
            continue

        current_epochRV = convert_data(
            frame.get_KW_value("DRS_RV"), new_units=kilometer_second, as_value=True
        )

        wavelengths_to_interpolate = np.zeros(wavelengths.shape, dtype=bool)

        # until now the mask has ones in the regions to remove
        blocks = build_blocks(np.where(~s2d_mask))

        for block in blocks:
            start = remove_RVshift(wavelengths[block[0]], current_epochRV)
            end = remove_RVshift(wavelengths[block[-1]], current_epochRV)
            interpolation_indexes = np.where(
                np.logical_and(
                    self.wavelengths[order] >= start,
                    self.wavelengths[order] <= end,
                ),
            )
            wavelengths_to_interpolate[interpolation_indexes] = True

        int_ord, int_err = frame.interpolate_spectrum_to_wavelength(
            order=order,
            new_wavelengths=self.wavelengths[order][wavelengths_to_interpolate],
            shift_RV_by=current_epochRV,
            RV_shift_mode="remove",
            include_invalid=False,
        )

        self.spectra[order][wavelengths_to_interpolate] += int_ord
        self.uncertainties[order][wavelengths_to_interpolate] += int_err**2

        # reject non-interpolated points
        self.spectral_mask.add_indexes_to_mask_order(
            order=order,
            indexes=np.where(wavelengths_to_interpolate == False),
            mask_type=MISSING_DATA,
        )

        # reject low-flux points
        self.spectral_mask.add_indexes_to_mask_order(
            order=order,
            indexes=np.where(
                self.spectra[order]
                < self._internal_configs["FLUX_threshold_for_template"]
            ),
            mask_type=MISSING_DATA,
        )

    self.uncertainties = np.sqrt(self.uncertainties) / len(self.used_fpaths)
    self.spectra /= len(self.used_fpaths)

perform_calculations(in_queue, out_queue, buffer_info, **kwargs)

Compute the stellar template from the input S2D data.

Accesses the data from shared memory arrays!

Source code in src/ASTRA/template_creation/stellar_templates/sum_stellar.py 
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
def perform_calculations(self, in_queue, out_queue, buffer_info, **kwargs):
    """Compute the stellar template from the input S2D data.

    Accesses the data from shared memory arrays!
    """
    DataClassProxy: DataClass = kwargs["dataClass"]
    frame_RV_map = kwargs["frame_RV_map"]

    shared_buffers = []
    (
        stellar_template,
        stellar_template_errors,
        stellar_template_wavelengths,
        counts,
        shared_buffers,
    ) = open_buffer(buffer_info, open_type="template", buffers=shared_buffers)

    pixels_in_order = stellar_template[0].size
    try:
        while True:
            data_in = in_queue.get()
            continue_computation = True
            if not isinstance(data_in, (list, tuple)):
                if not np.isfinite(data_in):
                    return None
                logger.critical("Wrong data format in the communication queue")
                raise InvalidConfiguration

            frameID, order = data_in

            try:
                (
                    wavelengths,
                    s2d_data,
                    s2d_uncerts,
                    s2d_mask,
                ) = DataClassProxy.get_frame_OBS_order(frameID, order)
            except BadOrderError:
                continue_computation = False

            if continue_computation:
                current_epochRV = convert_data(
                    frame_RV_map[frameID], new_units=kilometer_second, as_value=True
                )

                wavelengths_to_interpolate = np.zeros(
                    stellar_template_wavelengths[order].shape, dtype=bool
                )

                # until now the mask has ones in the regions to remove
                blocks = build_blocks(np.where(~s2d_mask))

                for block in blocks:
                    start = remove_RVshift(wavelengths[block[0]], current_epochRV)
                    end = remove_RVshift(wavelengths[block[-1]], current_epochRV)
                    interpolation_indexes = np.where(
                        np.logical_and(
                            stellar_template_wavelengths[order] >= start,
                            stellar_template_wavelengths[order] <= end,
                        ),
                    )
                    wavelengths_to_interpolate[interpolation_indexes] = True

                template_indices = wavelengths_to_interpolate

                try:
                    interp_ord, interp_err = DataClassProxy.interpolate_frame_order(
                        frameID=frameID,
                        order=order,
                        new_wavelengths=stellar_template_wavelengths[order][
                            template_indices
                        ],
                        shift_RV_by=current_epochRV,
                        RV_shift_mode="remove",
                        include_invalid=False,
                    )

                except Exception as e:
                    logger.critical("Interpolation failed due to: {}", e)
                    raise e

                stellar_template[order][wavelengths_to_interpolate] += interp_ord
                stellar_template_errors[order][wavelengths_to_interpolate] += (
                    interp_err**2
                )
                a = counts[order]
                a[wavelengths_to_interpolate] = a[wavelengths_to_interpolate] + 1
                counts[order] = a

                valid_pixels = np.sum(wavelengths_to_interpolate)
            else:
                valid_pixels = 0
            out_queue.put(
                (frameID, order, (pixels_in_order - valid_pixels) / pixels_in_order)
            )
    except Exception as e:
        # TODO: fix the procedure for when the workers die

        print(traceback.print_tb(e.__traceback__))
        close_buffers(shared_buffers)

        out_queue.put(np.inf)
        print(f"Template creation dead due to {e}, in order {order}")
        return INTERNAL_ERROR

Telluric templates

Handle Telluric models.

TelluricModel

Bases: TemplateFramework

The TelluricModel is responsible for the creation of the telluric template for each sub-Instrument.

This object supports the following user parameters:

User parameters:

================================ ================ ================ Parameter name Mandatory Default Value ================================ ================ ================ CREATION_MODE True ---- APPLICATION_MODE False removal EXTENSION_MODE False lines ================================ ================ ================

Note: Also check the User parameters of the parent classes for further customization options of SBART

Source code in src/ASTRA/template_creation/TelluricModel.py 
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
class TelluricModel(TemplateFramework):
    # noinspection LongLine
    """The TelluricModel is responsible for the creation of the telluric template for each sub-Instrument.

    This object supports the following user parameters:

    **User parameters:**

    ================================ ================ ================
    Parameter name                      Mandatory      Default Value
    ================================ ================ ================
    CREATION_MODE                       True            ----
    APPLICATION_MODE                    False           removal
    EXTENSION_MODE                      False           lines
    ================================ ================ ================


    *Note:* Also check the **User parameters** of the parent classes for further customization options of SBART

    """

    model_type = "Telluric"

    template_map: ClassVar[dict[TELLURIC_CREATION_MODE, TelluricTemplate]] = {
        TELLURIC_CREATION_MODE.telfit: TelfitTelluric,
        TELLURIC_CREATION_MODE.OHemission: OHemissionTelluric,
    }

    _default_params = TemplateFramework._default_params + DefaultValues(
        CREATION_MODE=UserParam(
            TELLURIC_CREATION_MODE.telfit,
            constraint=ValueFromIterable(TELLURIC_CREATION_MODE),
            mandatory=False,
        ),
        APPLICATION_MODE=UserParam(
            TELLURIC_APPLICATION_MODE.removal,
            constraint=ValueFromIterable(TELLURIC_APPLICATION_MODE),
        ),
        EXTENSION_MODE=UserParam(
            TELLURIC_EXTENSION.LINES,
            constraint=ValueFromIterable(TELLURIC_EXTENSION),
        ),
    )

    def __init__(self, usage_mode: str, root_folder_path: UI_PATH, user_configs: UI_DICT) -> None:
        """Instantiate of the object.

        Parameters
        ----------
        usage_mode : str
            How to telluric templates will be applied to the data. If 'individual' -> template applied to its own
            subInstrument. If 'merged', we use a template built from merging the templates from all subInstruments
        root_folder_path: Union [str, pathlib.Path]
            Path to the folder inside which SBART will store its outputs
        user_configs: Optional[Dict[str, Any]]
            Dictionary with the keys and values of the user parameters that have been described above

        """
        super().__init__(mode="", root_folder_path=root_folder_path, user_configs=user_configs)

        logger.info("Starting Telluric Model")

        # can be 'individual' -> template applied to its own subInstrument
        # or  'merged' and we use a template built from merging the templates from all subInstruments
        self._usage_mode = usage_mode

    def load_templates_from_disk(self) -> None:  # noqa: D102
        super().load_templates_from_disk()
        for value in self.templates.values():
            value.update_extension_mode(self._internal_configs["EXTENSION_MODE"])

    def request_template(self, subInstrument: str) -> Type[BaseTemplate]:
        """Return the template for a given sub-Instrument.

        Parameters
        ----------
        subInstrument: str
            "name" of the subInstrument

        Returns
        -------
            Requested telluric Template

        """
        logger.debug("Serving {} template to subInstrument {}", self._usage_mode, subInstrument)
        if self._usage_mode == "":
            return self.templates["merged"]
        if self._usage_mode == "individual":
            return self.templates[subInstrument]
        raise InvalidConfiguration()

    def Generate_Model(
        self,
        dataClass: DataClass,
        telluric_configs: dict,
        force_computation: bool = False,
        store_templates: bool = True,
    ) -> None:
        """Generate a telluric model for all subInstruments with data, as defined in the parent implementation.

        Afterwards, allow to combine the telluric model of all sub-Instruments into a single telluric binary
         model which will then be used for all available observations.

        [**Warning**: the combination is yet to be implemented]

        Parameters
        ----------
        dataClass : :class:`~ASTRAdata_objects.DataClass`
            DataClass with the observations
        telluric_configs : dict
            Dictionary with the parameters needed to control the creation of the telluric template, following the
            specifications of the templates.
        force_computation : bool
            If False, it will attempt to lead the telluric template from disk before trying to compute them. If True,
            always compute telluric template, even if it exists
        store_templates: bool
            If True (default) store to disk

        """
        super().Generate_Model(
            dataClass=dataClass,
            template_configs=telluric_configs,
            attempt_to_load=not force_computation,
            store_templates=False,
        )

        if self._usage_mode == "merged":
            logger.info("Telluric model merging the templates from all epochs")
            self.merge_templates()

        if store_templates:
            self.store_templates_to_disk()

    def merge_templates(self) -> None:
        """Merge the telluric template of all sub-Instruments to create a master telluric binary template.

        Raises
        ------
        NotImplementedError
            The method is yet to be implemented

        """
        logger.info("Merging templates to create a global one")
        raise NotImplementedError

    # Internal Usage:

    def _find_templates_from_disk(self, which: TELLURIC_CREATION_MODE) -> list[str]:
        loading_path = self._internalPaths.get_path_to(self.__class__.model_type)
        logger.info("Loading {} template from disk inside directory", self.__class__.model_type)
        logger.info("\t" + loading_path)

        available_templates = []

        if not os.path.exists(loading_path):
            logger.warning(f"Could not find template to load in {loading_path}")
            raise TemplateNotExistsError()

        for fname in os.listdir(loading_path):
            if which.value in fname and fname.endswith("fits"):
                available_templates.append(fname)
        logger.info(
            "Found {} available templates: {} of type {}",
            len(available_templates),
            available_templates,
            which,
        )
        if len(available_templates) == 0:
            logger.critical("Could not find templates to load!")
            raise TemplateNotExistsError()
        return [os.path.join(loading_path, i) for i in available_templates]

    def _compute_template(self, data, subInstrument: str, user_configs: dict) -> TelluricTemplate:
        creation_mode = self._internal_configs["CREATION_MODE"]
        logger.info("Using template of type: {}", creation_mode)

        if creation_mode == "none":
            pass
        elif creation_mode == TELLURIC_CREATION_MODE.telfit:
            chosen_template = TelfitTelluric
        elif creation_mode == TELLURIC_CREATION_MODE.tapas:
            chosen_template = TapasTelluric
        elif creation_mode == TELLURIC_CREATION_MODE.OHemission:
            chosen_template = OHemissionTelluric
        else:
            raise InvalidConfiguration()

        if self.work_mode == WORKING_MODE.ONE_SHOT or self.templates[subInstrument] is None:
            tell_template = chosen_template(
                subInst=subInstrument,
                user_configs=user_configs,
                extension_mode=self._internal_configs["EXTENSION_MODE"],
                application_mode=self._internal_configs["APPLICATION_MODE"],
            )

            tell_template.load_information_from_DataClass(data)

            if self.is_for_removal:
                tell_template.create_telluric_template(dataClass=data)
            else:
                logger.debug("Telluric template in removal mode. Fitting from inside dataClass")

        else:
            # If this is a rolling mode, we just need to add the new information
            tell_template = self.templates[subInstrument]
            tell_template.ingest_new_rolling_observations(dataClass=data)

            if not self.is_for_removal:
                msg = "The rolling template can only be used to remove features, not correct "
                raise InvalidConfiguration(msg)

        return tell_template

    @property
    def is_for_removal(self) -> bool:
        """True if the template will be used to remove telluric features."""
        return self._internal_configs["APPLICATION_MODE"] == TELLURIC_APPLICATION_MODE.removal

    @property
    def is_for_correction(self) -> bool:
        """True if the template will be used to correct telluric features."""
        return self._internal_configs["APPLICATION_MODE"] == TELLURIC_APPLICATION_MODE.correction

model_type = 'Telluric' class-attribute instance-attribute

template_map: dict[TELLURIC_CREATION_MODE, TelluricTemplate] = {TELLURIC_CREATION_MODE.telfit: TelfitTelluric, TELLURIC_CREATION_MODE.OHemission: OHemissionTelluric} class-attribute

is_for_removal: bool property

True if the template will be used to remove telluric features.

is_for_correction: bool property

True if the template will be used to correct telluric features.

load_templates_from_disk() -> None

Source code in src/ASTRA/template_creation/TelluricModel.py 
96
97
98
99
def load_templates_from_disk(self) -> None:  # noqa: D102
    super().load_templates_from_disk()
    for value in self.templates.values():
        value.update_extension_mode(self._internal_configs["EXTENSION_MODE"])

request_template(subInstrument: str) -> Type[BaseTemplate]

Return the template for a given sub-Instrument.

Parameters:

Name Type Description Default
subInstrument str

"name" of the subInstrument

 required

Returns:

Type Description
Requested telluric Template
Source code in src/ASTRA/template_creation/TelluricModel.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def request_template(self, subInstrument: str) -> Type[BaseTemplate]:
    """Return the template for a given sub-Instrument.

    Parameters
    ----------
    subInstrument: str
        "name" of the subInstrument

    Returns
    -------
        Requested telluric Template

    """
    logger.debug("Serving {} template to subInstrument {}", self._usage_mode, subInstrument)
    if self._usage_mode == "":
        return self.templates["merged"]
    if self._usage_mode == "individual":
        return self.templates[subInstrument]
    raise InvalidConfiguration()

Generate_Model(dataClass: DataClass, telluric_configs: dict, force_computation: bool = False, store_templates: bool = True) -> None

Generate a telluric model for all subInstruments with data, as defined in the parent implementation.

Afterwards, allow to combine the telluric model of all sub-Instruments into a single telluric binary model which will then be used for all available observations.

[Warning: the combination is yet to be implemented]

Parameters:

Name Type Description Default
dataClass :class:`~ASTRAdata_objects.DataClass`

DataClass with the observations

 required
telluric_configs dict

Dictionary with the parameters needed to control the creation of the telluric template, following the specifications of the templates.

 required
force_computation bool

If False, it will attempt to lead the telluric template from disk before trying to compute them. If True, always compute telluric template, even if it exists

 False
store_templates bool

If True (default) store to disk

 True
Source code in src/ASTRA/template_creation/TelluricModel.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def Generate_Model(
    self,
    dataClass: DataClass,
    telluric_configs: dict,
    force_computation: bool = False,
    store_templates: bool = True,
) -> None:
    """Generate a telluric model for all subInstruments with data, as defined in the parent implementation.

    Afterwards, allow to combine the telluric model of all sub-Instruments into a single telluric binary
     model which will then be used for all available observations.

    [**Warning**: the combination is yet to be implemented]

    Parameters
    ----------
    dataClass : :class:`~ASTRAdata_objects.DataClass`
        DataClass with the observations
    telluric_configs : dict
        Dictionary with the parameters needed to control the creation of the telluric template, following the
        specifications of the templates.
    force_computation : bool
        If False, it will attempt to lead the telluric template from disk before trying to compute them. If True,
        always compute telluric template, even if it exists
    store_templates: bool
        If True (default) store to disk

    """
    super().Generate_Model(
        dataClass=dataClass,
        template_configs=telluric_configs,
        attempt_to_load=not force_computation,
        store_templates=False,
    )

    if self._usage_mode == "merged":
        logger.info("Telluric model merging the templates from all epochs")
        self.merge_templates()

    if store_templates:
        self.store_templates_to_disk()

merge_templates() -> None

Merge the telluric template of all sub-Instruments to create a master telluric binary template.

Raises:

Type Description
NotImplementedError

The method is yet to be implemented
Source code in src/ASTRA/template_creation/TelluricModel.py 
163
164
165
166
167
168
169
170
171
172
173
def merge_templates(self) -> None:
    """Merge the telluric template of all sub-Instruments to create a master telluric binary template.

    Raises
    ------
    NotImplementedError
        The method is yet to be implemented

    """
    logger.info("Merging templates to create a global one")
    raise NotImplementedError

Different models for the creation of telluric masks from synthetic spectra of Earth's transmittance profile.

Telluric_Template

TelluricTemplate

Bases: BaseTemplate

BaseClass of a telluric template, with all the necessary functionalities.

Inherits from the BaseTemplate, similarly to the stellar templates. The telluric template does not have the same shape as the S2D spectra. As we only need the wavelength regions in which it has detected a feature, we can preserve the original shape of the transmittance spectra.

User parameters:

================================ ================ ================ ================ ================ Parameter name Mandatory Default Value Valid Values Comment ================================ ================ ================ ================ ================ continuum_percentage_drop False 1 Inside [0, 100] [1] force_download False False boolean [2] ================================ ================ ================ ================ ================

  • [1]: Minimum drop (in relation to continuum level) needed to flag a point as a telluric feature
  • [2]: Force the download of all necessary data products, even if not necessary

Note: Also check the User parameters of the parent classes for further customization options of SBART

Subclassed by:
Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
class TelluricTemplate(BaseTemplate):
    """BaseClass of a telluric template, with all the necessary functionalities.

    Inherits from the BaseTemplate, similarly to the stellar templates. The telluric template does not
    have the same shape as the S2D spectra. As we only need the wavelength regions in which it has
    detected a feature, we can preserve the original shape of the transmittance spectra.

    **User parameters:**

    ================================ ================ ================ ================ ================
    Parameter name                      Mandatory      Default Value    Valid Values    Comment
    ================================ ================ ================ ================ ================
    continuum_percentage_drop           False              1            Inside [0, 100]  [1]
    force_download                      False              False         boolean         [2]
    ================================ ================ ================ ================ ================

    - [1]: Minimum drop (in relation to continuum level) needed to flag a point as a telluric feature
    - [2]: Force the download of all necessary data products, even if not necessary

    *Note:* Also check the **User parameters** of the parent classes for further customization options of SBART

    """

    _name = "Telluric"
    _default_params = BaseTemplate._default_params + DefaultValues(
        continuum_percentage_drop=UserParam(
            1, constraint=ValueInInterval([0, 100], include_edges=True)
        ),
        force_download=UserParam(False, constraint=BooleanValue),
        inverse_mask=UserParam(
            False,
            constraint=BooleanValue,
            description="If True, rejects everything that is not flagged as a telluric. Default is False.",
        ),
    )
    template_type = "Telluric"
    method_name = "Base"

    def __init__(
        self,
        subInst: str,
        extension_mode: TELLURIC_EXTENSION,
        user_configs: Union[None, dict] = None,
        application_mode: TELLURIC_APPLICATION_MODE = TELLURIC_APPLICATION_MODE.removal,
        loaded: bool = False,
    ):
        """Construct base telluric template.

        Parameters
        ----------
        extension_mode : str, optional
            How to handle with Earth motion during the year. If 'lines', compute the template in each possible location
            for the available observations and only remove the. If "window", each telluric feature will be extended by a
            region equal to the maximum BERV. By default: 'lines'

        """
        super().__init__(subInst, user_configs, loaded)

        self._associated_BERV = 0

        self._application_mode = application_mode
        self._extension_mode = extension_mode
        self._base_mask: list[list[float]] = []

        self._masked_wavelengths: list[list[float, float]] = []
        self._computed_wave_blocks = False

        self.transmittance_wavelengths = None
        self.transmittance_spectra = None
        self._continuum_level = None

        self._fitModel = Model(params_of_model=[])
        self.BERVS: list[float] = []
        self.MAXBERV = None
        self._reference_frameID = None  # used for telluric removal!

        if self.for_feature_correction:
            # Not the best way of doing this, but ....
            self.add_to_status(SUCCESS)

        self._loaded_dataclass_info = False

        self.template: np.ndarray = None

        self._metric_selection_conditions = Empty_condition()

        self.use_approximated_BERV_correction = False
        if self._internal_configs["inverse_mask"]:
            if self._extension_mode != TELLURIC_EXTENSION.WINDOW:
                raise custom_exceptions.InternalError(
                    "Can't inverse the telluric mask without using a window extension mode",
                )

    # Workaround to avoid calling spectra to the binary telluric template (to make things internally consistent)
    @property
    def spectra(self):
        return self.template

    @spectra.setter
    def spectra(self, updated_value):
        self.internal_val = updated_value

    def get_data_from_spectral_order(self, order: int, include_invalid: bool = False):
        uncertainties = (
            self.uncertainties[order] if self.for_feature_correction else None
        )
        return self.wavelengths[order], self.template[order], uncertainties

    def get_data_from_full_spectrum(self):
        uncertainties = self.uncertainties if self.for_feature_correction else None
        return self.wavelengths, self.template, uncertainties

    #####

    def _generate_model_parameters(self, dataClass):
        logger.debug("Generating the parameters of the model applied in {}", self.name)
        self._fitModel.generate_priors(dataClass)

    def load_information_from_DataClass(self, dataClass: DataClass) -> None:
        """Load the necessary information from the dataClass to generate the models.

        Loads:
            - BERV information (for the telluric removal)
            - Generates the model information (i.e. populates priors and bounds of self._fitModel

        Parameters
        ----------
        dataClass


        """
        logger.debug("Template {} loading data from the dataClass", self.name)

        self._load_BERV_info(dataClass)

        frames = dataClass.get_valid_frameIDS()
        zeroth_frame = dataClass.get_frame_by_ID(frames[0])

        self.use_approximated_BERV_correction = (
            zeroth_frame.use_approximated_BERV_correction
        )
        if not self.is_valid:
            return

        self._generate_model_parameters(dataClass)
        self._loaded_dataclass_info = True

    # Internal Data loading routines
    def _load_BERV_info(self, DataClass) -> None:
        """Load the BERV and MaxBERV of all observations from the subInstrument associated with this template.

        Parameters
        ----------
        DataClass : [type]
            [description]

        Returns
        -------
        Tuple[List[float], List[float]]
            [description]

        """
        if not self.is_valid:
            self.MAXBERV = np.nan * kilometer_second
            return

        if (
            self._associated_subInst
            not in DataClass.get_subInstruments_with_valid_frames()
        ):
            logger.warning(
                "{} has no valid observations. Not computing telluric template",
                self._associated_subInst,
            )
            self.add_to_status(MISSING_DATA)
            self.MAXBERV = np.nan * kilometer_second
            return

        BERVS = DataClass.collect_KW_observations(
            KW="BERV",
            subInstruments=[self._associated_subInst],
            include_invalid=False,
        )
        max_bervs = DataClass.collect_KW_observations(
            KW="MAX_BERV",
            subInstruments=[self._associated_subInst],
            include_invalid=False,
        )

        # It seems that numpy does not like lists of astropy.units elements
        unitless_max_bervs = [i.value for i in max_bervs]

        # if self.work_mode == WORKING_MODE.ROLLING:
        #     # No need to repeat BERVS
        #     self.BERVS.extend(BERVS)
        #     self.BERVS = list(set(self.BERVS))
        #     new_berv_max = max_bervs[np.argmax(unitless_max_bervs)]
        #     self.MAXBERV = max(new_berv_max, self.MAXBERV)
        # else:
        self.BERVS.extend(BERVS)
        self.MAXBERV = max_bervs[np.argmax(np.abs(unitless_max_bervs))]

    def _search_reference_frame(self, dataclass: DataClass) -> Union[int, float]:
        """Select the frame that will be used to construct the telluric template.

        By default, select the one with the highest relative humidity. If that
        keyword is not loaded, then uses the one with the highest airmass. If
        there are no valid frames in the associated subINstrument, returns
        np.nan

        Parameters
        ----------
        dataclass : DataClass
            [description]

        Returns
        -------
        int
            [description]

        """
        try:
            valid_frame_ids = dataclass.get_frameIDs_from_subInst(
                self._associated_subInst
            )
        except NoDataError:
            msg = f"{self._associated_subInst} has no valid observations. Not computing telluric template"

            logger.warning(msg)
            self.add_to_status(MISSING_DATA(msg))
            return np.nan

        # Placing upper limit of temperature at 50ºC
        self._metric_selection_conditions += KEYWORD_condition(
            "ambient_temperature", [[None, 323.15]]
        )

        if self.__class__.method_name.lower() == "telfit":
            # 1 December 2014, because no GDAS profile for telfit
            # ! may blow up if it removes all observations :(
            # add condition so that the reference observation is more than week
            # ago, to guarantee the GDAS profile already exists
            one_week_ago = int(Time.now().jd - 7)
            self._metric_selection_conditions += KEYWORD_condition(
                "BJD", [[2453340, one_week_ago]]
            )

        logger.debug(
            "Using Relative humidity as the selection criterion for reference observation"
        )
        metric_to_select = dataclass.collect_KW_observations(
            KW="relative_humidity",
            subInstruments=[self._associated_subInst],
            include_invalid=False,
            conditions=self._metric_selection_conditions,
        )

        # due to the conditions, some elements may be None
        metric_to_select = [-1 if m is None else m for m in metric_to_select]

        if not np.isfinite(metric_to_select[0]):
            logger.warning(
                "Relative humidity keyword was not loaded. Using airmass to select the reference observation for {}",
                self._associated_subInst,
            )
            metric_to_select = dataclass.collect_KW_observations(
                KW="airmass",
                subInstruments=[self._associated_subInst],
                include_invalid=False,
                conditions=self._metric_selection_conditions,
            )

            # due to the conditions, some elements may be None
            metric_to_select = [-1 if m is None else m for m in metric_to_select]

        chosen_frameID = valid_frame_ids[np.argmax(metric_to_select)]
        self._associated_BERV = dataclass.get_frame_by_ID(chosen_frameID).get_KW_value(
            "BERV"
        )

        logger.info(
            "Telluric Template from {} using {} as the reference",
            self._associated_subInst,
            dataclass.get_frame_by_ID(chosen_frameID),
        )

        return chosen_frameID

    #########################################
    #  Creation of the telluric binary mask #
    #########################################

    def create_telluric_template(
        self, dataClass, custom_frameID: Optional[int] = None
    ) -> None:
        logger.info(
            "Starting creation of {} template from {}",
            self.__class__.template_type,
            self._associated_subInst,
        )

        if custom_frameID is not None:
            logger.info(
                "Creation of Telluric Template is based on a custom frameID: {}",
                custom_frameID,
            )
            self._reference_frameID = custom_frameID
        else:
            self._reference_frameID = self._search_reference_frame(dataClass)

        if self.for_feature_correction:
            raise custom_exceptions.InvalidConfiguration(
                "No need to generate binary mask for a correction model"
            )

        self._base_checks_for_template_creation()

        if not self._loaded_dataclass_info:
            raise custom_exceptions.InvalidConfiguration(
                f"{self.name} did not load dataClass information"
            )

        logger.info(
            "Starting telluric template creation, with reference ID = {}",
            self._reference_frameID,
        )

        self._fitModel.disable_full_model()

    def update_extension_mode(self, new_mode: TELLURIC_EXTENSION) -> None:
        """Update the telluric extension mode.

        Args:
            new_mode (TELLURIC_EXTENSION): New mode.

        Raises:
            custom_exceptions.InvalidConfiguration: Mode is not recognized

        """
        if new_mode not in TELLURIC_EXTENSION:
            msg = f"{new_mode} not recognized"
            raise custom_exceptions.InvalidConfiguration(msg)

        self._computed_wave_blocks = False
        self._masked_wavelengths = []
        self._extension_mode = new_mode

    def ingest_new_rolling_observations(self, dataClass):
        """Ingest new observations into a rolling mode template."""
        if self.work_mode != WORKING_MODE.ROLLING:
            msg = "Can't access ROLLING mode functions without being in it"
            raise custom_exceptions.InvalidConfiguration(msg)

        self._computed_wave_blocks = False
        self._load_BERV_info(DataClass=dataClass)
        self._masked_wavelengths = []

    ###
    #   transmittance post-processing
    ###
    def _compute_wave_blocks(self) -> None:
        """Go from the binary template to a list of wavelength regions in which the template is non-zero.

        If the template was not computed, add nothing and raise a warning.
        """
        if not self.is_valid:
            msg = "Creating list of wavelength region (with tellurics) from a non-computed Telluric template. No telluric features identified"
            logger.warning(msg)
            raise custom_exceptions.InternalError(msg)

        # TODO:  optimize the list of blocked features! It will have a very bad scaling with N_{obs}
        logger.info("Creating list of blocked features due to tellurics")
        if not self.was_loaded:
            logger.info(
                "Extending telluric features with the mode: <{}>", self._extension_mode
            )

        if len(self._base_mask) == 0:
            self.build_blocks()

        updated_block = []

        for interval in self._base_mask:
            # first overlap search, to take advantage of the smaller list size in here (when compared against the "global" one)
            # The intervals are the start and end of each block
            updated_block.extend(
                find_overlaps(self._extend_detections([interval[0], interval[-1]]))
            )

        self._masked_wavelengths = find_overlaps(updated_block)
        self._computed_wave_blocks = True

        if self._internal_configs["inverse_mask"]:
            new_blocks = [[0, self._masked_wavelengths[0][0]]]
            for index in range(len(self._masked_wavelengths) - 1):
                new_blocks.append(
                    (
                        self._masked_wavelengths[index][1],
                        self._masked_wavelengths[index + 1][0],
                    ),
                )
            new_blocks.append(
                (
                    self._masked_wavelengths[-1][1],
                    self._masked_wavelengths[-1][1] * 1000,
                ),
            )

            self._masked_wavelengths = new_blocks

    def _extend_detections(self, telluric_block: list[list]) -> list[list[float]]:
        """Extend each block of telluric detection based on the self._extension_mode that was selected by the user.

        Independently of the extension mode, the telluric features are widened by +- 15m/s

        Args:
            telluric_block (list[list]): Start and end of a wavelength block that has a telluric feature!

        Returns:
            list[list[float]]: Updated position of the feature

        """
        updated_block = []

        berv_function = (
            apply_approximated_BERV_correction
            if self.use_approximated_BERV_correction
            else apply_BERV_correction
        )

        if self._extension_mode == TELLURIC_EXTENSION.LINES:
            for berv_value in self.BERVS:
                berv = berv_value.to(kilometer_second).value
                offset = 15 / 1000  # expand each line by 15 m/s to account for wings

                lowest_wavelength = berv_function(telluric_block[0], BERV=berv - offset)
                highest_wavelength = berv_function(
                    telluric_block[1], BERV=berv + offset
                )
                updated_block.append([lowest_wavelength, highest_wavelength])

        elif self._extension_mode == TELLURIC_EXTENSION.WINDOW:
            berv = self.MAXBERV.to(kilometer_second).value + 15 / 1000
            lowest_wavelength = berv_function(telluric_block[0], BERV=-berv)
            highest_wavelength = berv_function(telluric_block[1], BERV=berv)
            updated_block.append([lowest_wavelength, highest_wavelength])
        else:
            raise custom_exceptions.InternalError(
                f"{self._extension_mode} not recognized"
            )
        return updated_block

    def create_binary_template(self, continuum_level) -> None:
        """Transform transmittance spectra to binary template, based on percentage deviations from the continuum level.

        Parameters
        ----------
        continuum_level : np.ndarray()
            Continuum level that is estimated by the children classes

        """
        logger.info("Converting from transmittance spectra to binary mask!")

        if not self.for_feature_removal:
            logger.warning(
                "Telluric Template will not be used to remove spectra. No need to create binary mask"
            )
            return

        telluric_mask = create_binary_template(
            transmittance=self.transmittance_spectra,
            continuum_level=continuum_level,
            percentage_drop=self._internal_configs["continuum_percentage_drop"],
        )
        # Find a decrease of 1% in relation to the continuum level; Positive
        # gains (against the continuum value) are not considered as tellurics
        self.template = telluric_mask

        self.build_blocks()
        self._compute_wave_blocks()

    def build_blocks(self) -> None:
        indexes = build_blocks(np.where(self.template != 0))
        for telluric_block in indexes:
            wavel = self.wavelengths[telluric_block]
            self._base_mask.append([wavel[0], wavel[-1]])

    #######################################
    #  Outside access to the properties   #
    #######################################

    @property
    def contaminated_regions(self) -> list:
        """List of contaminated regions."""
        if not self.for_feature_removal:
            raise custom_exceptions.InvalidConfiguration(
                "{} is not a template constructed for telluric removal!"
            )
        if not self._computed_wave_blocks:
            logger.debug("No previous computation of wavelength blocks. Doing it now!")
            self._compute_wave_blocks()
        return self._masked_wavelengths

    @property
    def storage_name(self) -> str:
        return f"{self.__class__.method_name}_{self.__class__.template_type}"

    @property
    def for_feature_removal(self) -> bool:
        return self._application_mode == TELLURIC_APPLICATION_MODE.removal

    @property
    def for_feature_correction(self) -> bool:
        return self._application_mode == TELLURIC_APPLICATION_MODE.correction

    def store_metrics(self): ...

    def trigger_data_storage(self, *args, **kwargs) -> NoReturn:
        try:
            super().trigger_data_storage(*args, **kwargs)
            self.store_metrics()
        except custom_exceptions.FailedStorage:
            return

    #################################
    #           Data storage        #
    #################################

    def store_template(self, clobber: bool):
        if self.for_feature_correction:
            logger.info(
                "{} is on correction mode. Not storing conventional data products to disk",
                self.name,
            )
            return

        try:
            super().store_template(clobber)
        except custom_exceptions.FailedStorage:
            return

        header = fits.Header()

        header["TYPE"] = self.__class__.template_type
        header["subInst"] = self._associated_subInst
        header["VERSION"] = __version__
        header["IS_VALID"] = self.is_valid
        header["HIERARCH APPROX BERV CORRECTION"] = (
            self.use_approximated_BERV_correction
        )
        header["MAX_BERV"] = convert_data(
            self.MAXBERV, new_units=kilometer_second, as_value=True
        )

        for key, config_val in self._internal_configs.items():
            if (
                "path" in key.lower()
                or "user_" in key
                or isinstance(config_val, (list, tuple))
            ):
                continue

            if "FIT" in key:
                continue

            if key in ["SAVE_DISK_SPACE"]:
                continue

            if key in ["WORKING_MODE"]:
                header[f"HIERARCH {key}"] = config_val.value
            else:
                header[f"HIERARCH {key}"] = config_val
        hdu = fits.PrimaryHDU(data=[], header=header)

        contam = self._base_mask
        contam_imge = np.zeros((len(contam), 2))
        for row, entry in enumerate(contam):
            contam_imge[row] = entry
        hdu_contam = fits.ImageHDU(data=contam_imge, header=header, name="CONTAM")

        hdu_berv = fits.ImageHDU(
            data=convert_data(self.BERVS, new_units=kilometer_second, as_value=True),
            header=header,
            name="BERVS",
        )

        hdus_cubes = [hdu, hdu_contam, hdu_berv]

        if self.disk_save_level != DISK_SAVE_MODE.EXTREME:
            hdu_wave = fits.ImageHDU(data=self.wavelengths, header=header, name="Wave")
            complete_template = np.zeros(self.template.shape)
            for pair in self._masked_wavelengths:
                indexes = np.where(
                    np.logical_and(
                        self.wavelengths >= pair[0], self.wavelengths <= pair[1]
                    )
                )
                complete_template[indexes] = 1

            hdu_temp = fits.ImageHDU(data=complete_template, header=header, name="Temp")

            for val in [hdu_wave, hdu_temp]:
                hdus_cubes.append(val)

            hdu_transWave = fits.ImageHDU(
                data=self.transmittance_wavelengths, header=header, name="TRANSMIT_WAVE"
            )
            hdu_transSpec = fits.ImageHDU(
                data=self.transmittance_spectra, header=header, name="TRANSMIT_SPECTRA"
            )
            hdus_cubes.extend([hdu_transWave, hdu_transSpec])

        hdul = fits.HDUList(hdus_cubes)

        filename = f"{self.storage_name}_{self._associated_subInst}.fits"
        logger.debug(
            "Storing template to {}", self._internalPaths.root_storage_path / filename
        )
        hdul.writeto(self._internalPaths.root_storage_path / filename, overwrite=True)
        if self.disk_save_level != DISK_SAVE_MODE.EXTREME:
            metrics_path = self._internalPaths.get_path_to("metrics", as_posix=False)
            fig, axis = plt.subplots()
            axis.plot(self.transmittance_wavelengths, self.transmittance_spectra)
            axis.set_xlabel(r"Wavelength [$\AA$]")
            axis.set_ylabel("Transmittance")
            fig.savefig(metrics_path / f"transmittance_{self._associated_subInst}.png")
            plt.close(fig)

    def load_from_file(self, root_path: Path, loading_path: str) -> None:
        """TODO: save and load the actual flag to disk!

        Parameters
        ----------
        root_path
        loading_path

        Returns
        -------

        """
        super().load_from_file(root_path, loading_path)

        with fits.open(loading_path) as hdulist:
            if hdulist[1].header.get("VERSION", "") != __version__:
                logger.warning(
                    "Loaded template was not created under the current SBART version. Possible problems may arise",
                )

            self._associated_subInst = hdulist["CONTAM"].header["subInst"]
            self.MAXBERV = hdulist["CONTAM"].header["MAX_BERV"] * kilometer_second
            self.BERVS = [i * kilometer_second for i in hdulist["BERVS"].data.tolist()]

            # By default we load with a window extension mode
            self._extension_mode = TELLURIC_EXTENSION.WINDOW

            try:
                waves = hdulist["Wave"].data
                template = hdulist["Temp"].data
            except:
                waves, template = None, None
            try:
                self.transmittance_wavelengths = hdulist["TRANSMIT_WAVE"].data
                self.transmittance_spectra = hdulist["TRANSMIT_SPECTRA"].data
            except:
                self.transmittance_spectra = None
                self.transmittance_wavelengths = None

            try:
                self.template = template
                self.wavelengths = waves
            except KeyError:
                self.template = None
                self.wavelengths = None

            try:
                self.use_approximated_BERV_correction = hdulist[1].header[
                    "HIERARCH APPROX BERV CORRECTION"
                ]
            except KeyError:
                logger.warning("Loading old telluric template with missing keywords")
                self.use_approximated_BERV_correction = False

            self._computed_wave_blocks = False
            self._base_mask = hdulist["CONTAM"].data.tolist()
        self.add_to_status(DISK_LOADED_DATA(f"Loaded data from {loading_path}"))

    def _finish_template_creation(self):
        if self.template is None:
            # This ensures that we don't create multiple times the binary mask
            self.create_binary_template(self._continuum_level)

        super()._finish_template_creation()

template_type = 'Telluric' class-attribute instance-attribute

method_name = 'Base' class-attribute instance-attribute

transmittance_wavelengths = None instance-attribute

transmittance_spectra = None instance-attribute

BERVS: list[float] = [] instance-attribute

MAXBERV = None instance-attribute

template: np.ndarray = None instance-attribute

use_approximated_BERV_correction = False instance-attribute

spectra property writable

contaminated_regions: list property

List of contaminated regions.

storage_name: str property

for_feature_removal: bool property

for_feature_correction: bool property

get_data_from_spectral_order(order: int, include_invalid: bool = False)

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
146
147
148
149
150
def get_data_from_spectral_order(self, order: int, include_invalid: bool = False):
    uncertainties = (
        self.uncertainties[order] if self.for_feature_correction else None
    )
    return self.wavelengths[order], self.template[order], uncertainties

get_data_from_full_spectrum()

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
152
153
154
def get_data_from_full_spectrum(self):
    uncertainties = self.uncertainties if self.for_feature_correction else None
    return self.wavelengths, self.template, uncertainties

load_information_from_DataClass(dataClass: DataClass) -> None

Load the necessary information from the dataClass to generate the models.

Loads: - BERV information (for the telluric removal) - Generates the model information (i.e. populates priors and bounds of self._fitModel

Parameters:

Name Type Description Default
dataClass DataClass
required
Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def load_information_from_DataClass(self, dataClass: DataClass) -> None:
    """Load the necessary information from the dataClass to generate the models.

    Loads:
        - BERV information (for the telluric removal)
        - Generates the model information (i.e. populates priors and bounds of self._fitModel

    Parameters
    ----------
    dataClass


    """
    logger.debug("Template {} loading data from the dataClass", self.name)

    self._load_BERV_info(dataClass)

    frames = dataClass.get_valid_frameIDS()
    zeroth_frame = dataClass.get_frame_by_ID(frames[0])

    self.use_approximated_BERV_correction = (
        zeroth_frame.use_approximated_BERV_correction
    )
    if not self.is_valid:
        return

    self._generate_model_parameters(dataClass)
    self._loaded_dataclass_info = True

create_telluric_template(dataClass, custom_frameID: Optional[int] = None) -> None

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
def create_telluric_template(
    self, dataClass, custom_frameID: Optional[int] = None
) -> None:
    logger.info(
        "Starting creation of {} template from {}",
        self.__class__.template_type,
        self._associated_subInst,
    )

    if custom_frameID is not None:
        logger.info(
            "Creation of Telluric Template is based on a custom frameID: {}",
            custom_frameID,
        )
        self._reference_frameID = custom_frameID
    else:
        self._reference_frameID = self._search_reference_frame(dataClass)

    if self.for_feature_correction:
        raise custom_exceptions.InvalidConfiguration(
            "No need to generate binary mask for a correction model"
        )

    self._base_checks_for_template_creation()

    if not self._loaded_dataclass_info:
        raise custom_exceptions.InvalidConfiguration(
            f"{self.name} did not load dataClass information"
        )

    logger.info(
        "Starting telluric template creation, with reference ID = {}",
        self._reference_frameID,
    )

    self._fitModel.disable_full_model()

update_extension_mode(new_mode: TELLURIC_EXTENSION) -> None

Update the telluric extension mode.

Args: new_mode (TELLURIC_EXTENSION): New mode.

Raises: custom_exceptions.InvalidConfiguration: Mode is not recognized

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
def update_extension_mode(self, new_mode: TELLURIC_EXTENSION) -> None:
    """Update the telluric extension mode.

    Args:
        new_mode (TELLURIC_EXTENSION): New mode.

    Raises:
        custom_exceptions.InvalidConfiguration: Mode is not recognized

    """
    if new_mode not in TELLURIC_EXTENSION:
        msg = f"{new_mode} not recognized"
        raise custom_exceptions.InvalidConfiguration(msg)

    self._computed_wave_blocks = False
    self._masked_wavelengths = []
    self._extension_mode = new_mode

ingest_new_rolling_observations(dataClass)

Ingest new observations into a rolling mode template.

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
391
392
393
394
395
396
397
398
399
def ingest_new_rolling_observations(self, dataClass):
    """Ingest new observations into a rolling mode template."""
    if self.work_mode != WORKING_MODE.ROLLING:
        msg = "Can't access ROLLING mode functions without being in it"
        raise custom_exceptions.InvalidConfiguration(msg)

    self._computed_wave_blocks = False
    self._load_BERV_info(DataClass=dataClass)
    self._masked_wavelengths = []

create_binary_template(continuum_level) -> None

Transform transmittance spectra to binary template, based on percentage deviations from the continuum level.

Parameters:

Name Type Description Default
continuum_level ndarray()

Continuum level that is estimated by the children classes

 required
Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
def create_binary_template(self, continuum_level) -> None:
    """Transform transmittance spectra to binary template, based on percentage deviations from the continuum level.

    Parameters
    ----------
    continuum_level : np.ndarray()
        Continuum level that is estimated by the children classes

    """
    logger.info("Converting from transmittance spectra to binary mask!")

    if not self.for_feature_removal:
        logger.warning(
            "Telluric Template will not be used to remove spectra. No need to create binary mask"
        )
        return

    telluric_mask = create_binary_template(
        transmittance=self.transmittance_spectra,
        continuum_level=continuum_level,
        percentage_drop=self._internal_configs["continuum_percentage_drop"],
    )
    # Find a decrease of 1% in relation to the continuum level; Positive
    # gains (against the continuum value) are not considered as tellurics
    self.template = telluric_mask

    self.build_blocks()
    self._compute_wave_blocks()

build_blocks() -> None

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
525
526
527
528
529
def build_blocks(self) -> None:
    indexes = build_blocks(np.where(self.template != 0))
    for telluric_block in indexes:
        wavel = self.wavelengths[telluric_block]
        self._base_mask.append([wavel[0], wavel[-1]])

store_metrics()

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
559
def store_metrics(self): ...

trigger_data_storage(*args, **kwargs) -> NoReturn

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
561
562
563
564
565
566
def trigger_data_storage(self, *args, **kwargs) -> NoReturn:
    try:
        super().trigger_data_storage(*args, **kwargs)
        self.store_metrics()
    except custom_exceptions.FailedStorage:
        return

store_template(clobber: bool)

Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
def store_template(self, clobber: bool):
    if self.for_feature_correction:
        logger.info(
            "{} is on correction mode. Not storing conventional data products to disk",
            self.name,
        )
        return

    try:
        super().store_template(clobber)
    except custom_exceptions.FailedStorage:
        return

    header = fits.Header()

    header["TYPE"] = self.__class__.template_type
    header["subInst"] = self._associated_subInst
    header["VERSION"] = __version__
    header["IS_VALID"] = self.is_valid
    header["HIERARCH APPROX BERV CORRECTION"] = (
        self.use_approximated_BERV_correction
    )
    header["MAX_BERV"] = convert_data(
        self.MAXBERV, new_units=kilometer_second, as_value=True
    )

    for key, config_val in self._internal_configs.items():
        if (
            "path" in key.lower()
            or "user_" in key
            or isinstance(config_val, (list, tuple))
        ):
            continue

        if "FIT" in key:
            continue

        if key in ["SAVE_DISK_SPACE"]:
            continue

        if key in ["WORKING_MODE"]:
            header[f"HIERARCH {key}"] = config_val.value
        else:
            header[f"HIERARCH {key}"] = config_val
    hdu = fits.PrimaryHDU(data=[], header=header)

    contam = self._base_mask
    contam_imge = np.zeros((len(contam), 2))
    for row, entry in enumerate(contam):
        contam_imge[row] = entry
    hdu_contam = fits.ImageHDU(data=contam_imge, header=header, name="CONTAM")

    hdu_berv = fits.ImageHDU(
        data=convert_data(self.BERVS, new_units=kilometer_second, as_value=True),
        header=header,
        name="BERVS",
    )

    hdus_cubes = [hdu, hdu_contam, hdu_berv]

    if self.disk_save_level != DISK_SAVE_MODE.EXTREME:
        hdu_wave = fits.ImageHDU(data=self.wavelengths, header=header, name="Wave")
        complete_template = np.zeros(self.template.shape)
        for pair in self._masked_wavelengths:
            indexes = np.where(
                np.logical_and(
                    self.wavelengths >= pair[0], self.wavelengths <= pair[1]
                )
            )
            complete_template[indexes] = 1

        hdu_temp = fits.ImageHDU(data=complete_template, header=header, name="Temp")

        for val in [hdu_wave, hdu_temp]:
            hdus_cubes.append(val)

        hdu_transWave = fits.ImageHDU(
            data=self.transmittance_wavelengths, header=header, name="TRANSMIT_WAVE"
        )
        hdu_transSpec = fits.ImageHDU(
            data=self.transmittance_spectra, header=header, name="TRANSMIT_SPECTRA"
        )
        hdus_cubes.extend([hdu_transWave, hdu_transSpec])

    hdul = fits.HDUList(hdus_cubes)

    filename = f"{self.storage_name}_{self._associated_subInst}.fits"
    logger.debug(
        "Storing template to {}", self._internalPaths.root_storage_path / filename
    )
    hdul.writeto(self._internalPaths.root_storage_path / filename, overwrite=True)
    if self.disk_save_level != DISK_SAVE_MODE.EXTREME:
        metrics_path = self._internalPaths.get_path_to("metrics", as_posix=False)
        fig, axis = plt.subplots()
        axis.plot(self.transmittance_wavelengths, self.transmittance_spectra)
        axis.set_xlabel(r"Wavelength [$\AA$]")
        axis.set_ylabel("Transmittance")
        fig.savefig(metrics_path / f"transmittance_{self._associated_subInst}.png")
        plt.close(fig)

load_from_file(root_path: Path, loading_path: str) -> None

TODO: save and load the actual flag to disk!

Parameters:

Name Type Description Default
root_path Path
required
loading_path str
required
Source code in src/ASTRA/template_creation/telluric_templates/Telluric_Template.py 
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
def load_from_file(self, root_path: Path, loading_path: str) -> None:
    """TODO: save and load the actual flag to disk!

    Parameters
    ----------
    root_path
    loading_path

    Returns
    -------

    """
    super().load_from_file(root_path, loading_path)

    with fits.open(loading_path) as hdulist:
        if hdulist[1].header.get("VERSION", "") != __version__:
            logger.warning(
                "Loaded template was not created under the current SBART version. Possible problems may arise",
            )

        self._associated_subInst = hdulist["CONTAM"].header["subInst"]
        self.MAXBERV = hdulist["CONTAM"].header["MAX_BERV"] * kilometer_second
        self.BERVS = [i * kilometer_second for i in hdulist["BERVS"].data.tolist()]

        # By default we load with a window extension mode
        self._extension_mode = TELLURIC_EXTENSION.WINDOW

        try:
            waves = hdulist["Wave"].data
            template = hdulist["Temp"].data
        except:
            waves, template = None, None
        try:
            self.transmittance_wavelengths = hdulist["TRANSMIT_WAVE"].data
            self.transmittance_spectra = hdulist["TRANSMIT_SPECTRA"].data
        except:
            self.transmittance_spectra = None
            self.transmittance_wavelengths = None

        try:
            self.template = template
            self.wavelengths = waves
        except KeyError:
            self.template = None
            self.wavelengths = None

        try:
            self.use_approximated_BERV_correction = hdulist[1].header[
                "HIERARCH APPROX BERV CORRECTION"
            ]
        except KeyError:
            logger.warning("Loading old telluric template with missing keywords")
            self.use_approximated_BERV_correction = False

        self._computed_wave_blocks = False
        self._base_mask = hdulist["CONTAM"].data.tolist()
    self.add_to_status(DISK_LOADED_DATA(f"Loaded data from {loading_path}"))

telluric_from_OHemission

RESOURCES_PATH = Path(__file__).parent.parent.parent / 'resources' module-attribute

OHemissionTelluric

Bases: TelluricTemplate

Telluric mask for OH emission lines.

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_OHemission.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
class OHemissionTelluric(TelluricTemplate):
    """Telluric mask for OH emission lines."""

    _default_params = TelluricTemplate._default_params + DefaultValues(
        SKYcalcPath=UserParam(
            default_value=None,
            constraint=ValueFromDtype((str, Path, type(None))),
            description=("Path to a fits file provided from ESO's skycalc tool for the correct MJD"),
        ),
    )

    method_name = choices.TELLURIC_CREATION_MODE.OHemission.value

    def __init__(
        self,
        subInst: str,
        user_configs: Optional[UI_DICT] = None,
        extension_mode: str = "lines",
        application_mode: str = "removal",
        loaded: bool = False,
    ):
        super().__init__(
            subInst=subInst,
            extension_mode=extension_mode,
            user_configs=user_configs,
            loaded=loaded,
            application_mode=application_mode,
        )

    @custom_exceptions.ensure_invalid_template
    def create_telluric_template(self, dataClass: DataClass, custom_frameID: Optional[int] = None) -> None:
        """Create model for OH emission."""
        try:
            super().create_telluric_template(dataClass, custom_frameID=custom_frameID)
        except custom_exceptions.StopComputationError:
            return

        with fits.open(self._internal_configs["SKYcalcPath"]) as hdu:
            datable = hdu[1].data

        wavelengths, tell_spectra = datable["lam"], datable["flux"]

        template = np.zeros_like(tell_spectra)
        template[np.where(tell_spectra != 0)] = 1
        self.template = template

        # ! no median filtering (might still be needed in the future)
        self._continuum_level = 1.0
        self.wavelengths = wavelengths * 10  # convert to the prevalent wavelength units

        self.transmittance_wavelengths, self.transmittance_spectra = (
            wavelengths,
            tell_spectra,
        )

        self.build_blocks()
        self._compute_wave_blocks()
        self._finish_template_creation()

    def store_metrics(self):
        super().store_metrics()

method_name = choices.TELLURIC_CREATION_MODE.OHemission.value class-attribute instance-attribute

create_telluric_template(dataClass: DataClass, custom_frameID: Optional[int] = None) -> None

Create model for OH emission.

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_OHemission.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
@custom_exceptions.ensure_invalid_template
def create_telluric_template(self, dataClass: DataClass, custom_frameID: Optional[int] = None) -> None:
    """Create model for OH emission."""
    try:
        super().create_telluric_template(dataClass, custom_frameID=custom_frameID)
    except custom_exceptions.StopComputationError:
        return

    with fits.open(self._internal_configs["SKYcalcPath"]) as hdu:
        datable = hdu[1].data

    wavelengths, tell_spectra = datable["lam"], datable["flux"]

    template = np.zeros_like(tell_spectra)
    template[np.where(tell_spectra != 0)] = 1
    self.template = template

    # ! no median filtering (might still be needed in the future)
    self._continuum_level = 1.0
    self.wavelengths = wavelengths * 10  # convert to the prevalent wavelength units

    self.transmittance_wavelengths, self.transmittance_spectra = (
        wavelengths,
        tell_spectra,
    )

    self.build_blocks()
    self._compute_wave_blocks()
    self._finish_template_creation()

store_metrics()

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_OHemission.py 
84
85
def store_metrics(self):
    super().store_metrics()

telluric_from_telfit

RESOURCES_PATH = Path(__file__).parent.parent.parent / 'resources' module-attribute

atmospheric_profiles_coords_dict = {'HARPS': '-70.7-29.3', 'NIRPS': '-70.7-29.3', 'ESPRESSO': '-70.4-24.6', 'POET': '-70.4-24.6', 'HARPSN': '-17.9+28.8', 'CARMENES': '-2.5+37.2', 'CALIB_CARMENES': '-2.5+37.2', 'MAROONX': '-155.5+19.8'} module-attribute

TelfitTelluric

Bases: TelluricTemplate

Create Earth's transmittance spectra from TelFit.

It will be configured with the night with the highest relative humidity.

Atmosphere profiles are downloaded to ensure that we get the best model possible

User parameters: ================================ ================ ================ ================ ================ Parameter name Mandatory Default Value Valid Values Comment ================================ ================ ================ ================ ================ TELFIT_HUMIDITY_THRESHOLD False None Numerical [1] ================================ ================ ================ ================ ================

.. note:: [1] - Value used to enforce a maximum value for the humidity. If this is set any negative value, use the maximum humidity from the loaded observations

.. note:: Also check the User parameters of the parent classes for further customization options of SBART

.. warning:: Telfit is not installed by default with SBART!

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
class TelfitTelluric(TelluricTemplate):
    """Create Earth's transmittance spectra from TelFit.

    It will be configured with the night with the highest   relative humidity.

    Atmosphere profiles are downloaded to ensure that we get the best model possible

    **User parameters:**
    ================================ ================ ================ ================ ================
    Parameter name                      Mandatory      Default Value    Valid Values    Comment
    ================================ ================ ================ ================ ================
    TELFIT_HUMIDITY_THRESHOLD            False          None            Numerical        [1]
    ================================ ================ ================ ================ ================


    .. note::
        [1] - Value used to enforce a maximum value for the humidity. If this is set any negative value, use the
        maximum humidity from the loaded observations


    .. note::
        Also check the **User parameters** of the parent classes for further customization options of SBART

    .. warning::
        Telfit is not installed by default with SBART!

    """

    _default_params = TelluricTemplate._default_params + DefaultValues(
        atmosphere_profile=UserParam(
            "GDAS", constraint=StringValue
        ),  # GDAS / default / path
        FIT_MODEL=UserParam(False, constraint=BooleanValue),
        TELFIT_HUMIDITY_THRESHOLD=UserParam(default_value=-1, constraint=NumericValue),
        FIT_WAVELENGTH_STEP_SIZE=UserParam(0.001, constraint=Positive_Value_Constraint),
        # step size for telluric model wavelengths
        PARAMS_TO_FIT=UserParam(
            ["pressure", "humidity"],
            constraint=ValueFromIterable(
                ["temperature", "pressure", "humidity", "co2", "ch4", "n2o"]
            ),
        ),
        USE_GRID_OF_TRANSMITTANCE=UserParam(
            default_value=False,
            constraint=BooleanValue,
            description=(
                "If True (default False), uses a grid of pre-computed Telfit transmittances to generate"
                "the telluric model. If the grid doesn't exist, it will generate a new one"
            ),
        ),
        GRID_MAIN_PATH=UserParam(
            default_value=None,
            constraint=ValueFromDtype((str, Path, type(None))),
            description=(
                "If not None, it will be used to store the transmittance grid. If None"
                "stores in s_BART 'resources' folder in the location of installation"
            ),
        ),
        IND_WATER_MASK_THRESHOLD=UserParam(  # Ensuring that things don't blow up when storing the fits files (inf will do that)
            default_value=1e8,
            constraint=Positive_Value_Constraint,
            description="Independent masking of water features, using a different threshold than for other molecules",
        ),
    )

    method_name = choices.TELLURIC_CREATION_MODE.telfit.value

    def __init__(
        self,
        subInst: str,
        user_configs: Optional[UI_DICT] = None,
        extension_mode: TELLURIC_EXTENSION = TELLURIC_EXTENSION.LINES,
        application_mode: TELLURIC_APPLICATION_MODE = TELLURIC_APPLICATION_MODE.removal,
        loaded: bool = False,
    ):
        super().__init__(
            subInst=subInst,
            extension_mode=extension_mode,
            user_configs=user_configs,
            loaded=loaded,
            application_mode=application_mode,
        )

        model_components = []

        default_molecule_dict = {
            "co2": 385.34,
            "ch4": 1819.0,
            "n2o": 325.1,
        }

        model_components.append(
            ModelComponent(
                name="temperature",
                initial_guess=283.15,
                bounds=[0, None],
                default_enabled=True,
            ),
        )
        model_components.append(
            ModelComponent(
                name="pressure",
                initial_guess=100,
                bounds=[0, None],
                default_enabled=True,
            ),
        )
        model_components.append(
            ModelComponent(
                name="humidity",
                initial_guess=100,
                bounds=[0, None],
                default_enabled=True,
            ),
        )

        for name, value in default_molecule_dict.items():
            param = ModelComponent(
                name=name,
                initial_guess=value,
                bounds=[0, None],
                default_enabled=True,
                param_type="molecules::extra",
            )
            model_components.append(param)

        for comp in model_components:
            self._fitModel.add_extra_param(comp)

        self.modeler = None

    def _prepare_GDAS_data(self, dataClass, selected_frame):  # pylint: disable=C0103
        logger.info("Preparing GDAS data load!")

        resources_folder = RESOURCES_PATH / "atmosphere_profiles"

        conn = DB_connection()

        if self._internal_configs["atmosphere_profile"] == "GDAS":
            logger.info("Launching GDAS profile downloader")
            instrument, date = selected_frame.inst_name, selected_frame.get_KW_value(
                "ISO-DATE"
            )

            # JD for '2004-12-01T00:00:00', which is the same as the first profile from GDAS
            self._metric_selection_conditions += KEYWORD_condition(
                "BJD", [2453340.5, np.inf]
            )

            logger.warning(
                "Iterating over other possible frames to search for a working reference"
            )
            failed_tests = [self._reference_frameID]
            for kw in ["relative_humidity", "airmass"]:
                metric_to_select, frameIDs = dataClass.collect_KW_observations(
                    KW=kw,
                    subInstruments=[self._associated_subInst],
                    include_invalid=False,
                    conditions=self._metric_selection_conditions,
                    return_frameIDs=True,
                )
                metric_to_select = np.asarray(metric_to_select, dtype=float)
                metric_to_select = list(
                    i if i is not None else np.nan for i in metric_to_select
                )
                if not any(np.isfinite(metric_to_select)):
                    logger.warning(
                        f"Metric {kw} is not finite. Can't use it to select observatioons"
                    )
                    continue
                # Stop if we have finite values!
                break
            else:
                logger.info(f"Using {kw} as the metric to select transmittance profile")

            sort_inds = np.argsort(metric_to_select)
            frames_to_search = np.asarray(frameIDs)[sort_inds][::-1]

            max_search_iterations = len(frameIDs)
            logger.info("Starting loop to retrive GDAS profile")

            found = True
            for attempt_nb in range(max_search_iterations):
                selected_ID = frames_to_search[0]

                date = dataClass.get_KW_from_frameID(frameID=selected_ID, KW="ISO-DATE")

                gdas_filename = construct_gdas_filename(
                    instrument=instrument, datetime=date
                )
                try:
                    data = conn.get_GDAS_profile(gdas_filename=gdas_filename)
                    logger.info("Using cached version of the GDAS profile")
                    break
                except FileNotFoundError:
                    pass

                try:
                    with get_atmospheric_profile(
                        instrument, date, resources_folder
                    ) as gdas:
                        data = np.loadtxt(gdas).copy()
                    self._reference_frameID = selected_ID
                    self._associated_BERV = dataClass.get_KW_from_frameID(
                        KW="BERV", frameID=selected_ID
                    )
                    conn.add_new_profile(
                        gdas_filename=gdas_filename, data=data, instrument=instrument
                    )
                    break  # If we found it, no need to continue
                except KeyError:
                    logger.info(
                        f"{date=} failed to find GDAS profile (Iter: {attempt_nb + 1} / {max_search_iterations})",
                    )
                    failed_tests.append(selected_ID)
                    frames_to_search = frames_to_search[1:]
            else:
                logger.warning(
                    "Couldn't download any of the GDAS profiles. Moving on for the default profile"
                )
                found = False

        if self._internal_configs["atmosphere_profile"] == "default" or not found:
            logger.warning("Using the default atmosphere profile!")
            atmos_profile_file = (
                resources_folder / f"{selected_frame.inst_name}_atmosphere_profile.txt"
            )
            data = np.loadtxt(atmos_profile_file)

        elif os.path.exists(self._internal_configs["atmosphere_profile"]) and not len(
            self._internal_configs["atmosphere_profile"] == 0,
        ):
            # what does this do???
            atmos_profile_file = self._internal_configs["atmosphere_profile"]
            data = np.loadtxt(atmos_profile_file)
        logger.info("Finished setup of GDAS profile")
        return data

    def configure_modeler(self, dataclass) -> None:
        """See https://www.eso.org/sci/software/pipelines/skytools/molecfit#gdas"""
        selected_frame = dataclass.get_frame_by_ID(self._reference_frameID)
        logger.info("Configuring the Telfit modeler for {}", selected_frame)
        data = self._prepare_GDAS_data(
            dataClass=dataclass, selected_frame=selected_frame
        )

        # hPa       m     K       %
        Pres, height, Temp, _ = data.T
        # to ºC
        Temp = Temp - 273.15
        # # calculate dew point
        # !not needed anymore...
        # a, b = 17.62, 243.12
        # α = lambda T, RH: (np.log(RH / 100) + a * T / (b + T))
        # dew = (b * α(Temp, relhum)) / (a - α(Temp, relhum))
        # dew[np.isnan(dew)] = -999

        # """Load the GDAS atmosphere profile following the recipe from
        # https://telfit.readthedocs.io/en/latest/GDAS_atmosphere.html
        # """
        # if self._configs["atmosphere_profile"] == "default":
        #     logger.info("Using default GDAS profile")
        #     atmosphere_filename = os.path.join(
        #         SBART_LOC,
        #         # f"resources/atmosphere_profiles/{self.instrument_info['name']}_atmosphere_profile.txt",
        #         f"resources/atmosphere_profiles/CARMENES_atmosphere_profile.txt",
        #     )
        # else:
        #     atmosphere_filename = self._configs["atmosphere_profile"]
        # if not os.path.exists(atmosphere_filename):
        #     raise FileNotFoundError(
        #         "GDAS profile does not exist in: {}".format(atmosphere_filename)
        #     )

        # logger.info("Loading GDAS atmosphere profile from {}".format(atmosphere_filename))

        # Pres, height, Temp, dew = np.loadtxt(
        #     atmosphere_filename, usecols=(0, 1, 2, 3), skiprows=4, unpack=True
        # )

        # Sort the arrays by height.
        sorter = np.argsort(height)
        height = height[sorter]
        Pres = Pres[sorter]
        Temp = Temp[sorter]
        # dew = dew[sorter]

        # Convert dew point temperature to ppmv
        Pw = 6.116441 * 10 ** (7.591386 * Temp / (Temp + 240.7263))
        h2o = Pw / (Pres - Pw) * 1e6

        # Unit conversion
        height = height / 1000.0
        Temp = Temp + 273.15
        self.modeler.EditProfile("Temperature", height, Temp)
        self.modeler.EditProfile("Pressure", height, Pres)
        self.modeler.EditProfile("H2O", height, h2o)

        logger.info("Finished configurating the modeler")

    def _generate_telluric_model(
        self,
        dataClass,
        model_parameters,
        OBS_properties,
        instrument: str,
        fixed_params=None,
        grid=None,
    ) -> tuple[np.ndarray, np.ndarray]:
        """COmpute the transmittance spectra using tellfit.

        Parameters
        ----------
        wavelengths : np.ndarray
            wavelengths for the reference observation

        Returns
        -------
        Tuple[np.ndarray, np.ndarray]
            Wavelengths and transmittance arrays

        """
        control_dict = dict(zip(self._fitModel.get_enabled_params(), model_parameters))
        if fixed_params is not None:
            control_dict = {**control_dict, **fixed_params}

        if grid is None:
            # ask for a model slightly larger than the wavelength coverage of the instrument
            lowfreq = 1e7 / (OBS_properties["wavelength_coverage"][1] + 50)
            highfreq = 1e7 / (OBS_properties["wavelength_coverage"][0] - 50)
        else:
            lowfreq = 1e7 / grid.max()
            highfreq = 1e7 / grid.min()

        generate_new_grid = False
        if self._internal_configs["USE_GRID_OF_TRANSMITTANCE"]:
            if self._internal_configs["GRID_MAIN_PATH"] is not None:
                loc = Path(self._internal_configs["GRID_MAIN_PATH"])
            else:
                loc = RESOURCES_PATH / "atmosphere_grid" / instrument
                if not loc.exists() or len(list(loc.iterdir())) == 0:
                    logger.warning(
                        "Folder with grids is not available. Generating it now"
                    )
                    loc.mkdir(exist_ok=True, parents=True)
                    # Not the greatest code, but this ensures that the modeler
                    # is ready to go when we launch the telfit grid computation

                    try:
                        import telfit

                    except ModuleNotFoundError as e:
                        msg = "Telfit is not istalled"
                        raise custom_exceptions.InternalError(msg) from e

                    if self.modeler is None:
                        self.modeler = telfit.Modeler(print_lblrtm_output=False)
                        self.configure_modeler(dataClass)

                    self._construct_grid(
                        loc,
                        OBS_properties=OBS_properties,
                        control_dict=control_dict,
                        freqs=[lowfreq, highfreq],
                        grid=None,  # TODO: see if we can reduce the number of points in here
                    )
            logger.debug(f"Using grid of trasmittance values, searching in {loc}")
            # This needs to be updated if the grid changes
            closest_temp = int(np.ceil(control_dict["temperature"] / 2) * 2)
            closest_rhum = int(np.ceil(control_dict["humidity"] / 5) * 5)
            closest_airmass = np.round(
                np.ceil(OBS_properties["airmass"] / 0.1) * 0.1, 2
            )

            desired_file = (
                loc
                / f"transmittance_{closest_airmass}_{closest_rhum}_{closest_temp}.txt"
            )
            if not desired_file.exists():
                params = (closest_temp, closest_rhum, closest_airmass)

                logger.critical(
                    f"Using a grid of transmittance and missing the relevant profiles ({params}). Moving onwards to nominal run"
                )
            else:
                logger.debug("Found pre-computed transmittance spectra, serving it")
                data = np.loadtxt(desired_file)
                return data[:, 0], data[:, 1]

        try:
            import telfit

        except ModuleNotFoundError as e:
            msg = "Telfit is not istalled"
            raise custom_exceptions.InternalError(msg) from e

        if self.modeler is None:
            self.modeler = telfit.Modeler(print_lblrtm_output=False)
            self.configure_modeler(dataClass)

        observatory = OBS_properties["EarthLocation"]

        m = self.modeler.MakeModel(
            lowfreq=lowfreq,
            highfreq=highfreq,
            vac2air=False,
            angle=np.rad2deg(np.arccos(1 / OBS_properties["airmass"])),
            lat=observatory.lat.value,
            alt=observatory.height.value / 1e3,
            resolution=OBS_properties["resolution"],
            wavegrid=grid,
            **control_dict,
        )
        return m.x, m.y

    def _construct_grid(self, loc, OBS_properties, control_dict, freqs, grid):
        # If we update the grid steps, we will need to also change the search for the closest file
        airmass_values = np.arange(start=1, stop=2.1, step=0.1)
        humidity_values = np.arange(start=0, stop=35, step=5, dtype=int)
        temperature_values = np.arange(
            278, 294.1, step=2, dtype=int
        )  # from 5 to 20 Celsius

        combinations = product(airmass_values, humidity_values, temperature_values)

        target = partial(
            fun,
            modeler=self.modeler,
            loc=loc,
            OBS_properties=OBS_properties,
            freqs=freqs,
            control_dict=control_dict,
            grid=grid,
        )

        with Pool(processes=4) as p:
            p.map(target, combinations)

    @custom_exceptions.ensure_invalid_template
    def create_telluric_template(
        self, dataClass: DataClass, custom_frameID: Optional[int] = None
    ) -> None:
        """Create a telluric template from a TelFit transmission spectra [1].

        The model is created for the date in which the reference observation was made.

        It estimates the continuum level and classifies each point that shows a
        decrease of 1% as a telluric line. Furthermore, it creates a window of
        6 points to each side of a detection, to attempt to pick up the wings of
        the telluric line.

        Parameters
        ----------
        dataClass:
            DataClass object
        custom_frameID :
            If Not None, does not search for the "optimal" frameID to use as a basis

        Returns
        -------
        numpy.ndarray
            Telluric (binary) spectrum, for the wavelengths present in the input
            array

        Notes
        -----
        [1] https://github.com/kgullikson88/Telluric-Fitter

        """
        try:
            super().create_telluric_template(dataClass, custom_frameID=custom_frameID)
        except custom_exceptions.StopComputationError:
            return

        OBS_properties = {
            "airmass": dataClass.get_KW_from_frameID(
                "airmass", self._reference_frameID
            ),
            **dataClass.get_instrument_information(),
        }

        # Note: by default, when removing features, the "fit results" are the initial guesses
        #    Later on, if needed, this can be easily changed to fit before removing!

        parameter_values = self._fitModel.get_fit_results_from_frameID(
            frameID=self._reference_frameID,
            allow_disabled=True,
        )
        instrument = dataClass.get_frame_by_ID(self._reference_frameID).inst_name

        names = self._fitModel.get_component_names(include_disabled=True)
        logger.debug(f"Parameters in use: {names}")
        logger.info(f"Using params: {parameter_values}")

        # This is the model with every molecule in the dataset. If there is no individual masking for water
        # This will work with the same threshold for every single element
        threshold = self._internal_configs["continuum_percentage_drop"]

        USE_INDEPENDENT_WATER = self._internal_configs["IND_WATER_MASK_THRESHOLD"] < 1e6
        if USE_INDEPENDENT_WATER:
            logger.debug("There is water")

            threshold = self._internal_configs["IND_WATER_MASK_THRESHOLD"]

        wavelengths, tell_spectra = self._generate_telluric_model(
            dataClass=dataClass,
            model_parameters={},
            OBS_properties=OBS_properties,
            fixed_params=dict(zip(names, parameter_values)),
            instrument=instrument,
        )

        self.template = create_binary_template(
            transmittance=tell_spectra,
            continuum_level=1.0,
            percentage_drop=threshold,
        )

        # ! no median filtering (might still be needed in the future)
        self._continuum_level = 1.0
        self.wavelengths = wavelengths * 10  # convert to the prevalent wavelength units

        self.transmittance_wavelengths, self.transmittance_spectra = (
            wavelengths,
            tell_spectra,
        )

        if USE_INDEPENDENT_WATER:
            # We wanted an individual mask for water, which means that now we have to
            # create the mask with the correct transmittance for the other elements
            logger.debug("Starting full spectra model")
            parameters = {}
            for key, value in zip(names, parameter_values):
                parameters[key] = value
                if key == "humidity":
                    parameters[key] = 0

            logger.warning(f"Original value: {tell_spectra.shape}")

            waves, tell_spectra = self._generate_telluric_model(
                dataClass=dataClass,
                model_parameters={},
                OBS_properties=OBS_properties,
                fixed_params=parameters,
                grid=wavelengths,
                instrument=instrument,
            )

            logger.warning(f"New value: {tell_spectra.shape}")
            self.template += create_binary_template(
                transmittance=tell_spectra,
                continuum_level=1.0,
                percentage_drop=self._internal_configs["continuum_percentage_drop"],
            )

            # Ensure that we don't have values grater than 1, to keep consistency
            self.template[np.where(self.template > 1)] = 1

        self.build_blocks()
        self._compute_wave_blocks()
        self._finish_template_creation()

    def _generate_model_parameters(self, dataClass):
        """Custom generation of the priors for the humidity, temperature and pressure

        Parameters
        ----------
        dataClass

        Returns
        -------

        """
        super()._generate_model_parameters(dataClass)
        for frameID in dataClass.get_frameIDs_from_subInst(
            self._associated_subInst, include_invalid=False
        ):
            frame = dataClass.get_frame_by_ID(frameID)
            initial_guess = {
                "humidity": frame.get_KW_value("relative_humidity"),
                "temperature": frame.get_KW_value("ambient_temperature"),
                "pressure": dataClass.get_instrument_information()["site_pressure"],
            }

            if not np.isfinite(initial_guess["humidity"]):
                initial_guess["humidity"] = 25
                logger.warning(
                    "Relative humidity is not finite. Using default value of {}%",
                    initial_guess["humidity"],
                )

            if (
                0
                <= self._internal_configs["TELFIT_HUMIDITY_THRESHOLD"]
                < initial_guess["humidity"]
            ):
                initial_guess["humidity"] = self._internal_configs[
                    "TELFIT_HUMIDITY_THRESHOLD"
                ]
                user_cap = self._internal_configs["TELFIT_HUMIDITY_THRESHOLD"]
                curr_val = initial_guess["humidity"]
                logger.warning(
                    f"Relative humidity ({curr_val}) is above the user-provided threshold ({user_cap}). Falling back to it",
                )

            if not np.isfinite(initial_guess["temperature"]):
                initial_guess["temperature"] = 290.5
                logger.warning(
                    "Ambient temperature is not finite. Using default value of {}K",
                    initial_guess["temperature"],
                )

            self._fitModel.update_params_initial_guesses(
                frameID=frame.frameID, guesses=initial_guess
            )

            if self.for_feature_removal:
                finals = [
                    comps.get_initial_guess(frameID, True)
                    for comps in self._fitModel.get_enabled_components()
                ]
                self._fitModel.store_frameID_results(
                    frameID, finals, result_flag=SUCCESS
                )

        for param_name in self._fitModel.get_enabled_params():
            if param_name not in self._internal_configs["PARAMS_TO_FIT"]:
                logger.info(
                    "{} not fitting {}. Fixing it to initial guess",
                    self.name,
                    param_name,
                )
                self._fitModel.disable_param(param_name)

    def store_metrics(self):
        super().store_metrics()

        if (
            self.for_feature_correction
            or self.disk_save_level == DISK_SAVE_MODE.DISABLED
        ):
            metrics_path = self._internalPaths.get_path_to("metrics", as_posix=False)
            parameter_values = self._fitModel.get_fit_results_from_frameID(
                frameID=self._reference_frameID,
                allow_disabled=True,
            )
            names = self._fitModel.get_component_names(include_disabled=True)

            with open(
                metrics_path / f"telfit_info_{self._associated_subInst}.txt", mode="w"
            ) as to_write:
                for nam, val in zip(names, parameter_values):
                    to_write.write(f"{nam}:  {val}\n")

method_name = choices.TELLURIC_CREATION_MODE.telfit.value class-attribute instance-attribute

modeler = None instance-attribute

configure_modeler(dataclass) -> None

See https://www.eso.org/sci/software/pipelines/skytools/molecfit#gdas

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
def configure_modeler(self, dataclass) -> None:
    """See https://www.eso.org/sci/software/pipelines/skytools/molecfit#gdas"""
    selected_frame = dataclass.get_frame_by_ID(self._reference_frameID)
    logger.info("Configuring the Telfit modeler for {}", selected_frame)
    data = self._prepare_GDAS_data(
        dataClass=dataclass, selected_frame=selected_frame
    )

    # hPa       m     K       %
    Pres, height, Temp, _ = data.T
    # to ºC
    Temp = Temp - 273.15
    # # calculate dew point
    # !not needed anymore...
    # a, b = 17.62, 243.12
    # α = lambda T, RH: (np.log(RH / 100) + a * T / (b + T))
    # dew = (b * α(Temp, relhum)) / (a - α(Temp, relhum))
    # dew[np.isnan(dew)] = -999

    # """Load the GDAS atmosphere profile following the recipe from
    # https://telfit.readthedocs.io/en/latest/GDAS_atmosphere.html
    # """
    # if self._configs["atmosphere_profile"] == "default":
    #     logger.info("Using default GDAS profile")
    #     atmosphere_filename = os.path.join(
    #         SBART_LOC,
    #         # f"resources/atmosphere_profiles/{self.instrument_info['name']}_atmosphere_profile.txt",
    #         f"resources/atmosphere_profiles/CARMENES_atmosphere_profile.txt",
    #     )
    # else:
    #     atmosphere_filename = self._configs["atmosphere_profile"]
    # if not os.path.exists(atmosphere_filename):
    #     raise FileNotFoundError(
    #         "GDAS profile does not exist in: {}".format(atmosphere_filename)
    #     )

    # logger.info("Loading GDAS atmosphere profile from {}".format(atmosphere_filename))

    # Pres, height, Temp, dew = np.loadtxt(
    #     atmosphere_filename, usecols=(0, 1, 2, 3), skiprows=4, unpack=True
    # )

    # Sort the arrays by height.
    sorter = np.argsort(height)
    height = height[sorter]
    Pres = Pres[sorter]
    Temp = Temp[sorter]
    # dew = dew[sorter]

    # Convert dew point temperature to ppmv
    Pw = 6.116441 * 10 ** (7.591386 * Temp / (Temp + 240.7263))
    h2o = Pw / (Pres - Pw) * 1e6

    # Unit conversion
    height = height / 1000.0
    Temp = Temp + 273.15
    self.modeler.EditProfile("Temperature", height, Temp)
    self.modeler.EditProfile("Pressure", height, Pres)
    self.modeler.EditProfile("H2O", height, h2o)

    logger.info("Finished configurating the modeler")

create_telluric_template(dataClass: DataClass, custom_frameID: Optional[int] = None) -> None

Create a telluric template from a TelFit transmission spectra [1].

The model is created for the date in which the reference observation was made.

It estimates the continuum level and classifies each point that shows a decrease of 1% as a telluric line. Furthermore, it creates a window of 6 points to each side of a detection, to attempt to pick up the wings of the telluric line.

Parameters:

Name Type Description Default
dataClass DataClass

DataClass object

 required
custom_frameID Optional[int]

If Not None, does not search for the "optimal" frameID to use as a basis

 None

Returns:

Type Description
ndarray

Telluric (binary) spectrum, for the wavelengths present in the input array
Notes

[1] https://github.com/kgullikson88/Telluric-Fitter

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
@custom_exceptions.ensure_invalid_template
def create_telluric_template(
    self, dataClass: DataClass, custom_frameID: Optional[int] = None
) -> None:
    """Create a telluric template from a TelFit transmission spectra [1].

    The model is created for the date in which the reference observation was made.

    It estimates the continuum level and classifies each point that shows a
    decrease of 1% as a telluric line. Furthermore, it creates a window of
    6 points to each side of a detection, to attempt to pick up the wings of
    the telluric line.

    Parameters
    ----------
    dataClass:
        DataClass object
    custom_frameID :
        If Not None, does not search for the "optimal" frameID to use as a basis

    Returns
    -------
    numpy.ndarray
        Telluric (binary) spectrum, for the wavelengths present in the input
        array

    Notes
    -----
    [1] https://github.com/kgullikson88/Telluric-Fitter

    """
    try:
        super().create_telluric_template(dataClass, custom_frameID=custom_frameID)
    except custom_exceptions.StopComputationError:
        return

    OBS_properties = {
        "airmass": dataClass.get_KW_from_frameID(
            "airmass", self._reference_frameID
        ),
        **dataClass.get_instrument_information(),
    }

    # Note: by default, when removing features, the "fit results" are the initial guesses
    #    Later on, if needed, this can be easily changed to fit before removing!

    parameter_values = self._fitModel.get_fit_results_from_frameID(
        frameID=self._reference_frameID,
        allow_disabled=True,
    )
    instrument = dataClass.get_frame_by_ID(self._reference_frameID).inst_name

    names = self._fitModel.get_component_names(include_disabled=True)
    logger.debug(f"Parameters in use: {names}")
    logger.info(f"Using params: {parameter_values}")

    # This is the model with every molecule in the dataset. If there is no individual masking for water
    # This will work with the same threshold for every single element
    threshold = self._internal_configs["continuum_percentage_drop"]

    USE_INDEPENDENT_WATER = self._internal_configs["IND_WATER_MASK_THRESHOLD"] < 1e6
    if USE_INDEPENDENT_WATER:
        logger.debug("There is water")

        threshold = self._internal_configs["IND_WATER_MASK_THRESHOLD"]

    wavelengths, tell_spectra = self._generate_telluric_model(
        dataClass=dataClass,
        model_parameters={},
        OBS_properties=OBS_properties,
        fixed_params=dict(zip(names, parameter_values)),
        instrument=instrument,
    )

    self.template = create_binary_template(
        transmittance=tell_spectra,
        continuum_level=1.0,
        percentage_drop=threshold,
    )

    # ! no median filtering (might still be needed in the future)
    self._continuum_level = 1.0
    self.wavelengths = wavelengths * 10  # convert to the prevalent wavelength units

    self.transmittance_wavelengths, self.transmittance_spectra = (
        wavelengths,
        tell_spectra,
    )

    if USE_INDEPENDENT_WATER:
        # We wanted an individual mask for water, which means that now we have to
        # create the mask with the correct transmittance for the other elements
        logger.debug("Starting full spectra model")
        parameters = {}
        for key, value in zip(names, parameter_values):
            parameters[key] = value
            if key == "humidity":
                parameters[key] = 0

        logger.warning(f"Original value: {tell_spectra.shape}")

        waves, tell_spectra = self._generate_telluric_model(
            dataClass=dataClass,
            model_parameters={},
            OBS_properties=OBS_properties,
            fixed_params=parameters,
            grid=wavelengths,
            instrument=instrument,
        )

        logger.warning(f"New value: {tell_spectra.shape}")
        self.template += create_binary_template(
            transmittance=tell_spectra,
            continuum_level=1.0,
            percentage_drop=self._internal_configs["continuum_percentage_drop"],
        )

        # Ensure that we don't have values grater than 1, to keep consistency
        self.template[np.where(self.template > 1)] = 1

    self.build_blocks()
    self._compute_wave_blocks()
    self._finish_template_creation()

store_metrics()

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
def store_metrics(self):
    super().store_metrics()

    if (
        self.for_feature_correction
        or self.disk_save_level == DISK_SAVE_MODE.DISABLED
    ):
        metrics_path = self._internalPaths.get_path_to("metrics", as_posix=False)
        parameter_values = self._fitModel.get_fit_results_from_frameID(
            frameID=self._reference_frameID,
            allow_disabled=True,
        )
        names = self._fitModel.get_component_names(include_disabled=True)

        with open(
            metrics_path / f"telfit_info_{self._associated_subInst}.txt", mode="w"
        ) as to_write:
            for nam, val in zip(names, parameter_values):
                to_write.write(f"{nam}:  {val}\n")

fun(params, modeler, loc, OBS_properties, freqs, control_dict, grid) -> None

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def fun(params, modeler, loc, OBS_properties, freqs, control_dict, grid) -> None:
    airmass, humidity, temperature = params

    fname = f"transmittance_{airmass:.1f}_{humidity:.0f}_{temperature:.0f}.txt"

    lowfreq, highfreq = freqs
    observatory = OBS_properties["EarthLocation"]

    control_dict["temperature"] = temperature
    control_dict["humidity"] = humidity

    m = modeler.MakeModel(
        lowfreq=lowfreq,
        highfreq=highfreq,
        vac2air=False,
        angle=np.rad2deg(np.arccos(1 / airmass)),
        lat=observatory.lat.value,
        alt=observatory.height.value / 1e3,
        resolution=OBS_properties["resolution"],
        wavegrid=grid,
        **control_dict,
    )
    wave, transmit = m.x, m.y
    np.savetxt(fname=loc / fname, X=np.c_[wave, transmit])

construct_gdas_filename(instrument, datetime)

Construct the filename as in GDAS archive

Args: instrument (str): Instrument name datetime (str): Date and time, as in the header of fits files

Raises: custom_exceptions.InvalidConfiguration: description

Returns: type: description

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def construct_gdas_filename(instrument, datetime):
    """Construct the filename as in GDAS archive

    Args:
        instrument (str): Instrument name
        datetime (str): Date and time, as in the header of fits files

    Raises:
        custom_exceptions.InvalidConfiguration: _description_

    Returns:
        _type_: _description_

    """
    try:
        coords = atmospheric_profiles_coords_dict[instrument]
    except KeyError as exc:
        raise custom_exceptions.InvalidConfiguration(
            f"Telfit template does not support {instrument}. Available instruments: {list(atmospheric_profiles_coords_dict.keys())}",
        ) from exc

    date, time = datetime.split("T")
    hour = int(time.split(":")[0])
    # GDAS profiles are available every 3 hours, find the closest one
    available = [0, 3, 6, 9, 12, 15, 18, 21]
    closest_hour = min(available, key=lambda list_value: abs(list_value - hour))

    return f"C{coords}D{date}T{closest_hour:02d}.gdas"

download_gdas_archive(archive_name, storage_path)

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
114
115
116
117
def download_gdas_archive(archive_name, storage_path):
    url = "https://ftp.eso.org/pub/dfs/pipelines/skytools/molecfit/gdas/"
    url += archive_name
    urllib.request.urlretrieve(url, storage_path)

get_atmospheric_profile(instrument, datetime, storage_folder)

Source code in src/ASTRA/template_creation/telluric_templates/telluric_from_telfit.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
@contextlib.contextmanager
def get_atmospheric_profile(instrument, datetime, storage_folder):
    try:
        coords = atmospheric_profiles_coords_dict[instrument]
    except KeyError as exc:
        raise custom_exceptions.InvalidConfiguration(
            f"Telfit template does not support {instrument}. Available instruments: {list(atmospheric_profiles_coords_dict.keys())}",
        ) from exc

    archive_name = f"gdas_profiles_C{coords}.tar.gz"
    profile_archive = os.path.join(storage_folder, f"gdas_profiles_C{coords}.tar.gz")
    # NOTE: these archives have thousands of files, maybe don't extract them
    exists = os.path.exists(profile_archive)

    if not exists or file_older_than(profile_archive, 7):  # !check if 7 is ok
        # download it from ESO
        logger.info("Downloading GDAS tar file from archive")
        download_gdas_archive(archive_name, profile_archive)

    msg = f"Looking for GDAS atmosphere profile in {profile_archive}"
    logger.info(msg)

    gdas_file = construct_gdas_filename(instrument=instrument, datetime=datetime)

    logger.info(f"Closest GDAS atmospheric profile is {gdas_file}")

    logger.info("Extracting file from the archive")
    tf = tarfile.open(profile_archive, "r:gz")

    try:
        yield tf.extractfile(gdas_file)
    finally:
        tf.close()