Add get_<>_image_area funcs to sarkit.crsd and sarkit.cphd

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - `compute_dwelltimes_using_poly` to `sarkit.cphd`
+- `get_channel_image_area`, `get_extended_image_area`, and `get_scene_image_area` to `sarkit.cphd` and `sarkit.crsd`
 
 ### Removed
 - Unused `_processing` module

sarkit/cphd/__init__.py

@@ -114,6 +114,15 @@
    llh_to_iac
    iac_to_llh
 
+Image Area
+
+.. autosummary::
+   :toctree: generated/
+
+   get_channel_image_area
+   get_extended_image_area
+   get_scene_image_area
+
 Reference Geometry Parameters
 
 .. autosummary::
@@ -201,6 +210,9 @@
 )
 from ._scenecoords import (
     ecf_to_iac,
+    get_channel_image_area,
+    get_extended_image_area,
+    get_scene_image_area,
     hae_iac_to_llh,
     hae_llh_to_iac,
     iac_to_ecf,
@@ -278,8 +290,11 @@
     "dtype_to_binary_format_string",
     "dtype_to_pvp_element",
     "ecf_to_iac",
+    "get_channel_image_area",
     "get_defined_pvp_dtype",
+    "get_extended_image_area",
     "get_pvp_dtype",
+    "get_scene_image_area",
     "hae_iac_to_llh",
     "hae_llh_to_iac",
     "iac_to_ecf",

sarkit/cphd/_scenecoords.py

@@ -295,3 +295,127 @@ def iac_to_llh_from_ew(
             sc_ew["ReferenceSurface"]["HAE"]["uIAYLL"],
         )
     raise ValueError("Could not determine ReferenceSurface")
+
+
+def get_image_area_vertices_from_ew(
+    imgarea_ew: sarkit.xmlhelp.ElementWrapper, *, use_polygon: bool = True
+) -> np.ndarray:
+    """Return the vertices of an ElementWrapped ImageArea.
+
+    Not intended for public API.
+
+    Parameters
+    ----------
+    imgarea_ew : sarkit.xmlhelp.ElementWrapper
+        Element-wrapped ImageArea with children "X1Y1", "X2Y2", and (optionally) "Polygon"
+    use_polygon : bool, optional
+        If ``True``, the polygon, if present, is considered.
+        If ``False``, the polygon is ignored.
+
+    Returns
+    -------
+    ndarray
+        If ``use_polygon == True`` and a polygon is present, its vertices are returned.
+        Otherwise, the four corners of the rectangle described by X1Y1 and X2Y2 are returned.
+    """
+    if use_polygon and (poly := imgarea_ew.get("Polygon", None)) is not None:
+        return poly
+    x1, y1 = imgarea_ew["X1Y1"]
+    x2, y2 = imgarea_ew["X2Y2"]
+    return np.array(
+        [
+            [x1, y1],
+            [x1, y2],
+            [x2, y2],
+            [x2, y1],
+        ]
+    )
+
+
+def get_channel_image_area(
+    cphd_xmltree: lxml.etree.ElementTree, ch_id: str, *, use_polygon: bool = True
+) -> np.ndarray:
+    """Return the vertices of the channel image area identified by ``ch_id``.
+
+    Parameters
+    ----------
+    cphd_xmltree : lxml.etree.ElementTree
+        CPHD XML
+    ch_id : str
+        Channel unique identifier
+    use_polygon : bool, optional
+        If ``True``, the polygon, if present, is considered.
+        If ``False``, the polygon is ignored.
+
+    Returns
+    -------
+    ndarray
+        Vertices of the channel image area in IAC coordinates with IAX, IAY components in meters in the last dimension.
+        If ``use_polygon == True`` and a polygon is present, its vertices are returned.
+        Otherwise, the four corners of the rectangle described by X1Y1 and X2Y2 are returned.
+    """
+    ew = cphd_xml.ElementWrapper(cphd_xmltree.getroot())
+    chan_param_ew = ew["Channel"].find("Parameters", Identifier=ch_id)
+    imgarea_ew = (
+        chan_param_ew["ImageArea"]
+        if "ImageArea" in chan_param_ew
+        else ew["SceneCoordinates"]["ImageArea"]
+    )
+    return get_image_area_vertices_from_ew(imgarea_ew, use_polygon=use_polygon)
+
+
+def get_scene_image_area(
+    cphd_xmltree: lxml.etree.ElementTree, *, use_polygon: bool = True
+) -> np.ndarray:
+    """Return the vertices of the scene image area.
+
+    Parameters
+    ----------
+    cphd_xmltree : lxml.etree.ElementTree
+        CPHD XML
+    use_polygon : bool, optional
+        If ``True``, the polygon, if present, is considered.
+        If ``False``, the polygon is ignored.
+
+    Returns
+    -------
+    ndarray
+        Vertices of the scene image area in IAC coordinates with IAX, IAY components in meters in the last dimension.
+        If ``use_polygon == True`` and a polygon is present, its vertices are returned.
+        Otherwise, the four corners of the rectangle described by X1Y1 and X2Y2 are returned.
+    """
+    ew = cphd_xml.ElementWrapper(cphd_xmltree.getroot())
+    return get_image_area_vertices_from_ew(
+        ew["SceneCoordinates"]["ImageArea"], use_polygon=use_polygon
+    )
+
+
+def get_extended_image_area(
+    cphd_xmltree: lxml.etree.ElementTree, *, use_polygon: bool = True
+) -> np.ndarray | None:
+    """Return the vertices of the extended image area or ``None`` if one is not defined.
+
+    Parameters
+    ----------
+    cphd_xmltree : lxml.etree.ElementTree
+        CPHD XML
+    use_polygon : bool, optional
+        If ``True``, the polygon, if present, is considered.
+        If ``False``, the polygon is ignored.
+
+    Returns
+    -------
+    ndarray or None
+        If an extended image area is defined, the Vertices of the extended image area in IAC coordinates with IAX, IAY
+        components in meters in the last dimension.
+        If ``use_polygon == True`` and a polygon is present, its vertices are returned.
+        Otherwise, the four corners of the rectangle described by X1Y1 and X2Y2 are returned.
+        If an extended image area is not defined, ``None`` is returned.
+    """
+    ew = cphd_xml.ElementWrapper(cphd_xmltree.getroot())
+    imgarea_ew = ew["SceneCoordinates"].get("ExtendedArea", None)
+    return (
+        None
+        if imgarea_ew is None
+        else get_image_area_vertices_from_ew(imgarea_ew, use_polygon=use_polygon)
+    )

sarkit/crsd/__init__.py

@@ -101,6 +101,15 @@
    llh_to_iac
    iac_to_llh
 
+Image Area
+
+.. autosummary::
+   :toctree: generated/
+
+   get_channel_image_area
+   get_extended_image_area
+   get_scene_image_area
+
 Receive Channel Parameters
 ==========================
 
@@ -210,6 +219,9 @@
 )
 from ._scenecoords import (
     ecf_to_iac,
+    get_channel_image_area,
+    get_extended_image_area,
+    get_scene_image_area,
     iac_to_ecf,
     iac_to_llh,
     llh_to_iac,
@@ -292,10 +304,13 @@
     "dtype_to_ppp_element",
     "dtype_to_pvp_element",
     "ecf_to_iac",
+    "get_channel_image_area",
     "get_defined_ppp_dtype",
     "get_defined_pvp_dtype",
+    "get_extended_image_area",
     "get_ppp_dtype",
     "get_pvp_dtype",
+    "get_scene_image_area",
     "iac_to_ecf",
     "iac_to_llh",
     "interpolate_support_array",

sarkit/crsd/_scenecoords.py

@@ -37,6 +37,59 @@ def iac_to_llh(
     return cphd_scenecoords.iac_to_llh_from_ew(sc_ew, pt_iac)
 
 
-for func in (ecf_to_iac, iac_to_ecf, llh_to_iac, iac_to_llh):
+def get_channel_image_area(
+    crsd_xmltree: lxml.etree.ElementTree, ch_id: str, *, use_polygon: bool = True
+) -> np.ndarray:
+    # docstring copied from CPHD version
+    if lxml.etree.QName(crsd_xmltree.getroot()).localname != "CRSDsar":
+        raise ValueError("Only CRSDsar products have channel image areas")
+    ew = crsd_xml.ElementWrapper(crsd_xmltree.getroot())
+    chan_param_ew = ew["Channel"].find("Parameters", Identifier=ch_id)
+    imgarea_ew = (
+        chan_param_ew["SARImage"]["ImageArea"]
+        if "ImageArea" in chan_param_ew["SARImage"]
+        else ew["SceneCoordinates"]["ImageArea"]
+    )
+    return cphd_scenecoords.get_image_area_vertices_from_ew(
+        imgarea_ew, use_polygon=use_polygon
+    )
+
+
+def get_scene_image_area(
+    crsd_xmltree: lxml.etree.ElementTree, *, use_polygon: bool = True
+) -> np.ndarray:
+    # docstring copied from CPHD version
+    ew = crsd_xml.ElementWrapper(crsd_xmltree.getroot())
+    return cphd_scenecoords.get_image_area_vertices_from_ew(
+        ew["SceneCoordinates"]["ImageArea"], use_polygon=use_polygon
+    )
+
+
+def get_extended_image_area(
+    crsd_xmltree: lxml.etree.ElementTree, *, use_polygon: bool = True
+) -> np.ndarray | None:
+    # docstring copied from CPHD version
+    if lxml.etree.QName(crsd_xmltree.getroot()).localname != "CRSDsar":
+        raise ValueError("Only CRSDsar products can have extended image areas")
+    ew = crsd_xml.ElementWrapper(crsd_xmltree.getroot())
+    imgarea_ew = ew["SceneCoordinates"].get("ExtendedArea", None)
+    return (
+        None
+        if imgarea_ew is None
+        else cphd_scenecoords.get_image_area_vertices_from_ew(
+            imgarea_ew, use_polygon=use_polygon
+        )
+    )
+
+
+for func in (
+    ecf_to_iac,
+    iac_to_ecf,
+    llh_to_iac,
+    iac_to_llh,
+    get_channel_image_area,
+    get_scene_image_area,
+    get_extended_image_area,
+):
     newdoc = getattr(getattr(cphd_scenecoords, func.__name__), "__doc__", "")
     func.__doc__ = newdoc.replace("cphd", "crsd").replace("CPHD", "CRSD")

sarkit/verification/_cphd_consistency.py

@@ -16,7 +16,7 @@
 import numpy as np
 import numpy.lib.recfunctions as rfn
 import numpy.polynomial.polynomial as npp
-import shapely.geometry as shg
+import shapely
 from lxml import etree
 
 import sarkit.cphd as skcphd
@@ -326,7 +326,7 @@ def get_polygon(self, polygon_node, check=False):
                 f"{_get_root_path(polygon_node)} size attribute matches the number of vertices"
             ):
                 assert size == len(vertex_nodes)
-            shg_polygon = shg.Polygon(polygon)
+            shg_polygon = shapely.Polygon(polygon)
             with self.need(f"{_get_root_path(polygon_node)} is simple"):
                 assert shg_polygon.is_simple
             with self.need(f"{_get_root_path(polygon_node)} is clockwise"):
@@ -510,19 +510,10 @@ def check_channel_dwell_polys(self, channel_id, channel_node):
         codtime_poly = self.xmlhelp.load_elem(cod_node.find("./{*}CODTimePoly"))
         dwelltime_poly = self.xmlhelp.load_elem(dwell_node.find("./{*}DwellTimePoly"))
 
-        def _get_image_area_polygon(image_area_elem):
-            if image_area_elem.find("./{*}Polygon") is not None:
-                return shg.Polygon(
-                    self.xmlhelp.load_elem(image_area_elem.find("./{*}Polygon"))
-                )
-            x1, y1 = self.xmlhelp.load_elem(image_area_elem.find("./{*}X1Y1"))
-            x2, y2 = self.xmlhelp.load_elem(image_area_elem.find("./{*}X2Y2"))
-            return shg.box(x1, y1, x2, y2)
-
-        image_area_elem = channel_node.find("./{*}ImageArea")
-        if image_area_elem is None:
-            image_area_elem = self.cphdroot.find("./{*}SceneCoordinates/{*}ImageArea")
-        image_area_polygon = _get_image_area_polygon(image_area_elem)
+        ia_vertices = skcphd.get_channel_image_area(
+            self.cphdroot.getroottree(), channel_id
+        )
+        image_area_polygon = shapely.Polygon(ia_vertices)
 
         def _get_points_in_polygon(polygon, grid_size=25):
             bounds = np.asarray(polygon.bounds).reshape(
@@ -535,7 +526,7 @@ def _get_points_in_polygon(polygon, grid_size=25):
                 ),
                 axis=-1,
             )
