From 787314f1a62c229ed07e0e1fa80a2ef4b7091606 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Thu, 16 Jan 2025 17:25:43 -0500 Subject: [PATCH 01/13] implement roundtrip wcs<->wcsinfo --- romancal/pipeline/mosaic_pipeline.py | 102 +++++++++++++++++++-------- 1 file changed, 73 insertions(+), 29 deletions(-) diff --git a/romancal/pipeline/mosaic_pipeline.py b/romancal/pipeline/mosaic_pipeline.py index 7294ef3dd..b912df9ac 100644 --- a/romancal/pipeline/mosaic_pipeline.py +++ b/romancal/pipeline/mosaic_pipeline.py @@ -14,6 +14,8 @@ from astropy.modeling import models from gwcs import WCS, coordinate_frames +from stcal.alignment import util as wcs_util + import romancal.datamodels.filetype as filetype from romancal.datamodels import ModelLibrary @@ -114,13 +116,13 @@ def process(self, input): # check to see if there exists a skycell on disk if not create it if not isfile(skycell_file_name): - # extract the wcs info from the record for generate_tan_wcs + # extract the wcs info from the record for skycell_to_wcs log.info( "Creating skycell image at ra: %f dec %f", float(skycell_record["ra_center"]), float(skycell_record["dec_center"]), ) - skycell_wcs = generate_tan_wcs(skycell_record) + skycell_wcs = skycell_to_wcs(skycell_record) # skycell_wcs.bounding_box = bounding_box # For resample to use an external grid we need to pass it the skycell gwcs object @@ -128,6 +130,7 @@ def process(self, input): wcs_tree = {"wcs": skycell_wcs} wcs_file = asdf.AsdfFile(wcs_tree) wcs_file.write_to("skycell_wcs.asdf") + self.resample.output_wcs = "skycell_wcs.asdf" self.resample.output_shape = ( int(skycell_record["nx"]), @@ -162,41 +165,82 @@ def process(self, input): return result -def generate_tan_wcs(skycell_record): - """extract the wcs info from the record for generate_tan_wcs - we need the scale, ra, dec, bounding_box""" +def skycell_to_wcs(skycell_record): + """From a skycell record, generate a GWCS + + Parameters + ---------- + skycell_record : dict + A skycell record, or row, from the skycell patches table. + + Returns + ------- + wcsobj : wcs.GWCS + The GWCS object from the skycell record. + """ + wcsinfo = dict() + + # The scale is given in arcseconds per pixel. Convert to degrees. + wcsinfo['pixel_scale'] = float(skycell_record["pixel_scale"]) / 3600.0 - scale = float(skycell_record["pixel_scale"]) - ra_center = float(skycell_record["ra_projection_center"]) - dec_center = float(skycell_record["dec_projection_center"]) - shiftx = float(skycell_record["x0_projection"]) - shifty = float(skycell_record["y0_projection"]) + # Remaining components of the wcsinfo block + wcsinfo['ra_ref'] = float(skycell_record["ra_projection_center"]) + wcsinfo['dec_ref'] = float(skycell_record["dec_projection_center"]) + wcsinfo['x_ref'] = float(skycell_record["x0_projection"]) + wcsinfo['y_ref'] = float(skycell_record["y0_projection"]) + wcsinfo['orientat'] = float(skycell_record['orientat_projection_center']) + wcsinfo['rotation_matrix'] = None + + # Bounding box of the skycell. Note that the center of the pixels are at (0.5, 0.5) bounding_box = ( (-0.5, -0.5 + skycell_record["nx"]), (-0.5, -0.5 + skycell_record["ny"]), ) - # components of the model - # shift = models.Shift(shiftx) & models.Shift(shifty) + wcsobj = wcsinfo_to_wcs(wcsinfo, bounding_box=bounding_box) + return wcsobj - # select a scale for the skycell image, this will come from INS and may - # be optimized for the different survey programs - scale_x = scale - scale_y = scale - # set the pixelsscale to 0.1 arcsec/pixel - pixelscale = models.Scale(scale_x / 3600.0) & models.Scale(scale_y / 3600.0) - pixelshift = models.Shift(-1.0 * shiftx) & models.Shift(-1.0 * shifty) +def wcsinfo_to_wcs(wcsinfo, bounding_box=None, name='wcsinfo'): + """Create a GWCS from the L3 wcsinfo meta + + Parameters + ---------- + wcsinfo : dict or MosaicModel.meta.wcsinfo + The L3 wcsinfo to create a GWCS from. + + bounding_box : None or 4-tuple + The bounding box in detector/pixel space. Form of input is: + (x_left, x_right, y_bottom, y_top) + + name : str + Value of the `name` attribute of the GWCS object. + + Returns + ------- + wcs : wcs.GWCS + The GWCS object created. + """ + pixelshift = models.Shift(-wcsinfo['x_ref'], name='crpix1') & models.Shift(-wcsinfo['y_ref'], name='crpix2') + pixelscale = models.Scale(wcsinfo['pixel_scale'], name='cdelt1') & models.Scale(wcsinfo['pixel_scale'], name='cdelt2') tangent_projection = models.Pix2Sky_TAN() - celestial_rotation = models.RotateNative2Celestial(ra_center, dec_center, 180.0) - det2sky = pixelshift | pixelscale | tangent_projection | celestial_rotation - detector_frame = coordinate_frames.Frame2D( - name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix) - ) - sky_frame = coordinate_frames.CelestialFrame( - reference_frame=coordinates.ICRS(), name="icrs", unit=(u.deg, u.deg) - ) - wcsobj = WCS([(detector_frame, det2sky), (sky_frame, None)]) - wcsobj.bounding_box = bounding_box + celestial_rotation = models.RotateNative2Celestial(wcsinfo['ra_ref'], wcsinfo['dec_ref'], 180.) + + matrix = wcsinfo.get('rotation_matrix', None) + if matrix: + matrix = np.array(matrix) + else: + orientat = wcsinfo.get('orientat', 0.) + matrix = wcs_util.calc_rotation_matrix(np.deg2rad(orientat), v3i_yangle=0., vparity=1) + matrix = np.reshape(matrix, (2, 2)) + rotation = models.AffineTransformation2D(matrix, name='pc_rotation_matrix') + det2sky = pixelshift | rotation | pixelscale | tangent_projection | celestial_rotation + + detector_frame = coordinate_frames.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)) + sky_frame = coordinate_frames.CelestialFrame(reference_frame=coordinates.ICRS(), name='icrs', unit=(u.deg, u.deg)) + wcsobj = WCS([(detector_frame, det2sky), (sky_frame, None)], name=name) + + if bounding_box: + wcsobj.bounding_box = bounding_box return wcsobj From 179d97112dbb3aa2dff0ef8d74532c3e15271cc2 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Fri, 17 Jan 2025 08:30:10 -0500 Subject: [PATCH 02/13] Update l3 wcsinfo to ignore bounding boxes --- romancal/resample/resample.py | 13 +++++++++---- romancal/resample/tests/test_resample.py | 4 ++-- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py index 0e9ad543c..795d4cc02 100644 --- a/romancal/resample/resample.py +++ b/romancal/resample/resample.py @@ -809,13 +809,13 @@ def gwcs_into_l3(model, wcs): l3_wcsinfo.projection = "TAN" l3_wcsinfo.pixel_shape = model.shape + # Fill out image-local information pixel_center = [(v - 1) / 2.0 for v in model.shape[::-1]] world_center = wcs(*pixel_center) l3_wcsinfo.ra_center = world_center[0] l3_wcsinfo.dec_center = world_center[1] l3_wcsinfo.pixel_scale_local = compute_scale(wcs, world_center) l3_wcsinfo.orientat_local = calc_pa(wcs, *world_center) - try: footprint = utils.create_footprint(wcs, model.shape) except Exception as excp: @@ -831,6 +831,7 @@ def gwcs_into_l3(model, wcs): l3_wcsinfo.dec_corn4 = footprint[3][1] l3_wcsinfo.s_region = compute_s_region_keyword(footprint) + # Fill out wcs-general information try: l3_wcsinfo.x_ref = -transform["crpix1"].offset.value l3_wcsinfo.y_ref = -transform["crpix2"].offset.value @@ -840,7 +841,11 @@ def gwcs_into_l3(model, wcs): ) l3_wcsinfo.x_ref = pixel_center[0] l3_wcsinfo.y_ref = pixel_center[1] - world_ref = wcs(l3_wcsinfo.x_ref, l3_wcsinfo.y_ref) + + # The world reference is calculated ignoring the bounding box. Normally, + # this should not be done. However, when the wcs is determined by a skycell, + # the reference point is almost always outside the SCA. + world_ref = wcs(l3_wcsinfo.x_ref, l3_wcsinfo.y_ref, with_bounding_box=False) l3_wcsinfo.ra_ref = world_ref[0] l3_wcsinfo.dec_ref = world_ref[1] l3_wcsinfo.pixel_scale = compute_scale(wcs, world_ref) @@ -875,9 +880,9 @@ def calc_pa(wcs, ra, dec): The position angle in degrees. """ - delta_pix = [v for v in wcs.world_to_pixel_values(ra, dec)] + delta_pix = [v for v in wcs.invert(ra, dec, with_bounding_box=False)] delta_pix[1] += 1 - delta_coord = wcs.pixel_to_world(*delta_pix) + delta_coord = SkyCoord(*wcs(*delta_pix, with_bounding_box=False), frame='icrs', unit='deg') coord = SkyCoord(ra, dec, frame="icrs", unit="deg") return coord.position_angle(delta_coord).degree diff --git a/romancal/resample/tests/test_resample.py b/romancal/resample/tests/test_resample.py index 2923a2625..306e5c8ef 100644 --- a/romancal/resample/tests/test_resample.py +++ b/romancal/resample/tests/test_resample.py @@ -1084,8 +1084,8 @@ def test_l3_wcsinfo(multiple_exposures): "dec_corn3": 0.005043059486913547, "ra_corn4": 9.998944914237953, "dec_corn4": 0.00038355958555111314, - "orientat_local": 9.826983262839223, - "orientat": 9.826978421513601, + "orientat_local": 20.999999978134802, + "orientat": 20.999999978134802, "projection": "TAN", "s_region": ( "POLYGON ICRS 10.005109345783163 -0.001982743978690467 10.006897960220385 " From 7acc8a3abe9c3000bef89ab1b42c2ce7c53350ab Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Sat, 18 Jan 2025 23:51:25 -0500 Subject: [PATCH 03/13] add test for wcsinfo_to_wcs --- .../pipeline/tests/test_mosaic_pipeline.py | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 romancal/pipeline/tests/test_mosaic_pipeline.py diff --git a/romancal/pipeline/tests/test_mosaic_pipeline.py b/romancal/pipeline/tests/test_mosaic_pipeline.py new file mode 100644 index 000000000..1664c0aeb --- /dev/null +++ b/romancal/pipeline/tests/test_mosaic_pipeline.py @@ -0,0 +1,36 @@ +"""Unit tests for the mosaic pipeline""" +import numpy as np + +import romancal.pipeline.mosaic_pipeline as mp + +def test_wcsinfo_to_wcs(): + """Test integrity of wcsinfo_to_wcs""" + wcsinfo = { + 'ra_ref': 269.83219987378925, + 'dec_ref': 66.04081466149024, + 'x_ref': 2069.0914958388985, + 'y_ref': 2194.658767532754, + 'rotation_matrix': [[-0.9999964196507396, -0.00267594575838714], [-0.00267594575838714, 0.9999964196507396]], + 'pixel_scale': 3.036307317109957e-05, + 'pixel_shape': [4389, 4138], + 'ra_center': 269.82284964811464, + 'dec_center': 66.0369888162117, + 'ra_corn1': 269.98694025887136, + 'dec_corn1': 65.97426875366378, + 'ra_corn2': 269.98687579251805, + 'dec_corn2': 66.09988065827382, + 'ra_corn3': 269.6596332616879, + 'dec_corn3': 65.97389321243348, + 'ra_corn4': 269.6579498847431, + 'dec_corn4': 66.099533603104, + 'orientat': 359.8466793994546 + } + + wcs = mp.wcsinfo_to_wcs(wcsinfo) + + assert np.allclose(wcs(wcsinfo['x_ref'], wcsinfo['y_ref']), (wcsinfo['ra_ref'], wcsinfo['dec_ref'])) + assert np.allclose(wcs(4389 / 2., 4138 / 2.), (wcsinfo['ra_center'], wcsinfo['dec_center'])) + assert np.allclose(wcs(0., 0.), (wcsinfo['ra_corn1'], wcsinfo['dec_corn1'])) + assert np.allclose(wcs(0., wcsinfo['pixel_shape'][1] - 1.), (wcsinfo['ra_corn2'], wcsinfo['dec_corn2'])) + assert np.allclose(wcs(wcsinfo['pixel_shape'][0], 0.), (wcsinfo['ra_corn3'], wcsinfo['dec_corn3'])) + assert np.allclose(wcs(wcsinfo['pixel_shape'][0], wcsinfo['pixel_shape'][1]), (wcsinfo['ra_corn4'], wcsinfo['dec_corn4'])) From 01a30b52a25f2474566e4ca03a88c467ea4b5260 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Sun, 19 Jan 2025 00:23:33 -0500 Subject: [PATCH 04/13] add test for skycell_to_wcs and fixed bug in other test --- .../pipeline/tests/test_mosaic_pipeline.py | 45 ++++++++++++++++--- 1 file changed, 38 insertions(+), 7 deletions(-) diff --git a/romancal/pipeline/tests/test_mosaic_pipeline.py b/romancal/pipeline/tests/test_mosaic_pipeline.py index 1664c0aeb..3d1cc4788 100644 --- a/romancal/pipeline/tests/test_mosaic_pipeline.py +++ b/romancal/pipeline/tests/test_mosaic_pipeline.py @@ -3,6 +3,37 @@ import romancal.pipeline.mosaic_pipeline as mp + +def test_skycell_to_wcs(): + """Test integrity of skycell_to_wcs""" + + skycell = np.void( + ('r274dp63x31y81', 269.7783307416819, 66.04965143695566, + 1781.5, 1781.5, 355.9788, 3564, 3564, 67715.5, -110484.5, + 269.6657957545588, 65.9968687812357, 269.6483032937494, + 66.09523979539262, 269.89132874168854, 66.10234971630734, + 269.9079118635897, 66.00394719483091, + 0.1, 274.2857142857143, 63.0, 0.0, 463181), + dtype=[('name', ' Date: Sun, 19 Jan 2025 06:17:39 +0000 Subject: [PATCH 05/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/pipeline/mosaic_pipeline.py | 51 +++--- .../pipeline/tests/test_mosaic_pipeline.py | 154 +++++++++++++----- romancal/resample/resample.py | 4 +- 3 files changed, 145 insertions(+), 64 deletions(-) diff --git a/romancal/pipeline/mosaic_pipeline.py b/romancal/pipeline/mosaic_pipeline.py index b912df9ac..fb5f2315f 100644 --- a/romancal/pipeline/mosaic_pipeline.py +++ b/romancal/pipeline/mosaic_pipeline.py @@ -13,7 +13,6 @@ from astropy import units as u from astropy.modeling import models from gwcs import WCS, coordinate_frames - from stcal.alignment import util as wcs_util import romancal.datamodels.filetype as filetype @@ -181,15 +180,15 @@ def skycell_to_wcs(skycell_record): wcsinfo = dict() # The scale is given in arcseconds per pixel. Convert to degrees. - wcsinfo['pixel_scale'] = float(skycell_record["pixel_scale"]) / 3600.0 + wcsinfo["pixel_scale"] = float(skycell_record["pixel_scale"]) / 3600.0 # Remaining components of the wcsinfo block - wcsinfo['ra_ref'] = float(skycell_record["ra_projection_center"]) - wcsinfo['dec_ref'] = float(skycell_record["dec_projection_center"]) - wcsinfo['x_ref'] = float(skycell_record["x0_projection"]) - wcsinfo['y_ref'] = float(skycell_record["y0_projection"]) - wcsinfo['orientat'] = float(skycell_record['orientat_projection_center']) - wcsinfo['rotation_matrix'] = None + wcsinfo["ra_ref"] = float(skycell_record["ra_projection_center"]) + wcsinfo["dec_ref"] = float(skycell_record["dec_projection_center"]) + wcsinfo["x_ref"] = float(skycell_record["x0_projection"]) + wcsinfo["y_ref"] = float(skycell_record["y0_projection"]) + wcsinfo["orientat"] = float(skycell_record["orientat_projection_center"]) + wcsinfo["rotation_matrix"] = None # Bounding box of the skycell. Note that the center of the pixels are at (0.5, 0.5) bounding_box = ( @@ -201,7 +200,7 @@ def skycell_to_wcs(skycell_record): return wcsobj -def wcsinfo_to_wcs(wcsinfo, bounding_box=None, name='wcsinfo'): +def wcsinfo_to_wcs(wcsinfo, bounding_box=None, name="wcsinfo"): """Create a GWCS from the L3 wcsinfo meta Parameters @@ -221,23 +220,37 @@ def wcsinfo_to_wcs(wcsinfo, bounding_box=None, name='wcsinfo'): wcs : wcs.GWCS The GWCS object created. """ - pixelshift = models.Shift(-wcsinfo['x_ref'], name='crpix1') & models.Shift(-wcsinfo['y_ref'], name='crpix2') - pixelscale = models.Scale(wcsinfo['pixel_scale'], name='cdelt1') & models.Scale(wcsinfo['pixel_scale'], name='cdelt2') + pixelshift = models.Shift(-wcsinfo["x_ref"], name="crpix1") & models.Shift( + -wcsinfo["y_ref"], name="crpix2" + ) + pixelscale = models.Scale(wcsinfo["pixel_scale"], name="cdelt1") & models.Scale( + wcsinfo["pixel_scale"], name="cdelt2" + ) tangent_projection = models.Pix2Sky_TAN() - celestial_rotation = models.RotateNative2Celestial(wcsinfo['ra_ref'], wcsinfo['dec_ref'], 180.) + celestial_rotation = models.RotateNative2Celestial( + wcsinfo["ra_ref"], wcsinfo["dec_ref"], 180.0 + ) - matrix = wcsinfo.get('rotation_matrix', None) + matrix = wcsinfo.get("rotation_matrix", None) if matrix: matrix = np.array(matrix) else: - orientat = wcsinfo.get('orientat', 0.) - matrix = wcs_util.calc_rotation_matrix(np.deg2rad(orientat), v3i_yangle=0., vparity=1) + orientat = wcsinfo.get("orientat", 0.0) + matrix = wcs_util.calc_rotation_matrix( + np.deg2rad(orientat), v3i_yangle=0.0, vparity=1 + ) matrix = np.reshape(matrix, (2, 2)) - rotation = models.AffineTransformation2D(matrix, name='pc_rotation_matrix') - det2sky = pixelshift | rotation | pixelscale | tangent_projection | celestial_rotation + rotation = models.AffineTransformation2D(matrix, name="pc_rotation_matrix") + det2sky = ( + pixelshift | rotation | pixelscale | tangent_projection | celestial_rotation + ) - detector_frame = coordinate_frames.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)) - sky_frame = coordinate_frames.CelestialFrame(reference_frame=coordinates.ICRS(), name='icrs', unit=(u.deg, u.deg)) + detector_frame = coordinate_frames.Frame2D( + name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix) + ) + sky_frame = coordinate_frames.CelestialFrame( + reference_frame=coordinates.ICRS(), name="icrs", unit=(u.deg, u.deg) + ) wcsobj = WCS([(detector_frame, det2sky), (sky_frame, None)], name=name) if bounding_box: diff --git a/romancal/pipeline/tests/test_mosaic_pipeline.py b/romancal/pipeline/tests/test_mosaic_pipeline.py index 3d1cc4788..79bcd2d9b 100644 --- a/romancal/pipeline/tests/test_mosaic_pipeline.py +++ b/romancal/pipeline/tests/test_mosaic_pipeline.py @@ -1,4 +1,5 @@ """Unit tests for the mosaic pipeline""" + import numpy as np import romancal.pipeline.mosaic_pipeline as mp @@ -8,60 +9,125 @@ def test_skycell_to_wcs(): """Test integrity of skycell_to_wcs""" skycell = np.void( - ('r274dp63x31y81', 269.7783307416819, 66.04965143695566, - 1781.5, 1781.5, 355.9788, 3564, 3564, 67715.5, -110484.5, - 269.6657957545588, 65.9968687812357, 269.6483032937494, - 66.09523979539262, 269.89132874168854, 66.10234971630734, - 269.9079118635897, 66.00394719483091, - 0.1, 274.2857142857143, 63.0, 0.0, 463181), - dtype=[('name', ' Date: Sun, 19 Jan 2025 01:21:01 -0500 Subject: [PATCH 06/13] update changelog --- changes/1585.mosaic_pipeline.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changes/1585.mosaic_pipeline.rst diff --git a/changes/1585.mosaic_pipeline.rst b/changes/1585.mosaic_pipeline.rst new file mode 100644 index 000000000..977888ff6 --- /dev/null +++ b/changes/1585.mosaic_pipeline.rst @@ -0,0 +1 @@ +Roundtrip L3 wcsinfo especially when skycell specifications are used From f3825fcbe48418777c94054d45cb734a2640c8f2 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Tue, 21 Jan 2025 23:02:28 -0500 Subject: [PATCH 07/13] add round trip tests --- romancal/regtest/test_mos_pipeline.py | 11 +++++++ romancal/regtest/test_mos_skycell_pipeline.py | 11 +++++++ romancal/regtest/util.py | 32 +++++++++++++++++++ 3 files changed, 54 insertions(+) create mode 100644 romancal/regtest/util.py diff --git a/romancal/regtest/test_mos_pipeline.py b/romancal/regtest/test_mos_pipeline.py index 8fba060af..c4dfbd3df 100644 --- a/romancal/regtest/test_mos_pipeline.py +++ b/romancal/regtest/test_mos_pipeline.py @@ -5,8 +5,10 @@ import pytest import roman_datamodels as rdm +from romancal.pipeline import mosaic_pipeline from romancal.pipeline.mosaic_pipeline import MosaicPipeline +from . import util from .regtestdata import compare_asdf # mark all tests in this module @@ -116,3 +118,12 @@ def test_added_background(output_model): def test_added_background_level(output_model): # DMS400 assert any(output_model.meta.individual_image_meta.background["level"] != 0) + + +def test_wcsinfo_wcs_roundtrip(output_model): + """Test that the contents of wcsinfo reproduces the wcs""" + wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo) + + ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo) + assert ra_mad < 5.0e-4 + assert dec_mad < 5.0e-4 diff --git a/romancal/regtest/test_mos_skycell_pipeline.py b/romancal/regtest/test_mos_skycell_pipeline.py index ce5953d96..0f4a1d918 100644 --- a/romancal/regtest/test_mos_skycell_pipeline.py +++ b/romancal/regtest/test_mos_skycell_pipeline.py @@ -1,8 +1,10 @@ import pytest import roman_datamodels as rdm +from romancal.pipeline import mosaic_pipeline from romancal.pipeline.mosaic_pipeline import MosaicPipeline +from . import util from .regtestdata import compare_asdf # mark all tests in this module @@ -56,3 +58,12 @@ def test_resample_ran(output_model): def test_location_name(output_model): # test that the location_name matches the skycell selected assert output_model.meta.basic.location_name == "r274dp63x31y81" + + +def test_wcsinfo_wcs_roundtrip(output_model): + """Test that the contents of wcsinfo reproduces the wcs""" + wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo) + + ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo) + assert ra_mad < 5.0e-4 + assert dec_mad < 5.0e-4 diff --git a/romancal/regtest/util.py b/romancal/regtest/util.py new file mode 100644 index 000000000..27ef1997a --- /dev/null +++ b/romancal/regtest/util.py @@ -0,0 +1,32 @@ +"""Test utilities""" +from astropy.stats import mad_std +import numpy as np + + +def comp_wcs_grids_arcs(wcs_a, wcs_b, npix=4088, interval=10): + """Compare world grids produced by the two wcs + + Parameters + ---------- + wcs_a, wcs_b : gwcs.WCS + The wcs object to compare. + + npix : int + The size of the grid to produce. + + interval : int + The interval to check over. + + Returns + ------- + mad_std : float + The numpy MAD_STD in arcseconds + """ + xx, yy = np.meshgrid(np.linspace(0, npix, interval), np.linspace(0, npix, interval)) + ra_a, dec_a = wcs_a(xx, yy, with_bounding_box=False) + ra_b, dec_b = wcs_b(xx, yy, with_bounding_box=False) + + ra_mad = mad_std(ra_a - ra_b, ignore_nan=True) * 60. * 60. * 1000. + dec_mad = mad_std(dec_a - dec_b, ignore_nan=True) * 60. * 60. * 1000. + + return ra_mad, dec_mad From 7bec6a635de919582102f0b4b1f183aba2f98548 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 22 Jan 2025 04:02:58 +0000 Subject: [PATCH 08/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/regtest/util.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/romancal/regtest/util.py b/romancal/regtest/util.py index 27ef1997a..31bf2605f 100644 --- a/romancal/regtest/util.py +++ b/romancal/regtest/util.py @@ -1,6 +1,7 @@ """Test utilities""" -from astropy.stats import mad_std + import numpy as np +from astropy.stats import mad_std def comp_wcs_grids_arcs(wcs_a, wcs_b, npix=4088, interval=10): @@ -26,7 +27,7 @@ def comp_wcs_grids_arcs(wcs_a, wcs_b, npix=4088, interval=10): ra_a, dec_a = wcs_a(xx, yy, with_bounding_box=False) ra_b, dec_b = wcs_b(xx, yy, with_bounding_box=False) - ra_mad = mad_std(ra_a - ra_b, ignore_nan=True) * 60. * 60. * 1000. - dec_mad = mad_std(dec_a - dec_b, ignore_nan=True) * 60. * 60. * 1000. + ra_mad = mad_std(ra_a - ra_b, ignore_nan=True) * 60.0 * 60.0 * 1000.0 + dec_mad = mad_std(dec_a - dec_b, ignore_nan=True) * 60.0 * 60.0 * 1000.0 return ra_mad, dec_mad From a5c4fdc46a9ab5f9dd8d209ad73e5c8ca2cbf7a3 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Wed, 22 Jan 2025 14:11:14 -0500 Subject: [PATCH 09/13] update tests to combine stats into a single measure --- romancal/regtest/test_mos_pipeline.py | 3 +-- romancal/regtest/test_mos_skycell_pipeline.py | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/romancal/regtest/test_mos_pipeline.py b/romancal/regtest/test_mos_pipeline.py index c4dfbd3df..f797efd68 100644 --- a/romancal/regtest/test_mos_pipeline.py +++ b/romancal/regtest/test_mos_pipeline.py @@ -125,5 +125,4 @@ def test_wcsinfo_wcs_roundtrip(output_model): wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo) ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo) - assert ra_mad < 5.0e-4 - assert dec_mad < 5.0e-4 + assert (ra_mad + dec_mad) / 2. < 1.0e-4 diff --git a/romancal/regtest/test_mos_skycell_pipeline.py b/romancal/regtest/test_mos_skycell_pipeline.py index 0f4a1d918..0df383371 100644 --- a/romancal/regtest/test_mos_skycell_pipeline.py +++ b/romancal/regtest/test_mos_skycell_pipeline.py @@ -65,5 +65,4 @@ def test_wcsinfo_wcs_roundtrip(output_model): wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo) ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo) - assert ra_mad < 5.0e-4 - assert dec_mad < 5.0e-4 + assert (ra_mad + dec_mad) / 2.0 < 1.0e-4 From 348f36d65fb85b793c58e6d21903a130b1fccb29 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 22 Jan 2025 19:11:41 +0000 Subject: [PATCH 10/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/regtest/test_mos_pipeline.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/romancal/regtest/test_mos_pipeline.py b/romancal/regtest/test_mos_pipeline.py index f797efd68..35fde6077 100644 --- a/romancal/regtest/test_mos_pipeline.py +++ b/romancal/regtest/test_mos_pipeline.py @@ -125,4 +125,4 @@ def test_wcsinfo_wcs_roundtrip(output_model): wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo) ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo) - assert (ra_mad + dec_mad) / 2. < 1.0e-4 + assert (ra_mad + dec_mad) / 2.0 < 1.0e-4 From 8c93508c08200d545dd8e0a855176dfa43b1c2d1 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Mon, 27 Jan 2025 08:59:42 -0500 Subject: [PATCH 11/13] Get scale explicitly --- romancal/resample/resample.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py index 3af7f9220..83a1a8791 100644 --- a/romancal/resample/resample.py +++ b/romancal/resample/resample.py @@ -842,13 +842,17 @@ def gwcs_into_l3(model, wcs): l3_wcsinfo.x_ref = pixel_center[0] l3_wcsinfo.y_ref = pixel_center[1] - # The world reference is calculated ignoring the bounding box. Normally, - # this should not be done. However, when the wcs is determined by a skycell, - # the reference point is almost always outside the SCA. world_ref = wcs(l3_wcsinfo.x_ref, l3_wcsinfo.y_ref, with_bounding_box=False) l3_wcsinfo.ra_ref = world_ref[0] l3_wcsinfo.dec_ref = world_ref[1] - l3_wcsinfo.pixel_scale = compute_scale(wcs, world_ref) + + try: + cdelt1 = transform['cdelt1'].factor.value + cdelt2 = transform['cdelt2'].factor.value + l3_wcsinfo.pixel_scale = (cdelt1 + cdelt2) / 2. + except IndexError: + l3_wcsinfo.pixel_scale = compute_scale(wcs, world_ref) + l3_wcsinfo.orientat = calc_pa(wcs, *world_ref) try: From fbfa7df0238bf21e906ce900e1465967d1812497 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Mon, 27 Jan 2025 09:23:24 -0500 Subject: [PATCH 12/13] make test more stringent --- romancal/regtest/test_mos_pipeline.py | 2 +- romancal/regtest/test_mos_skycell_pipeline.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/romancal/regtest/test_mos_pipeline.py b/romancal/regtest/test_mos_pipeline.py index 35fde6077..dfa2247b7 100644 --- a/romancal/regtest/test_mos_pipeline.py +++ b/romancal/regtest/test_mos_pipeline.py @@ -125,4 +125,4 @@ def test_wcsinfo_wcs_roundtrip(output_model): wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo) ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo) - assert (ra_mad + dec_mad) / 2.0 < 1.0e-4 + assert (ra_mad + dec_mad) / 2.0 < 1.0e-5 diff --git a/romancal/regtest/test_mos_skycell_pipeline.py b/romancal/regtest/test_mos_skycell_pipeline.py index 0df383371..eddbc3a31 100644 --- a/romancal/regtest/test_mos_skycell_pipeline.py +++ b/romancal/regtest/test_mos_skycell_pipeline.py @@ -65,4 +65,4 @@ def test_wcsinfo_wcs_roundtrip(output_model): wcs_from_wcsinfo = mosaic_pipeline.wcsinfo_to_wcs(output_model.meta.wcsinfo) ra_mad, dec_mad = util.comp_wcs_grids_arcs(output_model.meta.wcs, wcs_from_wcsinfo) - assert (ra_mad + dec_mad) / 2.0 < 1.0e-4 + assert (ra_mad + dec_mad) / 2.0 < 1.0e-5 From 78770f6c3dcd9f4ce9aa3ce05c08af4329907575 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 27 Jan 2025 14:23:58 +0000 Subject: [PATCH 13/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/resample/resample.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py index 83a1a8791..48749ebea 100644 --- a/romancal/resample/resample.py +++ b/romancal/resample/resample.py @@ -847,9 +847,9 @@ def gwcs_into_l3(model, wcs): l3_wcsinfo.dec_ref = world_ref[1] try: - cdelt1 = transform['cdelt1'].factor.value - cdelt2 = transform['cdelt2'].factor.value - l3_wcsinfo.pixel_scale = (cdelt1 + cdelt2) / 2. + cdelt1 = transform["cdelt1"].factor.value + cdelt2 = transform["cdelt2"].factor.value + l3_wcsinfo.pixel_scale = (cdelt1 + cdelt2) / 2.0 except IndexError: l3_wcsinfo.pixel_scale = compute_scale(wcs, world_ref)