Merge branch 'L2SS-761-add-geo-coords' into 'master'

L2SS-761: Add geo (lat/long) coordinates to antennas to be able to draw them on a map

Closes L2SS-761

See merge request !310

tangostationcontrol/requirements.txt

@@ -13,3 +13,4 @@ docker >= 5.0.3 # Apache 2
 python-logstash-async >= 2.3.0 # MIT
 python-casacore >= 3.3.1 # LGPLv3
 etrs-itrs@git+https://github.com/brentjens/etrs-itrs # license pending
+lofarantpos >= 0.5.0 # Apache 2

tangostationcontrol/tangostationcontrol/beam/geo.py

 import etrsitrs
+import lofarantpos.geo
 import numpy
+import math
 
 """
    LOFAR station positions are measured in ETRS89, which are the coordinates of the position as it would be in 1989.
@@ -16,18 +18,34 @@ import numpy
    The ETRSitrs package does all the transformation calculations for us.
 """
 
+def _apply_fn_on_one_element_or_array(fn, array: numpy.array) -> numpy.array:
+    if array.ndim == 1:
+       # convert a single coordinate triple
+       return fn(array)
+    else:
+       # convert each coordinate triple
+       return numpy.apply_along_axis(fn, 1, array)
+
 def ETRS_to_ITRF(ETRS_coordinates: numpy.array, ITRF_reference_frame: str = "ITRF2005", ITRF_reference_epoch: float = 2015.5) -> numpy.array:
     """ Convert an array of coordinate triples from ETRS to ITRF, in the given reference frame and epoch. """
 
     # fetch converter
     ETRS_to_ITRF_fn = etrsitrs.convert_fn("ETRF2000", ITRF_reference_frame, ITRF_reference_epoch)
 
-    if ETRS_coordinates.ndim == 1:
-       # convert a single coordinate triple
-       ITRF_coordinates = ETRS_to_ITRF_fn(ETRS_coordinates)
-    else:
-       # convert each coordinate triple
-       ITRF_coordinates = numpy.apply_along_axis(ETRS_to_ITRF_fn, 1, ETRS_coordinates)
+    return _apply_fn_on_one_element_or_array(ETRS_to_ITRF_fn, ETRS_coordinates)
+
+def ETRS_to_GEO(ETRS_coordinates: numpy.array) -> numpy.array:
+    """ Convert an array of coordinate triples from ETRS to latitude/longitude (degrees). """
+
+    def ETRS_to_GEO_fn(etrs_coords):
+        geo_coords = lofarantpos.geo.geographic_from_xyz(etrs_coords)
+
+        return numpy.array([
+            geo_coords['lat_rad'] * 180 / math.pi,
+            geo_coords['lon_rad'] * 180 / math.pi
+            ])
+
+    return _apply_fn_on_one_element_or_array(ETRS_to_GEO_fn, ETRS_coordinates)
 
-    # return computed ITRF coordinates
-    return ITRF_coordinates
+# Geo coordinates are only used for rough positioning. The difference between ITRF and ETRS matters little here
+ITRF_to_GEO = ETRS_to_GEO

tangostationcontrol/tangostationcontrol/devices/recv.py

@@ -22,7 +22,7 @@ from math import pi
 
 # Additional import
 from tangostationcontrol.beam.hba_tile import HBATAntennaOffsets
-from tangostationcontrol.beam.geo import ETRS_to_ITRF
+from tangostationcontrol.beam.geo import ETRS_to_ITRF, ITRF_to_GEO
 from tangostationcontrol.common.entrypoint import entry
 from tangostationcontrol.common.lofar_logging import device_logging_to_python
 from tangostationcontrol.clients.attribute_wrapper import attribute_wrapper
@@ -234,6 +234,10 @@ class RECV(opcua_device):
         doc='Absolute reference position of antenna field, in ITRF',
         dtype=(numpy.float,), max_dim_x=3)
 
+    Antenna_Field_Reference_GEO_R = attribute(access=AttrWriteType.READ,
+        doc='Absolute reference position of antenna field, in latitude/longitude (degrees)',
+        dtype=(numpy.float,), max_dim_x=2)
+
     HBAT_antenna_ITRF_offsets_R = attribute(access=AttrWriteType.READ,
         doc='Offsets of the antennas within a tile, in ITRF ("iHBADeltas"). True shape: 96x16x3.',
         dtype=((numpy.float,),), max_dim_x=48, max_dim_y=96)
