diff --git a/tangostationcontrol/requirements.txt b/tangostationcontrol/requirements.txt
index 62cf05034ea7b40c479f6b75eb4c8d2641f2dcf9..b65be92168032276a4120804d581d4bc95e6c028 100644
--- a/tangostationcontrol/requirements.txt
+++ b/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
diff --git a/tangostationcontrol/tangostationcontrol/beam/geo.py b/tangostationcontrol/tangostationcontrol/beam/geo.py
index 033a6c4e4293573d05cb3fa09bcab5071c88a97b..2c8195796811237750becb29f28fcf904f34f03d 100644
--- a/tangostationcontrol/tangostationcontrol/beam/geo.py
+++ b/tangostationcontrol/tangostationcontrol/beam/geo.py
@@ -1,5 +1,7 @@
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
diff --git a/tangostationcontrol/tangostationcontrol/devices/recv.py b/tangostationcontrol/tangostationcontrol/devices/recv.py
index 93233f86f4a0dd0291205dea5f30e1ab9ae8c5fc..1254fbdc35e4975751c34836da248778a3468cd0 100644
--- a/tangostationcontrol/tangostationcontrol/devices/recv.py
+++ b/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:
@@ -250,6 +258,9 @@ class RECV(opcua_device):
# calculate them from ETRS coordinates if not, using the configured ITRF reference
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)
@@ -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
# ----------
diff --git a/tangostationcontrol/tangostationcontrol/test/beam/test_geo.py b/tangostationcontrol/tangostationcontrol/test/beam/test_geo.py
index 858b3f32e954d19271f8a0dc6fc3cba7b92f47e2..5694376be684b96e74a6197e5e8c3d13f3df6d39 100644
--- a/tangostationcontrol/tangostationcontrol/test/beam/test_geo.py
+++ b/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)