Dl3 reader for LST data #7

mstrzys · 2023-06-25T06:22:10Z

DL3 reader; first for LST data only; to add a DL3 reader also for MAGIC data the filetype check needs a rework, which will go into a separate PR.

- added an event class for LST dl3 files - tested that class can read LST dl3 files - entire chain not tested - DL3 format for other instruments (e.g. MAGIC) can be slightly different and mya require a different reader ToDo: - MAGIC DL3 reader - test entire chain with LST DL3 files

Found and corrected the following issues, some concern also previous commits to the main branch: - issue of the unti of the event time variable - wrong unit of the ra event variable - dl2 file reader had an error if the ra variable of the pointing position was not precomputed - wrong overwriting of the setter function for the bkgmap maker in the bkg classes

to set attribute bkg_maps instead of the correct _bkg_maps

src/pybkgmodel/camera.py

src/pybkgmodel/processing.py

IevgenVovk · 2023-07-14T09:36:15Z

In general PR looks okey save the inserted docstring placeholders (useless unless replaced with actual docs). I would encourage in the future not to mix topical PRs like this with bug and style fixes in the unrelated classes and modules.

src/pybkgmodel/data.py

- Modified the reader so it can read DL3 file regardless the instrument. - Data units are read from hdu header as long as provided. - Tested with LST and MAGIC DL3 files. TODO: - Observatory location not part of the mandatory header keywords. May need to be provided in input file by the user. Will be part of future pull request - Drift method implemented but not tested due to lack of corresponding data.

performance - added some functionalities suggested by Lukas Nickel, but adapted them to the current version - changed the function for retrieving the Obs ID from filename to reading it from the hdu header - changed the is_compatible, so it also identifies compressed files - added the erfaAstromInterpolator for the coordinate transform for speed up - added docstrings to all functions related to the DL3 reader

modes are not to well defined e.g. the H.E.S.S. public data use WOBBLE as the obs mode. Hence, added it to support such files until the allowed modes are better defined.

mstrzys · 2024-03-21T09:55:30Z

The latest changes allow the DL3 reader to not only read the LST files, but DL3 files in general. I tested it with LST, H.E.S.S. (public data release), and MAGIC DL3 files.

Furthermore added the suggestions from Lukas regarding the interpolator for the coordinate transform and the changes regarding the Obs ID retrieval; changed the compatibility check to not reject compressed files as well.

Current version is ready for review.

src/pybkgmodel/data.py

