Merge pull request #1247 from cta-observatory/update_lst1_location

Update LST1 location from ctapipe_io_lst

lstchain/datachecks/containers.py

@@ -7,11 +7,12 @@
 from astropy import units as u
 from astropy.time import Time
 from astropy.coordinates import SkyCoord, AltAz
-from lstchain.reco.utils import location
+
 import warnings
 from ctapipe.core import Container, Field
 from ctapipe.utils import get_bright_stars
 from ctapipe.coordinates import CameraFrame, TelescopeFrame
+from ctapipe_io_lst.constants import LST1_LOCATION
 
 __all__ = [
     'DL1DataCheckContainer',
@@ -131,9 +132,9 @@ def fill_event_wise_info(self, subrun_index, table, mask, geom,
         self.trigger_type = \
             count_trig_types(table['trigger_type'][mask])
         if 'ucts_jump' in table.columns:
-            # After one (or n) genuine UCTS jumps in a run, the first event (or n events) 
-            # of every subsequent subrun file (if analyzed on its own) will have ucts_jump=True, 
-            # but these are not new jumps, just the ones from previous subruns, so they should 
+            # After one (or n) genuine UCTS jumps in a run, the first event (or n events)
+            # of every subsequent subrun file (if analyzed on its own) will have ucts_jump=True,
+            # but these are not new jumps, just the ones from previous subruns, so they should
             # not be counted.
             uj = table['ucts_jump'].data.copy()
             # find the first False value, and set to False also all the earlier ones:
@@ -156,7 +157,7 @@ def fill_event_wise_info(self, subrun_index, table, mask, geom,
         telescope_pointing = SkyCoord(alt=self.mean_alt_tel,
                                       az=self.mean_az_tel,
                                       frame=AltAz(obstime=time_utc,
-                                                  location=location))
+                                                  location=LST1_LOCATION))
         self.tel_ra = telescope_pointing.icrs.ra.to(u.deg)
         self.tel_dec = telescope_pointing.icrs.dec.to(u.deg)
 
@@ -324,7 +325,7 @@ def fill_pixel_wise_info(self, table, mask, histogram_binnings,
         sampled_times = self.dragon_time
         obstime = Time(sampled_times[int(len(sampled_times)/2)],
                        scale='utc', format='unix')
-        horizon_frame = AltAz(location=location, obstime=obstime)
+        horizon_frame = AltAz(location=LST1_LOCATION, obstime=obstime)
         pointing = SkyCoord(az=self.mean_az_tel,
                             alt=self.mean_alt_tel,
                             frame=horizon_frame)
@@ -336,8 +337,8 @@ def fill_pixel_wise_info(self, table, mask, histogram_binnings,
         camera_frame = CameraFrame(telescope_pointing=pointing,
                                    focal_length=focal_length*relative_shift,
                                    obstime=obstime,
-                                   location=location)
-        telescope_frame = TelescopeFrame(obstime=obstime, location=location)
+                                   location=LST1_LOCATION)
+        telescope_frame = TelescopeFrame(obstime=obstime, location=LST1_LOCATION)
 
         # radius around star within which we consider the pixel may be affected
         # (hence we will later not raise a flag if e.g. its pedestal std dev is
@@ -427,7 +428,7 @@ class DL1DataCheckHistogramBins(Container):
         """
         Histogram bins for the DL1 Datacheck
         """
-        
+
         # delta_t between consecutive events (ms)
         hist_delta_t = Field(
             default_factory=lambda: np.linspace(-1.e-2, 2., 200),

lstchain/high_level/hdu_table.py

@@ -18,7 +18,9 @@
 from astropy.time import Time
 
 from lstchain.__init__ import __version__
-from lstchain.reco.utils import camera_to_altaz, location
+from lstchain.reco.utils import camera_to_altaz
+from ctapipe_io_lst.constants import LST1_LOCATION
+
 
 __all__ = [
     "add_icrs_position_params",
@@ -272,7 +274,7 @@ def get_pointing_params(data, source_pos, time_utc):
     pnt_icrs = SkyCoord(
         alt=pointing_alt[0],
         az=pointing_az[0],
-        frame=AltAz(obstime=time_utc[0], location=location),
+        frame=AltAz(obstime=time_utc[0], location=LST1_LOCATION),
     ).transform_to(frame="icrs")
 
     source_pointing_diff = source_pos.separation(pnt_icrs)
@@ -294,7 +296,7 @@ def add_icrs_position_params(data, source_pos, time_utc):
     reco_az = data["reco_az"]
 
     reco_altaz = SkyCoord(
-        alt=reco_alt, az=reco_az, frame=AltAz(obstime=time_utc, location=location)
+        alt=reco_alt, az=reco_az, frame=AltAz(obstime=time_utc, location=LST1_LOCATION)
     )
 
     with erfa_astrom.set(ErfaAstromInterpolator(300 * u.s)):
@@ -313,14 +315,15 @@ def fill_reco_altaz_w_expected_pos(data):
     for source-dependent analysis. 
 
     Note: This is just a trick to easily extract ON/OFF events in gammapy
-    analysis. For source-dependent analysis, gammaness and alpha cut are 
+    analysis. For source-dependent analysis, gammaness and alpha cut are
     already applied when creating DL3 file and there is no need to apply additional
     cuts in higher analysis (e.g. on/background region cut in gammapy).
     This function fills the same reconstructed position (AltAz, ICRS frame,
-    derived from the first event) for all events. It is recommended to use 
+    derived from the first event) for all events. It is recommended to use
     `WobbleRegionsFinder` in gammapy to define the source and background region.
     """
-    # Compute the expected source position for the first event 
+    # Compute the expected source position for the first event
+
     obstime = Time(data["dragon_time"][0], scale="utc", format="unix")
     expected_src_x = data["expected_src_x"][0] * u.m
     expected_src_y = data["expected_src_y"][0] * u.m
@@ -342,7 +345,7 @@ def fill_reco_altaz_w_expected_pos(data):
 
     reco_altaz = SkyCoord(
         alt=data["reco_alt"][0], az=data["reco_az"][0],
-        frame=AltAz(obstime=obstime, location=location)
+        frame=AltAz(obstime=obstime, location=LST1_LOCATION)
     )
 
     with erfa_astrom.set(ErfaAstromInterpolator(300 * u.s)):
@@ -387,7 +390,7 @@ def create_event_list(
     tel_list = np.unique(data["tel_id"])
 
     time_params, time_utc = get_timing_params(data)
-    
+
     if not 'RA' in data.colnames:
         data = add_icrs_position_params(data, source_pos, time_utc)
     reco_icrs = SkyCoord(ra=data["RA"], dec=data["Dec"], unit="deg")
@@ -469,6 +472,19 @@ def create_event_list(
     ev_header["DEC_OBJ"] = source_pos.dec.to_value()
     ev_header["FOVALIGN"] = "RADEC"
 
+    ev_header["GEOLON"] = (
+        LST1_LOCATION.lon.to_value(u.deg),
+        "Geographic longitude of telescope (deg)",
+    )
+    ev_header["GEOLAT"] = (
+        LST1_LOCATION.lat.to_value(u.deg),
+        "Geographic latitude of telescope (deg)",
+    )
+    ev_header["ALTITUDE"] = (
+        round(LST1_LOCATION.height.to_value(u.m), 2),
+        "Geographic latitude of telescope (m)",
+    )
+
     # GTI table metadata
     gti_header = DEFAULT_HEADER.copy()
     gti_header["CREATED"] = Time.now().utc.iso
@@ -491,20 +507,12 @@ def create_event_list(
     pnt_header["MJDREFF"] = ev_header["MJDREFF"]
     pnt_header["TIMEUNIT"] = ev_header["TIMEUNIT"]
     pnt_header["TIMESYS"] = ev_header["TIMESYS"]
-    pnt_header["OBSGEO-L"] = (
-        location.lon.to_value(u.deg),
-        "Geographic longitude of telescope (deg)",
-    )
-    pnt_header["OBSGEO-B"] = (
-        location.lat.to_value(u.deg),
-        "Geographic latitude of telescope (deg)",
-    )
-    pnt_header["OBSGEO-H"] = (
-        round(location.height.to_value(u.m), 2),
-        "Geographic latitude of telescope (m)",
-    )
-
     pnt_header["TIMEREF"] = ev_header["TIMEREF"]
+
+    pnt_header["GEOLON"] = ev_header["GEOLON"]
+    pnt_header["GEOLAT"] = ev_header["GEOLAT"]
+    pnt_header["ALTITUDE"] = ev_header["ALTITUDE"]
+
     pnt_header["MEAN_ZEN"] = str(data_pars["ZEN_PNT"])
     pnt_header["MEAN_AZ"] = str(data_pars["AZ_PNT"])
     pnt_header["B_DELTA"] = str(data_pars["B_DELTA"])

lstchain/reco/utils.py

@@ -14,10 +14,11 @@
 import astropy.units as u
 import numpy as np
 import pandas as pd
-from astropy.coordinates import AltAz, SkyCoord, EarthLocation
+from astropy.coordinates import AltAz, SkyCoord
 from astropy.time import Time
 from ctapipe.coordinates import CameraFrame
 from ctapipe_io_lst import OPTICS
+from ctapipe_io_lst.constants import LST1_LOCATION
 
 from . import disp
 
@@ -50,10 +51,8 @@
     "get_events_in_GTI"
 ]
 
-# position of the LST1
-location = EarthLocation.from_geodetic(-17.89139 * u.deg, 28.76139 * u.deg, 2184 * u.m)
 obstime = Time("2018-11-01T02:00")
-horizon_frame = AltAz(location=location, obstime=obstime)
+horizon_frame = AltAz(location=LST1_LOCATION, obstime=obstime)
 
 # Geomagnetic parameters for the LST1 as per
 # https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml?#igrfwmm and
@@ -296,7 +295,7 @@ def camera_to_altaz(pos_x, pos_y, focal, pointing_alt, pointing_az, obstime=None
     """
     if not obstime:
         logging.info("No time given. To be use only for MC data.")
-    horizon_frame = AltAz(location=location, obstime=obstime)
+    horizon_frame = AltAz(location=LST1_LOCATION, obstime=obstime)
 
     pointing_direction = SkyCoord(
         alt=clip_alt(pointing_alt), az=pointing_az, frame=horizon_frame
@@ -362,7 +361,7 @@ def radec_to_camera(sky_coordinate, obstime, pointing_alt, pointing_az, focal):
     camera frame: `astropy.coordinates.sky_coordinate.SkyCoord`
     """
 
-    horizon_frame = AltAz(location=location, obstime=obstime)
+    horizon_frame = AltAz(location=LST1_LOCATION, obstime=obstime)
 
     pointing_direction = SkyCoord(
         alt=clip_alt(pointing_alt), az=pointing_az, frame=horizon_frame
@@ -372,7 +371,7 @@ def radec_to_camera(sky_coordinate, obstime, pointing_alt, pointing_az, focal):
         focal_length=focal,
         telescope_pointing=pointing_direction,
         obstime=obstime,
-        location=location,
+        location=LST1_LOCATION,
     )
 
     camera_pos = sky_coordinate.transform_to(camera_frame)
@@ -671,7 +670,7 @@ def get_effective_time(events):
     # might indicate the DAQ was stopped (e.g. if the table contains more
     # than one run). We set 0.01 s as limit to decide a "break" occurred:
     t_elapsed = np.sum(time_diff[time_diff < 0.01 * u.s])
-    
+
     # delta_t is the time elapsed since the previous triggered event.
     # We exclude the null values that might be set for the first even in a file.
     # Same as the elapsed time, we exclude events with delta_t larger than 0.01 s.
@@ -837,7 +836,7 @@ def get_events_in_GTI(events, CatB_cal_table):
 
     Parameters
     ----------
-    events : pandas DataFrame or astropy.table.QTable 
+    events : pandas DataFrame or astropy.table.QTable
         Data frame or table of DL1 or DL2 events.
     CatB_cal_table: table of CatB calibration applied to the events (dl1_mon_tel_CatB_cal_key)
 
@@ -850,4 +849,4 @@ def get_events_in_GTI(events, CatB_cal_table):
 
     gti_mask = gti[events['calibration_id']]
 
-    return events[gti_mask]
+    return events[gti_mask]

lstchain/scripts/lstchain_dl2_add_sourcepos.py

@@ -7,19 +7,19 @@
 from astropy.coordinates import SkyCoord, AltAz
 from astropy.coordinates.erfa_astrom import ErfaAstromInterpolator, erfa_astrom
 from astropy.time import Time
-from lstchain.reco.utils import location
 from ctapipe.coordinates import CameraFrame
 from ctapipe_io_lst import LSTEventSource
+from ctapipe_io_lst.constants import LST1_LOCATION
 
 parser = argparse.ArgumentParser(description="Source position adder. Calculates x,y coordinates of a source in CameraFrame, event-wise.")
 parser.add_argument('-f', '--dl2-file', dest='file', required=True,
                     type=str, help='DL2 file to be processed')
-parser.add_argument('-s', '--source-name', dest='source_name', 
+parser.add_argument('-s', '--source-name', dest='source_name',
                     type=str, default='Crab pulsar',
                     help='Source name (known to astropy; default: %(default)s)')
 
-# NOTE: the position "Crab pulsar" is better defined than "Crab" or 
-# "Crab Nebula", for which astropy.coordinates.SkyCoord.form_name returns 
+# NOTE: the position "Crab pulsar" is better defined than "Crab" or
+# "Crab Nebula", for which astropy.coordinates.SkyCoord.form_name returns
 # different values depending on the catalog it uses...
 
 
@@ -30,12 +30,12 @@ def main():
     tablename = "/dl2/event/telescope/parameters/LST_LSTCam"
     new_table_name = 'source_position'
 
-    # Effective focal length, i.e. including the effect of aberration 
+    # Effective focal length, i.e. including the effect of aberration
     # (which affects all image parameters and hence the reconstructed source
     # position in CameraFrame):
     subarray_info = LSTEventSource.create_subarray(tel_id=1)
     focal = subarray_info.tel[1].optics.effective_focal_length
-    # By using this effective focal the calculated nominal source position 
+    # By using this effective focal the calculated nominal source position
     # will be consistent with the reconstructed event directions in CameraFrame
 
     print('Selected source name:', args.source_name)
@@ -48,17 +48,17 @@ def main():
     pointing_az  = np.array(table['az_tel']) * u.rad
 
     time_utc = Time(table["dragon_time"], format="unix", scale="utc")
-    telescope_pointing = SkyCoord(alt=pointing_alt, az=pointing_az, 
-                                  frame=AltAz(obstime=time_utc, 
-                                              location=location))
+    telescope_pointing = SkyCoord(alt=pointing_alt, az=pointing_az,
+                                  frame=AltAz(obstime=time_utc,
+                                              location=LST1_LOCATION))
 
     # CameraFrame is terribly slow without the erfa interpolator below...
     with erfa_astrom.set(ErfaAstromInterpolator(5 * u.min)):
-        camera_frame = CameraFrame(focal_length=focal, 
-                                   telescope_pointing=telescope_pointing, 
-                                   location=location, obstime=time_utc)
+        camera_frame = CameraFrame(focal_length=focal,
+                                   telescope_pointing=telescope_pointing,
+                                   location=LST1_LOCATION, obstime=time_utc)
         source_pos_camera = source_pos.transform_to(camera_frame)
-    
+
     # Units: m (like reco_src_x, y)
     table['src_x'] = source_pos_camera.data.x.to_value(u.m)
     table['src_y'] = source_pos_camera.data.y.to_value(u.m)