Repository: cta-observatory/ctapipe
Pull request #2167 (2022-12-14 14:28:53)

Implement transformation from GroundFrame to EarthLocation


Patch:
diff --git a/ctapipe/coordinates/ground_frames.py b/ctapipe/coordinates/ground_frames.py
index 079267f4..634c2fc6 100644
--- a/ctapipe/coordinates/ground_frames.py
+++ b/ctapipe/coordinates/ground_frames.py
@@ -16,11 +16,13 @@ TODO:
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import (
+    ITRS,
     AffineTransform,
     AltAz,
     BaseCoordinateFrame,
     CartesianRepresentation,
     CoordinateAttribute,
+    EarthLocationAttribute,
     FunctionTransform,
     RepresentationMapping,
     frame_transform_graph,
@@ -35,6 +37,23 @@ __all__ = [
 ]
 
 
+def _altaz_to_earthlocation(altaz):
+    local_itrs = altaz.transform_to(ITRS(location=altaz.location))
+    itrs = ITRS(local_itrs.cartesian + altaz.location.get_itrs().cartesian)
+    return itrs.earth_location
+
+
+def _earthlocation_to_altaz(location, reference_location):
+    # See
+    # https://docs.astropy.org/en/stable/coordinates/common_errors.html#altaz-calculations-for-earth-based-objects
+    # for why this is necessary and we cannot just do
+    # `get_itrs().transform_to(AltAz())`
+    itrs_cart = location.get_itrs().cartesian
+    itrs_ref_cart = reference_location.get_itrs().cartesian
+    local_itrs = ITRS(itrs_cart - itrs_ref_cart, location=reference_location)
+    return local_itrs.transform_to(AltAz(location=reference_location))
+
+
 class GroundFrame(BaseCoordinateFrame):
     """Ground coordinate frame.  The ground coordinate frame is a simple
     cartesian frame describing the 3 dimensional position of objects
@@ -48,6 +67,44 @@ class GroundFrame(BaseCoordinateFrame):
     """
 
     default_representation = CartesianRepresentation
+    reference_location = EarthLocationAttribute()
+
+    def to_earth_location(self):
+        """Convert this GroundFrame coordinate into an `astropy.coordinates.EarthLocation`
+
+        This requires that the ``reference_location`` is set.
+        """
+        # in astropy, x points north, y points east, so we need a minus for y.
+        cart = CartesianRepresentation(self.x, -self.y, self.z)
+        altaz = AltAz(cart, location=self.reference_location)
+        return _altaz_to_earthlocation(altaz)
+
+    @classmethod
+    def from_earth_location(cls, location, reference_location):
+        """
+        Convert `~astropy.coordinates.EarthLocation` into Groundframe.
+
+        Parameters
+        ----------
+        location : astropy.coordinates.EarthLocation
+            The location to convert
+        reference_location : astropy.coordinates.EarthLocation
+            Reference location for the GroundFrame
+
+        Returns
+        -------
+        ground_frame : GroundFrame
+            EarthLocation converted to GroundFrame
+        """
+        altaz = _earthlocation_to_altaz(location, reference_location)
+        x, y, z = altaz.cartesian.xyz
+        # in astropy, x points north, y points east, so we need a minus for y.
+        return GroundFrame(x=x, y=-y, z=z, reference_location=reference_location)
+
+    @property
+    def observation_level(self):
+        """Height of the reference location"""
+        return self.reference_location.height
 
 
 class EastingNorthingFrame(BaseCoordinateFrame):
@@ -59,6 +116,7 @@ class EastingNorthingFrame(BaseCoordinateFrame):
     """
 
     default_representation = CartesianRepresentation
+    reference_location = EarthLocationAttribute()
 
     frame_specific_representation_info = {
         CartesianRepresentation: [
@@ -68,6 +126,28 @@ class EastingNorthingFrame(BaseCoordinateFrame):
         ]
     }
 
+    def to_earth_location(self):
+        # in astropy, x points north, y points east
+        cart = CartesianRepresentation(self.northing, self.easting, self.height)
+        altaz = AltAz(cart, location=self.reference_location)
+        return _altaz_to_earthlocation(altaz)
+
+    @classmethod
+    def from_earth_location(cls, location, reference_location):
+        altaz = _earthlocation_to_altaz(location, reference_location)
+        x, y, z = altaz.cartesian.xyz
+        # in astropy, x points north, y points east
+        return GroundFrame(
+            northing=x,
+            easting=y,
+            height=z,
+            reference_location=reference_location,
+        )
+
+    @property
+    def observation_level(self):
+        return self.reference_location.height
+
 
 class TiltedGroundFrame(BaseCoordinateFrame):
     """Tilted ground coordinate frame.  The tilted ground coordinate frame
diff --git a/docs/changes/2167.feature.rst b/docs/changes/2167.feature.rst
new file mode 100644
index 00000000..fc64d4ec
--- /dev/null
+++ b/docs/changes/2167.feature.rst
@@ -0,0 +1,4 @@
+It's now possible to transform between ``GroundFrame`` coordinates
+and ``astropy.coordinates.EarthLocation``, enabling the conversion
+between relative array coordinates (used in the simulation) and
+absolute real-world coordinates.


Test patch:
diff --git a/ctapipe/coordinates/tests/test_ground_frame.py b/ctapipe/coordinates/tests/test_ground_frame.py
new file mode 100644
index 00000000..6c5c9ab2
--- /dev/null
+++ b/ctapipe/coordinates/tests/test_ground_frame.py
@@ -0,0 +1,57 @@
+"""
+
+"""
+import astropy.units as u
+import pytest
+from astropy.coordinates import EarthLocation
+
+# LST-1, MAGIC-1, MAGIC-2
+locations = [
+    EarthLocation(lon=-17.89149701 * u.deg, lat=28.76152611 * u.deg, height=2185 * u.m),
+    EarthLocation(
+        lon=-17.8900665372 * u.deg, lat=28.7619439944 * u.deg, height=2184.763 * u.m
+    ),
+    EarthLocation(
+        lon=-17.8905595556 * u.deg, lat=28.7613094808 * u.deg, height=2186.764 * u.m
+    ),
+]
+
+# dish-area weighted center
+reference_location = EarthLocation(
+    lon=-17.890879 * u.deg,
+    lat=28.761579 * u.deg,
+    height=2199 * u.m,
+)
+
+# these are the current positions in the current simtel-config,
+# they used a slightly different inputs, so results are not expected to be
+# exactly equal
+mc_positions = [
+    [-6.336, 60.405, 12.5],
+    [41.054, -79.275, 10.75],
+    [-29.456001, -31.295, 11.92],
+] * u.m
+
+
+@pytest.mark.parametrize("location,expected", zip(locations, mc_positions))
+def test_ground_frame_to_earth_location(location, expected):
+    from ctapipe.coordinates import GroundFrame
+
+    ground_frame = GroundFrame.from_earth_location(location, reference_location)
+    x, y, z = ground_frame.cartesian.xyz
+    ex, ey, _ = expected
+
+    # the calculation for the mc position used slightly different inputs
+    # so we accept a bit of a tolerance here
+    assert u.isclose(x, ex, atol=0.6 * u.m)
+    assert u.isclose(y, ey, atol=0.6 * u.m)
+
+    # the mc positions are that of the center of the dish and were shifted
+    # by simtel, so we compare our expectation here and not to the simtel values
+    assert u.isclose(z, location.height - reference_location.height, atol=0.1 * u.m)
+
+    # test back_trafo
+    back = ground_frame.to_earth_location()
+    assert u.isclose(back.lon, location.lon)
+    assert u.isclose(back.lat, location.lat)
+    assert u.isclose(back.height, location.height)