From 7ade8dd94a8982f64ff6dc8e903ef10dae4f8421 Mon Sep 17 00:00:00 2001 From: Daniel Pressler Date: Tue, 22 Sep 2026 15:35:32 -0700 Subject: [PATCH] Add sicd_help.recompute_deltak --- CHANGELOG.md | 3 ++ sarkit_processing/sicd_help.py | 82 ++++++++++++++++++++++++++++++++++ tests/test_sicd_help.py | 80 +++++++++++++++++++++++++++++++++ 3 files changed, 165 insertions(+) create mode 100644 sarkit_processing/sicd_help.py create mode 100644 tests/test_sicd_help.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 194ea59..a78260d 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 +- `sicd_help` module + ## [0.2.0] - 2026-08-27 diff --git a/sarkit_processing/sicd_help.py b/sarkit_processing/sicd_help.py new file mode 100644 index 0000000..3562b6d --- /dev/null +++ b/sarkit_processing/sicd_help.py @@ -0,0 +1,82 @@ +"""Module for miscellaneous SICD functions.""" + +import copy +import functools + +import lxml.etree +import numpy as np +import numpy.polynomial.polynomial as npp +import sarkit.sicd as sksicd +import shapely + + +def _get_samples_in_poly(poly: shapely.Polygon, grid_size: int = 11) -> np.ndarray: + """Return samples that intersect a polygon.""" + bounds = np.asarray(poly.bounds).reshape(2, 2) # [[xmin, ymin], [xmax, ymax]] + mesh = np.stack( + np.meshgrid( + np.linspace(bounds[0, 0], bounds[1, 0], grid_size), + np.linspace(bounds[0, 1], bounds[1, 1], grid_size), + ), + axis=-1, + ) + inner_mesh = shapely.get_coordinates(poly.intersection(shapely.multipoints(mesh))) + poly_vertices = shapely.get_coordinates(poly.exterior)[:-1] + return np.concatenate( + [inner_mesh, poly_vertices], + axis=0, + ) + + +def recompute_deltak(sicd_xmltree: lxml.etree.ElementTree) -> lxml.etree.ElementTree: + """Return a copy of a SICD XML with recomputed Grid//DeltaK1 and Grid//DeltaK2.""" + sicd_xmltree = copy.deepcopy(sicd_xmltree) + ew = sksicd.ElementWrapper(sicd_xmltree.getroot()) + + if sicd_xmltree.find("{*}Grid/*/{*}DeltaKCOAPoly") is not None: + nrows = ew["ImageData"]["NumRows"] + ncols = ew["ImageData"]["NumCols"] + r0 = ew["ImageData"]["FirstRow"] + c0 = ew["ImageData"]["FirstCol"] + nrows_fi = ew["ImageData"]["FullImage"]["NumRows"] + ncols_fi = ew["ImageData"]["FullImage"]["NumCols"] + + polygons_rowcol = [ + shapely.box(-0.5, -0.5, nrows_fi - 0.5, ncols_fi - 0.5), + shapely.box(r0 - 0.5, c0 - 0.5, r0 + nrows - 0.5, c0 + ncols - 0.5), + ] + validdata = ew["ImageData"].get("ValidData", None) + if validdata is not None: + polygons_rowcol.append(shapely.Polygon(validdata).buffer(0.5)) + polygons_xrowycol = [ + shapely.transform( + p, functools.partial(sksicd.rowcol_to_xrowycol, sicd_xmltree) + ) + for p in polygons_rowcol + ] + polygon_to_sample = shapely.intersection_all(polygons_xrowycol) + pts_xrowycol = _get_samples_in_poly(polygon_to_sample) + + for rowcol in ("Row", "Col"): + griddir = ew["Grid"][rowcol] + dkcoapoly = griddir.get("DeltaKCOAPoly", None) + if dkcoapoly is not None: + dkcoa = npp.polyval2d(pts_xrowycol[:, 0], pts_xrowycol[:, 1], dkcoapoly) + dkcoa_min = dkcoa.min() + dkcoa_max = dkcoa.max() + else: + dkcoa_min = 0.0 + dkcoa_max = 0.0 + dk1 = dkcoa_min - (griddir["ImpRespBW"] / 2.0) + dk2 = dkcoa_max + (griddir["ImpRespBW"] / 2.0) + + # saturate aliased spectrum + dk_nyq = 0.5 / griddir["SS"] + if (dk1 < -dk_nyq) or (dk2 > dk_nyq): + dk1 = -dk_nyq + dk2 = dk_nyq + + griddir["DeltaK1"] = dk1 + griddir["DeltaK2"] = dk2 + + return sicd_xmltree diff --git a/tests/test_sicd_help.py b/tests/test_sicd_help.py new file mode 100644 index 0000000..8d75fb2 --- /dev/null +++ b/tests/test_sicd_help.py @@ -0,0 +1,80 @@ +import pathlib + +import lxml.etree +import numpy as np +import sarkit.sicd as sksicd + +import sarkit_processing.sicd_help as skp_sicdhelp + +sicd_xml_path = ( + pathlib.Path(__file__).absolute().parents[1] / "data/example-sicd-1.3.0.xml" +) + + +def test_recompute_deltak(): + # Contrive a test case to hit the various optional metadata paths + sicdxml = lxml.etree.parse(sicd_xml_path) + schema = lxml.etree.XMLSchema( + file=sksicd.VERSION_INFO[lxml.etree.QName(sicdxml.getroot()).namespace][ + "schema" + ] + ) + schema.assertValid(sicdxml) + ew = sksicd.ElementWrapper(sicdxml.getroot()) + # ValidData < Image < FullImage + ew["ImageData"].from_dict( + { + "NumRows": 5, + "NumCols": 7, + "FirstRow": 2, + "FirstCol": 2, + "FullImage": {"NumRows": 8, "NumCols": 10}, + "SCPPixel": [0, 0], + "ValidData": [[3, 3], [3, 7], [5, 7], [5, 3]], + } + ) + ew["Grid"]["Row"].from_dict( + {"SS": 1.0, "ImpRespBW": 0.1, "DeltaKCOAPoly": [[0.0, 0.0], [1e-2, 0.0]]} + ) + ew["Grid"]["Col"].from_dict({"SS": 2.0, "ImpRespBW": 0.2}) + del ew["Grid"]["Col"]["DeltaKCOAPoly"] + + def recompute(): + newxml = skp_sicdhelp.recompute_deltak(sicdxml) + schema.assertValid(newxml) + gridew = sksicd.ElementWrapper(newxml.find("{*}Grid")) + return [(gridew[rc]["DeltaK1"], gridew[rc]["DeltaK2"]) for rc in ("Row", "Col")] + + def within(a, b): + return (a[0] > b[0]) and (a[1] < b[1]) + + rowspan_nyq = np.array([-0.5, 0.5]) / ew["Grid"]["Row"]["SS"] + rowspan_nyq = np.array([-0.5, 0.5]) / ew["Grid"]["Col"]["SS"] + rowspan_orig, colspan_orig = recompute() + assert within(rowspan_orig, rowspan_nyq) + assert within(colspan_orig, rowspan_nyq) + + # remove ValidData + del ew["ImageData"]["ValidData"] + rowspan_novalid, colspan_novalid = recompute() + assert within(rowspan_orig, rowspan_novalid) + assert np.array_equal(colspan_novalid, colspan_orig) + + # make Image = FullImage + ew["ImageData"].from_dict( + { + "NumRows": ew["ImageData"]["FullImage"]["NumRows"], + "NumCols": ew["ImageData"]["FullImage"]["NumCols"], + "FirstRow": 0, + "FirstCol": 0, + } + ) + rowspan_fullimg, colspan_fullimg = recompute() + assert within(rowspan_novalid, rowspan_fullimg) + assert np.array_equal(colspan_fullimg, colspan_orig) + + # alias case + ew["Grid"]["Col"]["ImpRespBW"] = 1.0 + rowspan_aliascol, colspan_aliascol = recompute() + assert rowspan_aliascol == rowspan_fullimg + assert np.array_equal(colspan_aliascol, rowspan_nyq)