diff --git a/drizzle/tests/test_utils.py b/drizzle/tests/test_utils.py index a8dde22..68dac2b 100644 --- a/drizzle/tests/test_utils.py +++ b/drizzle/tests/test_utils.py @@ -111,65 +111,11 @@ def test_estimate_pixel_scale_ratio(): input_file = os.path.join(DATA_DIR, 'j8bt06nyq_flt.fits') with fits.open(input_file) as h: - hdr = h[1].header - assert abs(hdr['CRPIX1'] / 512 - 1) < 1e-14 - assert abs(hdr['CRPIX2'] / 512 - 1) < 1e-14 - assert abs(hdr['CRVAL1'] / 6.027148333333000000000 - 1) < 1e-14 - assert abs(hdr['CRVAL2'] / 72.08351111111000000000 + 1) < 1e-14 - assert abs(hdr['CD1_1'] / 7.360500000000000000000E-06 + 1) < 1e-14 - assert abs(hdr['CD1_2'] / 1.845600000000000000000E-06 - 1) < 1e-14 - assert abs(hdr['CD2_1'] / 2.868580000000000000000E-06 - 1) < 1e-14 - assert abs(hdr['CD2_2'] / 6.649460000000000000000E-06 - 1) < 1e-14 w = wcs.WCS(h[1].header) - assert np.allclose(w.wcs.crval, [6.027148333333000000000, -72.08351111111000000000], rtol=1e-14) - assert np.allclose(w.wcs.crpix, [512.0, 512.0], rtol=1e-14) - assert np.allclose( - w.wcs.cd, - [ - [-7.360500000000000000000E-06, 1.845600000000000000000E-06], - [2.868580000000000000000E-06, 6.649460000000000000000E-06] - ], - rtol=1e-14 - ) - - refpix = np.array([0., 0.]) - l1, phi1 = w.pixel_to_world_values(*(refpix - 0.5)) - l2, phi2 = w.pixel_to_world_values(*(refpix + [-0.5, 0.5])) - l3, phi3 = w.pixel_to_world_values(*(refpix + 0.5)) - l4, phi4 = w.pixel_to_world_values(*(refpix + [0.5, -0.5])) - assert abs(l1 - 6.0363204188670885) < 1e-14 - assert abs(l2 - 6.036326416554582) < 1e-14 - assert abs(l3 - 6.036302482421591) < 1e-14 - assert abs(l4 - 6.036296484726444) < 1e-14 - assert abs(phi1 - -72.08837937371479) < 1e-14 - assert abs(phi2 - -72.08837272397372) < 1e-14 - assert abs(phi3 - -72.08836985651423) < 1e-14 - assert abs(phi4 - -72.08837650625458) < 1e-14 - #pscale1 = _estimate_pixel_scale(w, refpix) - #assert abs(pscale1 - 1.285368350331663e-07) < 1.e-21 - - refpix = np.array([512., 512.]) - l1, phi1 = w.pixel_to_world_values(*(refpix - 0.5)) - l2, phi2 = w.pixel_to_world_values(*(refpix + [-0.5, 0.5])) - l3, phi3 = w.pixel_to_world_values(*(refpix + 0.5)) - l4, phi4 = w.pixel_to_world_values(*(refpix + [0.5, -0.5])) - assert abs(l1 - 6.027139369821101) < 1e-14 - assert abs(l2 - 6.027145369226525) < 1e-14 - assert abs(l3 - 6.027121442811089) < 1e-14 - assert abs(l4 - 6.027115443398036) < 1e-14 - assert abs(phi1 - -72.08350635208977) < 1e-13 - assert abs(phi2 - -72.08349970262996) < 1e-13 - assert abs(phi3 - -72.08349683404813) < 1e-13 - assert abs(phi4 - -72.08350348350724) < 1e-13 - #pscale2 = _estimate_pixel_scale(w, refpix) - #assert abs(pscale2 - 1.285368361004753e-07) < 1.e-21 pscale = estimate_pixel_scale_ratio(w, w, w.wcs.crpix, (0, 0)) - #assert (pscale1 / pscale2 - 0.9999999916964737) < 1e-14 - - # assert abs(pscale - 0.9999999916967218) < 1e-14 # if using numpy in estimate_pixel_scale_ratio - assert abs(pscale - 0.9999999916964737) < 1e-14 + assert abs(pscale - 0.9999999916964737) < 1.0e-9 def test_estimate_pixel_scale_no_refpix(): diff --git a/drizzle/utils.py b/drizzle/utils.py index 6ad6656..a3362af 100644 --- a/drizzle/utils.py +++ b/drizzle/utils.py @@ -127,13 +127,12 @@ def estimate_pixel_scale_ratio(wcs_from, wcs_to, refpix_from=None, refpix_to=Non """ pscale_ratio = (_estimate_pixel_scale(wcs_to, refpix_to) / _estimate_pixel_scale(wcs_from, refpix_from)) - return pscale_ratio.astype("float") + return pscale_ratio def _estimate_pixel_scale(wcs, refpix): # estimate pixel scale (in rad) using approximate algorithm # from https://trs.jpl.nasa.gov/handle/2014/40409 - if refpix is None: if hasattr(wcs, 'bounding_box') and wcs.bounding_box is not None: refpix = np.mean(wcs.bounding_box, axis=-1) @@ -146,17 +145,17 @@ def _estimate_pixel_scale(wcs, refpix): else: refpix = np.asarray(refpix) - l1, phi1 = np.array(wcs.pixel_to_world_values(*(refpix - 0.5)), dtype='longdouble') - l2, phi2 = np.array(wcs.pixel_to_world_values(*(refpix + [-0.5, 0.5])), dtype='longdouble') - l3, phi3 = np.array(wcs.pixel_to_world_values(*(refpix + 0.5)), dtype='longdouble') - l4, phi4 = np.array(wcs.pixel_to_world_values(*(refpix + [0.5, -0.5])), dtype='longdouble') + l1, phi1 = wcs.pixel_to_world_values(*(refpix - 0.5)) + l2, phi2 = wcs.pixel_to_world_values(*(refpix + [-0.5, 0.5])) + l3, phi3 = wcs.pixel_to_world_values(*(refpix + 0.5)) + l4, phi4 = wcs.pixel_to_world_values(*(refpix + [0.5, -0.5])) area = _DEG2RAD * abs( 0.5 * ( - (l4 - l2) * (np.sin(_DEG2RAD * phi1) - np.sin(_DEG2RAD * phi3)) + - (l1 - l3) * (np.sin(_DEG2RAD * phi2) - np.sin(_DEG2RAD * phi4)) + (l4 - l2) * (math.sin(_DEG2RAD * phi1) - math.sin(_DEG2RAD * phi3)) + + (l1 - l3) * (math.sin(_DEG2RAD * phi2) - math.sin(_DEG2RAD * phi4)) ) ) - return np.sqrt(area) + return math.sqrt(area) def decode_context(context, x, y):