@@ -242,6 +246,10 @@ class RECV(opcua_device):
         doc='Absolute reference position of each tile, in ITRF',
         dtype=((numpy.float,),), max_dim_x=3, max_dim_y=96)
 
+    HBAT_reference_GEO_R = attribute(access=AttrWriteType.READ,
+        doc='Absolute reference position of each tile, in latitude/longitude (degrees)',
+        dtype=((numpy.float,),), max_dim_x=2, max_dim_y=96)
+
     def read_Antenna_Field_Reference_ITRF_R(self):
         # provide ITRF field coordinates if they were configured
         if self.Antenna_Field_Reference_ITRF:
@@ -251,6 +259,9 @@ class RECV(opcua_device):
         ETRS_coordinates = numpy.array(self.Antenna_Field_Reference_ETRS).reshape(3)
         return ETRS_to_ITRF(ETRS_coordinates, self.ITRF_Reference_Frame, self.ITRF_Reference_Epoch)
 
+    def read_Antenna_Field_Reference_GEO_R(self):
+        return ITRF_to_GEO(self.read_Antenna_Field_Reference_ITRF_R())
+    
     def read_HBAT_antenna_ITRF_offsets_R(self):
         base_antenna_offsets        = numpy.array(self.HBAT_base_antenna_offsets).reshape(16,3)
         PQR_to_ETRS_rotation_matrix = numpy.array(self.HBAT_PQR_to_ETRS_rotation_matrix).reshape(3,3)
@@ -274,6 +285,9 @@ class RECV(opcua_device):
         ETRS_coordinates = numpy.array(self.HBAT_reference_ETRS).reshape(96,3)
         return ETRS_to_ITRF(ETRS_coordinates, self.ITRF_Reference_Frame, self.ITRF_Reference_Epoch)
 
+    def read_HBAT_reference_GEO_R(self):
+        return ITRF_to_GEO(self.read_HBAT_reference_ITRF_R())
+
     # ----------
     # Summarising Attributes
     # ----------

tangostationcontrol/tangostationcontrol/test/beam/test_geo.py

@@ -7,7 +7,7 @@
 # Distributed under the terms of the APACHE license.
 # See LICENSE.txt for more info.
 
-from tangostationcontrol.beam.geo import ETRS_to_ITRF
+from tangostationcontrol.beam.geo import ETRS_to_ITRF, ETRS_to_GEO
 
 from tangostationcontrol.test import base
 
@@ -41,3 +41,32 @@ class TestETRS_to_ITRF(base.TestCase):
         LOFAR1_CS001_LBA_ITRF = [3826923.50275, 460915.488115, 5064643.517]
 
         numpy.testing.assert_almost_equal(CS001_LBA_ITRF, LOFAR1_CS001_LBA_ITRF, decimal=1.5)
+
+class TestETRS_to_GEO(base.TestCase):
+    def test_convert_single_coordinate(self):
+        """ Convert a single coordinate. """
+        ETRS_coords = numpy.array([1.0, 1.0, 1.0])
+        GEO_coords = ETRS_to_GEO(ETRS_coords)
+
+        self.assertEqual((2,), GEO_coords.shape)
+
+    def test_convert_array(self):
+        """ Convert an array of coordinates. """
+        ETRS_coords = numpy.array([ [1.0, 1.0, 1.0], [2.0, 2.0, 2.0], [3.0, 3.0, 3.0] ])
+        GEO_coords = ETRS_to_GEO(ETRS_coords)
+
+        self.assertEqual((3,2), GEO_coords.shape)
+
+    def test_verify_CS001_LBA(self):
+        """ Verify if the calculated CS001LBA phase center matches those calculated in LOFAR1. """
+
+        # See CLBA in MAC/Deployment/data/Coordinates/ETRF_FILES/CS001/CS001-antenna-positions-ETRS.csv 
+        CS001_LBA_ETRS = [3826923.942, 460915.117, 5064643.229]
+
+        # Convert to GEO
+        CS001_LBA_GEO = ETRS_to_GEO(numpy.array(CS001_LBA_ETRS))
+
+        # verify against actual position
+        LOFAR1_CS001_LBA_GEO = [52.911, 6.868]
+
+        numpy.testing.assert_almost_equal(CS001_LBA_GEO, LOFAR1_CS001_LBA_GEO, decimal=3)