[WIP] SML plotting code

src/poligrain/plot_map.py

@@ -2,6 +2,8 @@
 
 from __future__ import annotations
 
+import math
+
 import matplotlib.axes
 import matplotlib.patheffects as pe
 import matplotlib.pyplot as plt
@@ -214,3 +216,89 @@ def plot_lines(
         vmax=vmax,
         cmap=cmap,
     )
+
+
+def calculate_azimuth(
+    site_0_lon: npt.ArrayLike | float,
+    site_1_lat: npt.ArrayLike | float,
+    satellite_lon: npt.ArrayLike | float,
+) -> float:
+    """Calculate SML azimuth angel.
+
+    Calculates the azimuth based on the location of a sensor and the
+    longitude of a geostationary satellite.
+
+    Parameters
+    ----------
+       site_0_lon : npt.ArrayLike | float
+           Longitude of sensor.
+       site_1_lat : npt.ArrayLike | float
+           Latitude of sensor.
+       satellite_lon : npt.ArrayLike | float
+           Longitude of geostationary satellite.
+
+    Returns
+    -------
+       float
+
+    """
+    # Convert degrees to radians
+    dish_lat_rad = math.radians(site_1_lat)
+    dish_lon_rad = math.radians(site_0_lon)
+    satellite_lon_rad = math.radians(satellite_lon)
+
+    # Calculate the difference in longitude
+    delta_lon = satellite_lon_rad - dish_lon_rad
+
+    # Calculate the azimuth angle
+    x = math.sin(delta_lon)
+    y = math.cos(dish_lat_rad) * math.tan(math.radians(0)) - math.sin(
+        dish_lat_rad
+    ) * math.cos(delta_lon)
+
+    azimuth_rad = math.atan2(x, y)
+
+    # Convert from radians to degrees and normalize to 0-360
+    azimuth_deg = math.degrees(azimuth_rad)
+    return (azimuth_deg + 360) % 360
+
+
+def calc_sml_line_on_ground(
+    x: npt.ArrayLike,
+    y: npt.ArrayLike,
+    elevation_angel: npt.ArrayLike,
+    azimuth_angel: npt.ArrayLike,
+    melting_layer_height: float,
+):
+    """Calculate the projection of one or many SML on the ground for plotting.
+
+    Parameters
+    ----------
+    x : npt.ArrayLike | float
+        x coordinate of the location of the SML antenna in projected coordinates
+        with the unit meters.
+    y : npt.ArrayLike | float
+        y coordinate of the location of the SML antenna in projected coordinates
+        with the unit meters.
+    elevation_angel: npt.ArrayLike | float
+        Antenna angel in degrees.
+    azimuth_angel: npt.ArrayLike | float
+        Azimuth angel in degrees.
+    melting_layer_height: float
+        Height of melting layer above sensor in meter.
+
+    Returns
+    -------
+    Tuple of the antenna coordinates and the coordinates of the SML beam
+    projected to the ground based on the melting layer height.
+
+    """
+    length_on_ground = (
+        melting_layer_height
+        * np.sin(np.radians(90))
+        / np.sin(np.radians(elevation_angel))
+    )
+    endy = y + length_on_ground * np.cos(np.radians(azimuth_angel))
+    endx = x + length_on_ground * np.sin(np.radians(azimuth_angel))
+
+    return (x, y, endx, endy)

tests/test_plot_map.py

@@ -158,3 +158,17 @@ def test_plot_lines_with_dataarray_raise_wrong_shape():
     ds_cmls = xr.open_dataset("tests/test_data/openMRG_CML_180minutes.nc")
     with pytest.raises(ValueError, match="has to be 1D"):
         plg.plot_map.plot_lines(ds_cmls.rsl.isel(sublink_id=0))
+
+
+def test_calc_sml_line_on_ground():
+    x0, y0, x1, y1 = plg.plot_map.calc_sml_line_on_ground(
+        x=7.576159e05,
+        y=6.454985e06,
+        azimuth_angel=159.862,
+        elevation_angel=36.950,
+        melting_layer_height=4000,
+    )
+
+    np.testing.assert_almost_equal(
+        (x0, y0, x1, y1), (757615.9, 6454985.0, 759906.8480931921, 6448737.532106993)
+    )