Use SkyGrid functions to write the HDF5 file and fix error in angular…

… convention (#5000)

* Using SkyGrid functions to write the HDF5 file and calculate time
delays

* Adding requested changes and changing the name of the instrument
attribute into detectors (doesn't work)

* Renaming the instruments argument as instruments and not detectors
after discussion

* Fixing indentations and name of variables, removing the calculation of
the time delays for now

* Adding extraline after imports

* Modifying the the conversion of cartesian coordinates into spherical
ones

* Add sanity checks on RA and DEC

* Correcting checks on ra and dec and ensure that they are np arrays

* Changing the way of importing SkyGrid class and correcting the import
of pycbc.conversions

---------

Co-authored-by: jacquot <jacquot@>

bin/pycbc_make_sky_grid

@@ -18,8 +18,8 @@ from scipy.spatial.transform import Rotation as R
 
 import pycbc
 from pycbc.detector import Detector
-from pycbc.io.hdf import HFile
 from pycbc.types import angle_as_radians
+from pycbc.tmpltbank.sky_grid import SkyGrid
 
 
 def spher_to_cart(sky_points):
@@ -36,6 +36,7 @@ def cart_to_spher(sky_points):
     spher = np.zeros((len(sky_points), 2))
     spher[:, 0] = np.arctan2(sky_points[:, 1], sky_points[:, 0])
     spher[:, 1] = np.arcsin(sky_points[:, 2])
+    spher[spher[:, 0] <0, 0] += 2 * np.pi
     return spher
 
 
@@ -163,25 +164,15 @@ rota = r.apply(cart)
 # Convert cartesian coordinates back to spherical coordinates
 spher = cart_to_spher(rota)
 
-# Calculate the time delays between the Earth center
-# and each detector for each sky point
-time_delays = [
-    [
-        detectors[i].time_delay_from_earth_center(
-            spher[j][0], spher[j][1], args.trigger_time
-        )
-        for j in range(len(spher))
-    ]
-    for i in range(len(detectors))
-]
+sky_grid = SkyGrid(
+    spher[:, 0], spher[:, 1], args.instruments, args.trigger_time
+)
+
+extra_attributes = {
+    'trigger_ra': args.ra,
+    'trigger_dec': args.dec,
+    'sky_error': args.sky_error,
+    'timing_uncertainty': args.timing_uncertainty
+}
 
-with HFile(args.output, 'w') as hf:
-    hf['ra'] = spher[:, 0]
-    hf['dec'] = spher[:, 1]
-    hf['trigger_ra'] = [args.ra]
-    hf['trigger_dec'] = [args.dec]
-    hf['sky_error'] = [args.sky_error]
-    hf['trigger_time'] = [args.trigger_time]
-    hf['timing_uncertainty'] = [args.timing_uncertainty]
-    hf['instruments'] = [d for d in args.instruments]
-    hf['time_delays'] = time_delays
+sky_grid.write_to_file(args.output, extra_attributes)

pycbc/tmpltbank/sky_grid.py

@@ -8,6 +8,7 @@
 import h5py
 
 from pycbc.detector import Detector
+from pycbc.conversions import ensurearray
 
 
 class SkyGrid:
@@ -34,6 +35,11 @@ def __init__(self, ra, dec, detectors, ref_gps_time):
         # We store the points in a 2D array internally, first dimension runs
         # over the list of points, second dimension is RA/dec.
         # Question: should we use Astropy sky positions instead?
+        ra, dec, _ = ensurearray(ra,dec)
+        if (ra < 0).any() or (ra > 2 * np.pi).any():
+            raise ValueError('RA must be in the range [0,2π]')
+        if (dec < -np.pi/2).any() or (dec > np.pi/2).any():
+            raise ValueError('DEC must be in the range [-π/2, π/2]')
         self.positions = np.vstack([ra, dec]).T
         self.detectors = sorted(detectors)
         self.ref_gps_time = ref_gps_time
@@ -91,13 +97,17 @@ def read_from_file(cls, path):
             ref_gps_time = hf.attrs['ref_gps_time']
         return cls(ra, dec, detectors, ref_gps_time)
 
-    def write_to_file(self, path):
+    def write_to_file(self, path, extra_attrs=None, extra_datasets=None):
         """Writes a sky grid to an HDF5 file."""
         with h5py.File(path, 'w') as hf:
             hf['ra'] = self.ras
             hf['dec'] = self.decs
             hf.attrs['detectors'] = self.detectors
             hf.attrs['ref_gps_time'] = self.ref_gps_time
+            for attribute in (extra_attrs or {}):
+                hf.attrs[attribute] = extra_attrs[attribute]
+            for dataset in (extra_datasets or {}):
+                hf[dataset] = extra_datasets[dataset]
 
     def calculate_antenna_patterns(self):
         """Calculate the antenna pattern functions at each point in the grid