from the corresponding PR (#7). - change name from Dl2 to DL2 to agree with the DL3 naming - changed to way the event times are read from DL3 files to be consistent with the 'mjd' naming of the dictionary - the functionality of the code was not changed

src/pybkgmodel/data.py

maxnoe · 2024-07-26T09:25:48Z

Hi,

I got a notification here and had a cursory glance, but what I saw made me very nervous about timing information.

You add and subtract raw mjd and unix values at different points in the code here instead of using astropy Time and TimeDelta objects.

This will result in timestamps that are off by the number of leap seconds between those two times. Never do that.

You are also completely neglecting the time scale.

maxnoe · 2024-07-26T09:26:34Z

src/pybkgmodel/data.py

+                            event_data[name] *= u.one
+
+                # Event times need to be converted from Instrument reference epoch
+                ref_epoch = astropy.time.Time(evt_head['MJDREFI']+evt_head['MJDREFF'], format='mjd')


Use the two argument form of astropy.time.Time, otherwise you will have precision issues:

ref_epoch = astropy.time.Time(evt_head['MJDREFI'], evt_head['MJDREFF'], format='mjd')

also, you need to read the time scale

added, though I found that in the header of the DL3 files the TIMESYS is given in all-caps whereas astropy prefers lower cases, is there a reason for why it is all-caps in the LST files?

Because that's what the FITS standard is prescribing

maxnoe · 2024-07-26T09:27:43Z

src/pybkgmodel/data.py

+                # Event times need to be converted from Instrument reference epoch
+                ref_epoch = astropy.time.Time(evt_head['MJDREFI']+evt_head['MJDREFF'], format='mjd')
+
+                event_data['mjd'] = astropy.time.Time((evt_data['TIME'].to_numpy() 


Do not do this! this results in times that are off by the number of leap seconds between the reference epoch and the event time.

The correct solution is (ref_epoch + evt_data["TIME"].quantity).mjd

maxnoe · 2024-07-26T09:28:53Z

src/pybkgmodel/data.py

+                                                      ).mjd * u.d
+
+                # Adapt
+                evt_time = astropy.time.Time(event_data['mjd'], format='mjd')


now you create a time object a second time.

Just to evt_time = ref_epoch + evt_data["TIME"].quantity and then evt_data["mjd"] = evt_time.mjd

All of this is also already implemented in Gammapy, e.g. here: https://github.com/gammapy/gammapy/blob/main/gammapy/utils/time.py

Simplified according to your suggestion.

maxnoe · 2024-07-26T09:35:11Z

src/pybkgmodel/data.py

+            try:
+                evt_hdu = input_file["EVENTS"]
+                evt_head = evt_hdu.header
+                evt_data = pandas.DataFrame(evt_hdu.data)


use an astropy table to have support for the quantities and metadata stored in the hdu

Thanks for that hint, simplified it a lot. Adapted.

- changed the time handling following Max' recommendations in the discussion of the PR - changed the way the fits file compatibility is checked following Michele's recommendations in the discussion of the PR - changed the setter for the bkg_map_makers following Ievgen's recommendation TODO: - Test with data

maxnoe · 2024-07-31T09:06:28Z

src/pybkgmodel/data.py

-                                                       + ref_epoch.unix),
-                                                      scale='utc',
-                                                      format='unix'
+                event_data['mjd'] = astropy.time.Time((evt_data['TIME'].quantity 


The outer time here is unnecessary.

The result of ref_epoch + evt_data["TIME"].quantity is a time, you are just making a copy again.

maxnoe · 2024-07-31T09:07:11Z

src/pybkgmodel/data.py

                # TODO: current observatory location only La Palma, no mandatory header keyword
                obs_loc  = EarthLocation(lat=28.761758*u.deg,
                                        lon=-17.890659*u.deg,
                                        height=2200*u.m)

                if evt_head['OBS_MODE'] in ('POINTING', 'WOBBLE'):

-                    alt_az_frame = AltAz(obstime=evt_time,
+                    alt_az_frame = AltAz(obstime=event_data['mjd'],


I don't think you can pass a raw mjd value to AltAz(obstime=), this needs to be a Time instance.

Indeed, corrected.

- simplified the reading of the data from the fits files with astropy.Tables follwing Max' suggestion - fixed the obstime object in the AltAz frame definition following Max' comment in the discussion of the PR

- fixed the setter function of the bkg_map_maker

- debugging the setter for the bkg_map_maker - current version not running

- intermediate commit for checking the updated setter

- ready for further review - setter for bkg_map_maker fixed - time scale added to ref_epoch - checked that the code runs with with LST-1 data

maxnoe · 2024-08-01T09:10:44Z

src/pybkgmodel/camera.py

@@ -375,7 +375,7 @@ def get_pixel_coords(self):

        x = (self.xedges[1:] + self.xedges[:-1]) / 2
        y = (self.yedges[1:] + self.yedges[:-1]) / 2
-        xx, yy = np.meshgrid(x, y, indexing='ij')
+        xx, yy = numpy.meshgrid(x, y, indexing='ij')


I think this is the wrong resolution of a merge conflict?

The import at the top is import numpy as np

As entire merge went wrong. I reset the merge. New version should not have this mistake.

maxnoe · 2024-08-01T09:10:51Z

src/pybkgmodel/data.py

@@ -312,11 +314,11 @@ def load_events(cls, file_name, cuts):
                        event_data[name] = data[key]

                    # Post processing
-                    event_data['true_zd'] = np.degrees(event_data['true_zd'])
-                    event_data['true_az'] = np.degrees(event_data['true_az'])
+                    event_data['true_zd'] = numpy.degrees(event_data['true_zd'])


same as above

As entire merge went wrong. I reset the merge. New version should not have this mistake.

maxnoe · 2024-08-01T09:13:03Z

@mstrzys I think something went seriously wrong in your last merge commit, all the changes related to reading LST DL3 data were removed

- changed the remaining numpy's to np, not covered by the merger

- removed the np.bool for bool inherited from the main branch - tested with LST1 data and ensure that it runs properly

maxnoe · 2024-09-03T10:20:30Z

src/pybkgmodel/data.py

+                                              )
+
+                evt_time = evt_data['TIME'].quantity + ref_epoch
+                event_data['mjd'] = evt_time.mjd * u.d


I am not sure, but I think this mjd will now depend on the timescale of the input file, and I guess you use this value internally and assume it's in UTC or some other fixed timescale?

Why are you storing the mjd in any case and not the full time object?

maxnoe · 2024-09-03T10:21:40Z

src/pybkgmodel/data.py

+                                                            * len(event_data['pointing_zd'])
+                                                            ) * u.deg
+
+                elif evt_head['OBS_MODE'] == 'DRIFT':


Just raise an error to let users know this is not supported instead of providing preliminary or invalid data (same for the next else block, this should be an actual error, not a print statement)

Marcel Strzys added 5 commits June 13, 2023 19:58

- just made the LstDl3event class usable by the higher classes

475cdba

bug fix: wrong attribute in the stacked maps class: get_maps() attempted

6a91533

to set attribute bkg_maps instead of the correct _bkg_maps

- just remove unnecessary comma

d0fa65a

mstrzys requested a review from IevgenVovk June 25, 2023 06:22

mstrzys self-assigned this Jun 25, 2023

IevgenVovk reviewed Jul 14, 2023

View reviewed changes

src/pybkgmodel/camera.py Show resolved Hide resolved

IevgenVovk reviewed Jul 14, 2023

View reviewed changes

src/pybkgmodel/processing.py Outdated Show resolved Hide resolved

maxnoe reviewed Jul 18, 2023

View reviewed changes

src/pybkgmodel/data.py Outdated Show resolved Hide resolved

LukasNickel force-pushed the dl3_reader branch from cc52930 to d0fa65a Compare July 18, 2023 15:11

Marcel Strzys added 4 commits March 19, 2024 18:42

<feat> adapted the model.py file to use the new DL3 reader.

0358b26

<feat> Added the Obs mode WOBBLE in addition to POINTING. Seems the Obs

61e2e9a

modes are not to well defined e.g. the H.E.S.S. public data use WOBBLE as the obs mode. Hence, added it to support such files until the allowed modes are better defined.

IevgenVovk reviewed Apr 19, 2024

View reviewed changes

HealthyPear reviewed Jun 12, 2024

View reviewed changes

src/pybkgmodel/data.py Outdated Show resolved Hide resolved

maxnoe reviewed Jul 26, 2024

View reviewed changes

maxnoe reviewed Jul 31, 2024

View reviewed changes

mstrzys added 2 commits July 31, 2024 18:49

FIX:

eab77dd

- simplified the reading of the data from the fits files with astropy.Tables follwing Max' suggestion - fixed the obstime object in the AltAz frame definition following Max' comment in the discussion of the PR

FIX:

98e48fc

- fixed the setter function of the bkg_map_maker

mstrzys added 3 commits July 31, 2024 19:25

FIX:

8966b98

- debugging the setter for the bkg_map_maker - current version not running

FIX:

84aeceb

- intermediate commit for checking the updated setter

FIX:

a5f98dd

- ready for further review - setter for bkg_map_maker fixed - time scale added to ref_epoch - checked that the code runs with with LST-1 data

IevgenVovk approved these changes Aug 1, 2024

View reviewed changes

mstrzys requested a review from maxnoe August 1, 2024 08:31

maxnoe reviewed Aug 1, 2024

View reviewed changes

mstrzys force-pushed the dl3_reader branch from ba680ce to a5f98dd Compare August 5, 2024 03:36

Marcel Strzys added 3 commits August 5, 2024 13:06

Merge branch 'main' into dl3_reader

64a3e9c

FIX:

af0e3a9

- changed the remaining numpy's to np, not covered by the merger

FIX:

33ea973

- removed the np.bool for bool inherited from the main branch - tested with LST1 data and ensure that it runs properly

mstrzys requested a review from maxnoe August 5, 2024 06:01

maxnoe reviewed Sep 3, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Dl3 reader for LST data #7

Dl3 reader for LST data #7

mstrzys commented Jun 25, 2023

IevgenVovk commented Jul 14, 2023

mstrzys commented Mar 21, 2024 •

edited

Loading

maxnoe commented Jul 26, 2024 •

edited

Loading

maxnoe Jul 26, 2024

maxnoe Jul 26, 2024

mstrzys Jul 31, 2024

maxnoe Aug 1, 2024

maxnoe Jul 26, 2024

mstrzys Jul 31, 2024

maxnoe Jul 26, 2024

maxnoe Jul 26, 2024

mstrzys Jul 31, 2024

maxnoe Jul 26, 2024

mstrzys Jul 31, 2024

maxnoe Jul 31, 2024

mstrzys Jul 31, 2024

maxnoe Jul 31, 2024

mstrzys Jul 31, 2024

maxnoe Aug 1, 2024

mstrzys Aug 5, 2024

maxnoe Aug 1, 2024

mstrzys Aug 5, 2024

maxnoe commented Aug 1, 2024

maxnoe Sep 3, 2024

maxnoe Sep 3, 2024

Dl3 reader for LST data #7

Are you sure you want to change the base?

Dl3 reader for LST data #7

Conversation

mstrzys commented Jun 25, 2023

IevgenVovk commented Jul 14, 2023

mstrzys commented Mar 21, 2024 • edited Loading

maxnoe commented Jul 26, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

maxnoe commented Aug 1, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mstrzys commented Mar 21, 2024 •

edited

Loading

maxnoe commented Jul 26, 2024 •

edited

Loading