From 2168535dd77c480fcc811bf0438b1b7a3cec6f33 Mon Sep 17 00:00:00 2001 From: Daniel Pressler Date: Thu, 6 Aug 2026 16:40:39 +0000 Subject: [PATCH 1/2] Add get_<>_image_area funcs to sarkit.cphd --- CHANGELOG.md | 3 + sarkit/cphd/__init__.py | 15 +++ sarkit/cphd/_scenecoords.py | 124 +++++++++++++++++++++++ sarkit/verification/_cphd_consistency.py | 43 +++----- tests/core/cphd/test_scenecoords.py | 70 +++++++++++++ 5 files changed, 224 insertions(+), 31 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e11385b..1e7b11e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- `get_channel_image_area`, `get_extended_image_area`, and `get_scene_image_area` to `sarkit.cphd` + ### Removed - Unused `_processing` module diff --git a/sarkit/cphd/__init__.py b/sarkit/cphd/__init__.py index 476d65d..89043f9 100644 --- a/sarkit/cphd/__init__.py +++ b/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:: @@ -192,6 +201,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, @@ -268,8 +280,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", diff --git a/sarkit/cphd/_scenecoords.py b/sarkit/cphd/_scenecoords.py index e2c2b88..679802e 100644 --- a/sarkit/cphd/_scenecoords.py +++ b/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) + ) diff --git a/sarkit/verification/_cphd_consistency.py b/sarkit/verification/_cphd_consistency.py index e560502..312f117 100644 --- a/sarkit/verification/_cphd_consistency.py +++ b/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): diff --git a/tests/core/cphd/test_scenecoords.py b/tests/core/cphd/test_scenecoords.py index 6a7a307..7cb33cb 100644 --- a/tests/core/cphd/test_scenecoords.py +++ b/tests/core/cphd/test_scenecoords.py @@ -1,8 +1,10 @@ +import copy import pathlib import lxml.etree import numpy as np import pytest +import shapely import sarkit.cphd as skcphd import sarkit.wgs84 @@ -115,3 +117,71 @@ def test_derived_tofrom_iac(surf_type, cphd_xmltree_func): skcphd.iac_to_llh(cphd_xmltree, pt_iacs[..., :2]), skcphd.iac_to_llh(cphd_xmltree, pt_iacs * [1, 1, 0]), ) + + +@pytest.mark.parametrize("has_polygon", (True, False)) +@pytest.mark.parametrize("use_polygon", (True, False)) +def test_image_area_funcs(has_polygon, use_polygon): + xmltree = lxml.etree.parse(DATAPATH / "example-cphd-1.1.0.xml") + + extended_poly = shapely.Polygon([[-10, -10], [0, 20], [10, -10]]) + scene_poly = shapely.buffer(extended_poly, -1) + ch0_poly = shapely.buffer(scene_poly, -1) + ch1_poly = shapely.buffer(ch0_poly, -1) + + def set_ia(ia_ew: skcphd.ElementWrapper, polygon: shapely.Polygon): + ia_ew["X1Y1"] = [polygon.bounds[0], polygon.bounds[1]] + ia_ew["X2Y2"] = [polygon.bounds[2], polygon.bounds[3]] + del ia_ew["Polygon"] + if has_polygon: + ia_ew["Polygon"] = shapely.get_coordinates(polygon)[:-1, :] + + ew = skcphd.ElementWrapper(xmltree.getroot()) + set_ia(ew["SceneCoordinates"]["ExtendedArea"], extended_poly) + set_ia(ew["SceneCoordinates"]["ImageArea"], scene_poly) + chpar0 = ew["Channel"]["Parameters"][0] + set_ia(chpar0["ImageArea"], ch0_poly) + chpar1 = copy.deepcopy(chpar0) + chpar1["Identifier"] = "chpar1_id" # assume this is unique + set_ia(chpar1["ImageArea"], ch1_poly) + ew["Channel"].add("Parameters", chpar1) + chpar2 = copy.deepcopy(chpar0) + chpar2["Identifier"] = "chpar2_id" # assume this is unique + del chpar2["ImageArea"] + ew["Channel"].add("Parameters", chpar2) + + def check_imgarea(actual, expected): + if has_polygon and use_polygon: + assert shapely.equals(shapely.Polygon(actual), expected) + else: + assert shapely.equals(shapely.Polygon(actual), expected.envelope) + + check_imgarea( + skcphd.get_channel_image_area( + xmltree, chpar0["Identifier"], use_polygon=use_polygon + ), + ch0_poly, + ) + check_imgarea( + skcphd.get_channel_image_area( + xmltree, chpar1["Identifier"], use_polygon=use_polygon + ), + ch1_poly, + ) + check_imgarea( + skcphd.get_channel_image_area( + xmltree, chpar2["Identifier"], use_polygon=use_polygon + ), + scene_poly, + ) + + check_imgarea( + skcphd.get_scene_image_area(xmltree, use_polygon=use_polygon), scene_poly + ) + + check_imgarea( + skcphd.get_extended_image_area(xmltree, use_polygon=use_polygon), extended_poly + ) + + del ew["SceneCoordinates"]["ExtendedArea"] + assert skcphd.get_extended_image_area(xmltree, use_polygon=use_polygon) is None From 865026dd7696be97189a80f41986ab7fb8c50767 Mon Sep 17 00:00:00 2001 From: Daniel Pressler Date: Fri, 7 Aug 2026 14:18:12 -0700 Subject: [PATCH 2/2] Add get_<>_image_area funcs to sarkit.crsd --- CHANGELOG.md | 2 +- sarkit/crsd/__init__.py | 15 ++++ sarkit/crsd/_scenecoords.py | 55 +++++++++++++- tests/conftest.py | 104 +------------------------- tests/core/crsd/test_scenecoords.py | 80 ++++++++++++++++++++ tests/utils.py | 112 +++++++++++++++++++++++++++- 6 files changed, 265 insertions(+), 103 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e7b11e..508c48f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added -- `get_channel_image_area`, `get_extended_image_area`, and `get_scene_image_area` 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 diff --git a/sarkit/crsd/__init__.py b/sarkit/crsd/__init__.py index 912d986..ffa4f2d 100644 --- a/sarkit/crsd/__init__.py +++ b/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", diff --git a/sarkit/crsd/_scenecoords.py b/sarkit/crsd/_scenecoords.py index 4fbd792..0b7f55f 100644 --- a/sarkit/crsd/_scenecoords.py +++ b/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") diff --git a/tests/conftest.py b/tests/conftest.py index aac3690..950cb9c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,4 +1,3 @@ -import copy import pathlib import numpy as np @@ -12,6 +11,8 @@ import sarkit.sidd as sksidd from sarkit import _constants +from . import utils + DATAPATH = pathlib.Path(__file__).parents[1] / "data" good_cphd_xml_path = DATAPATH / "example-cphd-1.1.0.xml" @@ -426,71 +427,13 @@ def compute_vh_pol(pos, acx, acy, ref_pt, ant_pol_ref, txrcv_pol_ref): yield tmp_crsd -def _remove(root, pattern): - if (elem := root.find(pattern)) is not None: - elem.getparent().remove(elem) - else: - print(f"Cannot find {pattern=}") - - -def _replace_error(crsd_etree, sensor_type): - sar_error = crsd_etree.find("{*}ErrorParameters/{*}SARImage") - elem_ns = etree.QName(sar_error).namespace - retval = copy.deepcopy(sar_error.find("{*}Monostatic")) - retval.tag = f"{{{elem_ns}}}{sensor_type}Sensor" - sar_error.addnext(retval) - helper = skcrsd.XmlHelper(crsd_etree) - ndx = {"Tx": 0, "Rcv": 1}[sensor_type] - helper.set_elem( - retval.find(".//{*}TimeFreqCov"), - skcrsd.MtxType((3, 3)).parse_elem(retval.find(".//{*}TimeFreqCov"))[ - [ndx, 2], : - ][:, [ndx, 2]], - ) - time_decorr = copy.deepcopy(retval.find(f".//{{*}}{{{sensor_type}}}TimeDecorr")) - if time_decorr is not None: - time_decorr.tag = f"{{{elem_ns}}}TimeDecorr" - _remove(retval, "{*}TxTimeDecorr") - _remove(retval, "{*}RcvTimeDecorr") - retval.find(".//{*}ClockFreqDecorr").addprevious(time_decorr) - sar_error.getparent().remove(sar_error) - - -def _repack_support_arrays(crsd_etree): - offset = 0 - for array in crsd_etree.findall("{*}Data/{*}Support/{*}SupportArray"): - array.find("{*}ArrayByteOffset").text = str(offset) - offset += ( - int(array.findtext("{*}NumRows")) - * int(array.findtext("{*}NumCols")) - * int(array.findtext("{*}BytesPerElement")) - ) - return offset - - @pytest.fixture(scope="session") def example_crsdtx(tmp_path_factory, example_crsdsar): with example_crsdsar.open("rb") as f, skcrsd.Reader(f) as cr: crsd_etree = cr.metadata.xmltree sequence_id = crsd_etree.findtext("{*}TxSequence/{*}Parameters/{*}Identifier") ppps = cr.read_ppps(sequence_id) - crsd_etree.find(".//{*}RefPulseIndex").text = crsd_etree.find( - ".//{*}RefVectorPulseIndex" - ).text - ns = etree.QName(crsd_etree.getroot()).namespace - crsd_etree.getroot().tag = f"{{{ns}}}CRSDtx" - _remove(crsd_etree, "{*}SARInfo") - _remove(crsd_etree, "{*}ReceiveInfo") - _remove(crsd_etree, "{*}Global/{*}Receive") - _remove(crsd_etree, "{*}SceneCoordinates/{*}ExtendedArea") - _remove(crsd_etree, "{*}SceneCoordinates/{*}ImageGrid") - _remove(crsd_etree, "{*}Data/{*}Receive") - _remove(crsd_etree, "{*}Channel") - _remove(crsd_etree, "{*}ReferenceGeometry/{*}SARImage") - _remove(crsd_etree, "{*}ReferenceGeometry/{*}RcvParameters") - _remove(crsd_etree, "{*}DwellPolynomials") - _remove(crsd_etree, "{*}PVP") - _replace_error(crsd_etree, "Tx") + utils.crsdsar_xml_to_crsdtx(crsd_etree) tmp_crsd = ( tmp_path_factory.mktemp("data") / good_crsd_xml_path.with_suffix(".crsd").name ) @@ -510,46 +453,7 @@ def example_crsdrcv(tmp_path_factory, example_crsdsar): channel_id = crsd_etree.findtext("{*}Channel/{*}Parameters/{*}Identifier") pvps = cr.read_pvps(channel_id) signal = cr.read_signal(channel_id) - ns = etree.QName(crsd_etree.getroot()).namespace - crsd_etree.getroot().tag = f"{{{ns}}}CRSDrcv" - _remove(crsd_etree, "{*}SARInfo") - _remove(crsd_etree, "{*}TransmitInfo") - _remove(crsd_etree, "{*}Global/{*}Transmit") - _remove(crsd_etree, "{*}SceneCoordinates/{*}ExtendedArea") - _remove(crsd_etree, "{*}SceneCoordinates/{*}ImageGrid") - _remove(crsd_etree, "{*}Data/{*}Transmit") - _remove(crsd_etree, "{*}TxSequence") - _remove(crsd_etree, "{*}Channel/{*}Parameters/{*}SARImage") - _remove(crsd_etree, "{*}ReferenceGeometry/{*}SARImage") - _remove(crsd_etree, "{*}ReferenceGeometry/{*}TxParameters") - _remove(crsd_etree, "{*}DwellPolynomials") - fx_ids = [ - x.text - for x in crsd_etree.findall("{*}SupportArray/{*}FxResponseArray/{*}Identifier") - ] - xm_ids = [ - x.text for x in crsd_etree.findall("{*}SupportArray/{*}XMArray/{*}Identifier") - ] - _remove(crsd_etree, "{*}SupportArray/{*}FxResponseArray") - _remove(crsd_etree, "{*}SupportArray/{*}XMArray") - for x in fx_ids + xm_ids: - _remove( - crsd_etree, - f"{{*}}Data/{{*}}Support/{{*}}SupportArray[{{*}}SAId='{x}']", - ) - nsa = crsd_etree.find("{*}Data/{*}Support/{*}NumSupportArrays") - nsa.text = str(int(nsa.text) - len(fx_ids + xm_ids)) - _repack_support_arrays(crsd_etree) - _remove(crsd_etree, "{*}PPP") - tx_pulse_index_offset = int(crsd_etree.findtext("{*}PVP/{*}TxPulseIndex/{*}Offset")) - _remove(crsd_etree, "{*}PVP/{*}TxPulseIndex") - for pvp_offset in crsd_etree.findall("{*}PVP/*/{*}Offset"): - if int(pvp_offset.text) > tx_pulse_index_offset: - pvp_offset.text = str(int(pvp_offset.text) - 1) - crsd_etree.find("{*}Data/{*}Receive/{*}NumBytesPVP").text = str( - int(crsd_etree.findtext("{*}Data/{*}Receive/{*}NumBytesPVP")) - 8 - ) - _replace_error(crsd_etree, "Rcv") + utils.crsdsar_xml_to_crsdrcv(crsd_etree) new_pvp_dtype = skcrsd.get_pvp_dtype(crsd_etree) new_pvps = np.zeros(pvps.shape, new_pvp_dtype) for field in new_pvp_dtype.fields: diff --git a/tests/core/crsd/test_scenecoords.py b/tests/core/crsd/test_scenecoords.py index 97c62f9..076da1a 100644 --- a/tests/core/crsd/test_scenecoords.py +++ b/tests/core/crsd/test_scenecoords.py @@ -1,8 +1,10 @@ +import copy import pathlib import lxml.etree import numpy as np import pytest +import shapely import sarkit.crsd as skcrsd import sarkit.wgs84 @@ -61,3 +63,81 @@ def test_derived_tofrom_iac(surf_type, xmltree_func): skcrsd.iac_to_llh(xmltree, pt_iacs[..., :2]), skcrsd.iac_to_llh(xmltree, pt_iacs * [1, 1, 0]), ) + + +@pytest.mark.parametrize("crsd_type", ("CRSDsar", "CRSDrcv")) +@pytest.mark.parametrize("use_polygon", (True, False)) +def test_image_area_funcs(crsd_type, use_polygon): + xmltree = lxml.etree.parse(DATAPATH / "example-crsd-1.0.xml") + if crsd_type == "CRSDrcv": + tests.utils.crsdsar_xml_to_crsdrcv(xmltree) + + extended_poly = shapely.Polygon([[-10, -10], [0, 20], [10, -10]]) + scene_poly = shapely.buffer(extended_poly, -1) + ch0_poly = shapely.buffer(scene_poly, -1) + ch1_poly = shapely.buffer(ch0_poly, -1) + + def set_ia(ia_ew: skcrsd.ElementWrapper, polygon: shapely.Polygon): + ia_ew["X1Y1"] = [polygon.bounds[0], polygon.bounds[1]] + ia_ew["X2Y2"] = [polygon.bounds[2], polygon.bounds[3]] + ia_ew["Polygon"] = shapely.get_coordinates(polygon)[:-1, :] + + ew = skcrsd.ElementWrapper(xmltree.getroot()) + set_ia(ew["SceneCoordinates"]["ImageArea"], scene_poly) + + def check_imgarea(actual, expected): + if use_polygon: + assert shapely.equals(shapely.Polygon(actual), expected) + else: + assert shapely.equals(shapely.Polygon(actual), expected.envelope) + + check_imgarea( + skcrsd.get_scene_image_area(xmltree, use_polygon=use_polygon), scene_poly + ) + + if crsd_type == "CRSDsar": + set_ia(ew["SceneCoordinates"]["ExtendedArea"], extended_poly) + chpar0 = ew["Channel"]["Parameters"][0] + set_ia(chpar0["SARImage"]["ImageArea"], ch0_poly) + chpar1 = copy.deepcopy(chpar0) + chpar1["Identifier"] = "chpar1_id" # assume this is unique + set_ia(chpar1["SARImage"]["ImageArea"], ch1_poly) + ew["Channel"].add("Parameters", chpar1) + chpar2 = copy.deepcopy(chpar0) + chpar2["Identifier"] = "chpar2_id" # assume this is unique + del chpar2["SARImage"]["ImageArea"] + ew["Channel"].add("Parameters", chpar2) + + check_imgarea( + skcrsd.get_channel_image_area( + xmltree, chpar0["Identifier"], use_polygon=use_polygon + ), + ch0_poly, + ) + check_imgarea( + skcrsd.get_channel_image_area( + xmltree, chpar1["Identifier"], use_polygon=use_polygon + ), + ch1_poly, + ) + check_imgarea( + skcrsd.get_channel_image_area( + xmltree, chpar2["Identifier"], use_polygon=use_polygon + ), + scene_poly, + ) + + check_imgarea( + skcrsd.get_extended_image_area(xmltree, use_polygon=use_polygon), + extended_poly, + ) + + del ew["SceneCoordinates"]["ExtendedArea"] + assert skcrsd.get_extended_image_area(xmltree, use_polygon=use_polygon) is None + else: + with pytest.raises(ValueError, match="Only CRSDsar products"): + skcrsd.get_channel_image_area( + xmltree, xmltree.findtext(".//{*}RefChId"), use_polygon=use_polygon + ) + with pytest.raises(ValueError, match="Only CRSDsar products"): + skcrsd.get_extended_image_area(xmltree, use_polygon=use_polygon) diff --git a/tests/utils.py b/tests/utils.py index 7a878d3..ceb1665 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -1,16 +1,19 @@ import asyncio import builtins import contextlib +import copy import queue import threading +import lxml.etree import numpy as np from aiohttp import web +import sarkit.crsd as skcrsd import sarkit.wgs84 -# Python's built in http.server does not support the Range header. aoihttp does +# Python's built in http.server does not support the Range header. aiohttp does def _run_aiohttp_server(app, loop, ready_event, stop_event, msg_queue): asyncio.set_event_loop(loop) runner = web.AppRunner(app) @@ -90,3 +93,110 @@ def replace_planar_with_hae(root_ew): (sarkit.wgs84.cartesian_to_geodetic(iarp_ecf + uiay) - iarp_llh)[:2] ) del sc_ew["ReferenceSurface"]["Planar"] + + +def _remove(root, pattern): + if (elem := root.find(pattern)) is not None: + elem.getparent().remove(elem) + else: + print(f"Cannot find {pattern=}") + + +def _replace_error(crsd_etree, sensor_type): + sar_error = crsd_etree.find("{*}ErrorParameters/{*}SARImage") + elem_ns = lxml.etree.QName(sar_error).namespace + retval = copy.deepcopy(sar_error.find("{*}Monostatic")) + retval.tag = f"{{{elem_ns}}}{sensor_type}Sensor" + sar_error.addnext(retval) + helper = skcrsd.XmlHelper(crsd_etree) + ndx = {"Tx": 0, "Rcv": 1}[sensor_type] + helper.set_elem( + retval.find(".//{*}TimeFreqCov"), + skcrsd.MtxType((3, 3)).parse_elem(retval.find(".//{*}TimeFreqCov"))[ + [ndx, 2], : + ][:, [ndx, 2]], + ) + time_decorr = copy.deepcopy(retval.find(f".//{{*}}{{{sensor_type}}}TimeDecorr")) + if time_decorr is not None: + time_decorr.tag = f"{{{elem_ns}}}TimeDecorr" + _remove(retval, "{*}TxTimeDecorr") + _remove(retval, "{*}RcvTimeDecorr") + retval.find(".//{*}ClockFreqDecorr").addprevious(time_decorr) + sar_error.getparent().remove(sar_error) + + +def _repack_support_arrays(crsd_etree): + offset = 0 + for array in crsd_etree.findall("{*}Data/{*}Support/{*}SupportArray"): + array.find("{*}ArrayByteOffset").text = str(offset) + offset += ( + int(array.findtext("{*}NumRows")) + * int(array.findtext("{*}NumCols")) + * int(array.findtext("{*}BytesPerElement")) + ) + return offset + + +def crsdsar_xml_to_crsdtx(crsd_etree: lxml.etree.ElementTree) -> None: + """Modify a CRSDsar ElementTree into a CRSDtx ElementTree in place.""" + crsd_etree.find(".//{*}RefPulseIndex").text = crsd_etree.find( + ".//{*}RefVectorPulseIndex" + ).text + ns = lxml.etree.QName(crsd_etree.getroot()).namespace + crsd_etree.getroot().tag = f"{{{ns}}}CRSDtx" + _remove(crsd_etree, "{*}SARInfo") + _remove(crsd_etree, "{*}ReceiveInfo") + _remove(crsd_etree, "{*}Global/{*}Receive") + _remove(crsd_etree, "{*}SceneCoordinates/{*}ExtendedArea") + _remove(crsd_etree, "{*}SceneCoordinates/{*}ImageGrid") + _remove(crsd_etree, "{*}Data/{*}Receive") + _remove(crsd_etree, "{*}Channel") + _remove(crsd_etree, "{*}ReferenceGeometry/{*}SARImage") + _remove(crsd_etree, "{*}ReferenceGeometry/{*}RcvParameters") + _remove(crsd_etree, "{*}DwellPolynomials") + _remove(crsd_etree, "{*}PVP") + _replace_error(crsd_etree, "Tx") + + +def crsdsar_xml_to_crsdrcv(crsd_etree: lxml.etree.ElementTree) -> None: + """Modify a CRSDsar ElementTree into a CRSDtx ElementTree in place.""" + ns = lxml.etree.QName(crsd_etree.getroot()).namespace + crsd_etree.getroot().tag = f"{{{ns}}}CRSDrcv" + _remove(crsd_etree, "{*}SARInfo") + _remove(crsd_etree, "{*}TransmitInfo") + _remove(crsd_etree, "{*}Global/{*}Transmit") + _remove(crsd_etree, "{*}SceneCoordinates/{*}ExtendedArea") + _remove(crsd_etree, "{*}SceneCoordinates/{*}ImageGrid") + _remove(crsd_etree, "{*}Data/{*}Transmit") + _remove(crsd_etree, "{*}TxSequence") + _remove(crsd_etree, "{*}Channel/{*}Parameters/{*}SARImage") + _remove(crsd_etree, "{*}ReferenceGeometry/{*}SARImage") + _remove(crsd_etree, "{*}ReferenceGeometry/{*}TxParameters") + _remove(crsd_etree, "{*}DwellPolynomials") + fx_ids = [ + x.text + for x in crsd_etree.findall("{*}SupportArray/{*}FxResponseArray/{*}Identifier") + ] + xm_ids = [ + x.text for x in crsd_etree.findall("{*}SupportArray/{*}XMArray/{*}Identifier") + ] + _remove(crsd_etree, "{*}SupportArray/{*}FxResponseArray") + _remove(crsd_etree, "{*}SupportArray/{*}XMArray") + for x in fx_ids + xm_ids: + _remove( + crsd_etree, + f"{{*}}Data/{{*}}Support/{{*}}SupportArray[{{*}}SAId='{x}']", + ) + nsa = crsd_etree.find("{*}Data/{*}Support/{*}NumSupportArrays") + nsa.text = str(int(nsa.text) - len(fx_ids + xm_ids)) + _repack_support_arrays(crsd_etree) + _remove(crsd_etree, "{*}PPP") + tx_pulse_index_offset = int(crsd_etree.findtext("{*}PVP/{*}TxPulseIndex/{*}Offset")) + _remove(crsd_etree, "{*}PVP/{*}TxPulseIndex") + for pvp_offset in crsd_etree.findall("{*}PVP/*/{*}Offset"): + if int(pvp_offset.text) > tx_pulse_index_offset: + pvp_offset.text = str(int(pvp_offset.text) - 1) + crsd_etree.find("{*}Data/{*}Receive/{*}NumBytesPVP").text = str( + int(crsd_etree.findtext("{*}Data/{*}Receive/{*}NumBytesPVP")) - 8 + ) + _replace_error(crsd_etree, "Rcv")