-            coords = shg.MultiPoint(
+            coords = shapely.MultiPoint(
                 np.concatenate(
                     [mesh.reshape(-1, 2), np.asarray(polygon.exterior.coords)[:-1, :]],
                     axis=0,
@@ -554,17 +545,7 @@ def _get_points_in_polygon(polygon, grid_size=25):
         with self.precondition():
             pvp = self._get_channel_pvps(channel_id)
             mask = np.isfinite(pvp["TxTime"])
-
-            def calc_tref(v):
-                r_xmt = np.linalg.norm(v["TxPos"] - v["SRPPos"])
-                r_rcv = np.linalg.norm(v["RcvPos"] - v["SRPPos"])
-                return v["TxTime"] + r_xmt / (r_xmt + r_rcv) * (
-                    v["RcvTime"] - v["TxTime"]
-                )
-
-            pvps_tref1 = calc_tref(pvp[mask][0])
-            pvps_tref2 = calc_tref(pvp[mask][-1])
-
+            pvps_tref1, pvps_tref2 = skcphd.compute_t_ref_from_pvps(pvp[mask][[0, -1]])
             with self.need(
                 "/Dwell/CODTime/CODTimePoly and /Dwell/DwellTime/DwellTimePoly supported by PVPs"
             ):
@@ -1539,7 +1520,7 @@ def check_image_area_corner_points(self):
             assert [int(x.attrib["index"]) for x in vertex_nodes] == list(
                 range(1, len(vertex_nodes) + 1)
             )
-            shg_polygon = shg.Polygon(polygon)
+            shg_polygon = shapely.Polygon(polygon)
             with self.need("Polygon is simple"):
                 assert shg_polygon.is_simple
             with self.need("Polygon is clockwise"):
@@ -1576,10 +1557,10 @@ def check_extended_imagearea_polygon(self):
                 with self.precondition():
                     assert polygon_node is not None
                     polygon = self.get_polygon(polygon_node)
-                    shg_extended = shg.Polygon(extended_area_polygon)
-                    shg_polygon = shg.Polygon(polygon)
+                    shg_extended = shapely.Polygon(extended_area_polygon)
+                    shg_polygon = shapely.Polygon(polygon)
                     with self.need("Extended area polygon covers image area polygon"):
-                        assert shg.Polygon(shg_extended).covers(shg_polygon)
+                        assert shg_extended.covers(shg_polygon)
 
     @per_channel
     def check_channel_imagearea_x1y1(self, channel_id, channel_node):