From 4c2b6422014e8e65b222a61dcae877c11954a1cc Mon Sep 17 00:00:00 2001 From: "Thomson, Alec (CASS, Kensington)" Date: Sat, 11 May 2024 19:22:31 +1000 Subject: [PATCH 1/6] Use skycoords --- arrakis/cutout.py | 89 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 69 insertions(+), 20 deletions(-) diff --git a/arrakis/cutout.py b/arrakis/cutout.py index 1d1acf14..f6fe5479 100644 --- a/arrakis/cutout.py +++ b/arrakis/cutout.py @@ -148,19 +148,25 @@ def cutout_image( padder: float = cube.header["BMAJ"] * u.deg * pad - xlo = Longitude(cutout_args.ra_low * u.deg) - Longitude(padder) - xhi = Longitude(cutout_args.ra_high * u.deg) + Longitude(padder) - ylo = Latitude(cutout_args.dec_low * u.deg) - Latitude(padder) - yhi = Latitude(cutout_args.dec_high * u.deg) + Latitude(padder) + top_right = SkyCoord(cutout_args.ra_low * u.deg, cutout_args.dec_high * u.deg) + bottom_left = SkyCoord(cutout_args.ra_high * u.deg, cutout_args.dec_low * u.deg) - xp_lo, yp_lo = skycoord_to_pixel(SkyCoord(xlo, ylo), cube.wcs) - xp_hi, yp_hi = skycoord_to_pixel(SkyCoord(xhi, yhi), cube.wcs) + top_right_off = top_right.spherical_offsets_by(d_lon=padder, d_lat=padder) + bottom_left_off = bottom_left.spherical_offsets_by(d_lon=-padder, d_lat=-padder) + # Only need the critical corners - but just in case: + # top_left = SkyCoord(cutout_args.ra_high * u.deg, cutout_args.dec_high * u.deg) + # bottom_right = SkyCoord(cutout_args.ra_low, cutout_args.dec_low * u.deg) + # top_left_off = top_left.spherical_offsets_by(d_lon=-padder, d_lat=padder) + # bottom_right_off = bottom_right.spherical_offsets_by(d_lon=padder, d_lat=-padder) + + x_left, y_bottom = skycoord_to_pixel(bottom_left_off, cube.wcs) + x_right, y_top = skycoord_to_pixel(top_right_off, cube.wcs) # Round for cutout - yp_lo_idx = int(np.floor(yp_lo)) - yp_hi_idx = int(np.ceil(yp_hi)) - xp_lo_idx = int(np.floor(xp_hi)) - xp_hi_idx = int(np.ceil(xp_lo)) + yp_lo_idx = int(np.floor(y_bottom)) + yp_hi_idx = int(np.ceil(y_top)) + xp_lo_idx = int(np.floor(x_right)) + xp_hi_idx = int(np.ceil(x_left)) # Use subcube for header transformation cutout_cube = cube[:, yp_lo_idx:yp_hi_idx, xp_lo_idx:xp_hi_idx] @@ -237,32 +243,75 @@ def get_args( ra_max = np.max(coords.ra) ra_i_max = np.argmax(coords.ra) ra_off = Longitude(majs[ra_i_max]) - ra_high = ra_max + ra_off ra_min = np.min(coords.ra) ra_i_min = np.argmin(coords.ra) ra_off = Longitude(majs[ra_i_min]) - ra_low = ra_min - ra_off dec_max = np.max(coords.dec) dec_i_max = np.argmax(coords.dec) - dec_off = Longitude(majs[dec_i_max]) - dec_high = dec_max + dec_off + dec_off = Latitude(majs[dec_i_max]) dec_min = np.min(coords.dec) dec_i_min = np.argmin(coords.dec) - dec_off = Longitude(majs[dec_i_min]) - dec_low = dec_min - dec_off + dec_off = Latitude(majs[dec_i_min]) + + # Use SkyCoords to account for poles and meridian + top_right = SkyCoord(ra_min, dec_max) + top_left = SkyCoord(ra_max, dec_max) + bottom_right = SkyCoord(ra_min, dec_min) + bottom_left = SkyCoord(ra_max, dec_min) + top_right_off = top_right.spherical_offsets_by(d_lon=ra_off, d_lat=dec_off) + top_left_off = top_left.spherical_offsets_by(d_lon=-ra_off, d_lat=dec_off) + bottom_right_off = bottom_right.spherical_offsets_by( + d_lon=ra_off, d_lat=-dec_off + ) + bottom_left_off = bottom_left.spherical_offsets_by( + d_lon=-ra_off, d_lat=-dec_off + ) + ra_high: float = np.max( + [ + top_right_off.ra.deg, + top_left_off.ra.deg, + bottom_right_off.ra.deg, + bottom_left_off.ra.deg, + ] + ) + ra_low: float = np.min( + [ + top_right_off.ra.deg, + top_left_off.ra.deg, + bottom_right_off.ra.deg, + bottom_left_off.ra.deg, + ] + ) + dec_high: float = np.max( + [ + top_right_off.dec.deg, + top_left_off.dec.deg, + bottom_right_off.dec.deg, + bottom_left_off.dec.deg, + ] + ) + dec_low: float = np.min( + [ + top_right_off.dec.deg, + top_left_off.dec.deg, + bottom_right_off.dec.deg, + bottom_left_off.dec.deg, + ] + ) + except Exception as e: logger.debug(f"coords are {coords=}") logger.debug(f"comps are {comps=}") raise e return CutoutArgs( - ra_high=ra_high.deg, - ra_low=ra_low.deg, - dec_high=dec_high.deg, - dec_low=dec_low.deg, + ra_high=ra_high, + ra_low=ra_low, + dec_high=dec_high, + dec_low=dec_low, outdir=outdir, ) From 3ad92c54b0115251b0ba5b9fdd7be843df6102a0 Mon Sep 17 00:00:00 2001 From: "Thomson, Alec (CASS, Kensington)" Date: Sun, 12 May 2024 15:22:10 +1000 Subject: [PATCH 2/6] Sanity checks --- arrakis/cutout.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/arrakis/cutout.py b/arrakis/cutout.py index f6fe5479..02910db6 100644 --- a/arrakis/cutout.py +++ b/arrakis/cutout.py @@ -22,6 +22,7 @@ from astropy.utils.exceptions import AstropyWarning from astropy.wcs.utils import skycoord_to_pixel from prefect import flow, task +from prefect.task_runners import ConcurrentTaskRunner from spectral_cube import SpectralCube from spectral_cube.utils import SpectralCubeWarning from tqdm.auto import tqdm @@ -153,20 +154,22 @@ def cutout_image( top_right_off = top_right.spherical_offsets_by(d_lon=padder, d_lat=padder) bottom_left_off = bottom_left.spherical_offsets_by(d_lon=-padder, d_lat=-padder) - # Only need the critical corners - but just in case: - # top_left = SkyCoord(cutout_args.ra_high * u.deg, cutout_args.dec_high * u.deg) - # bottom_right = SkyCoord(cutout_args.ra_low, cutout_args.dec_low * u.deg) - # top_left_off = top_left.spherical_offsets_by(d_lon=-padder, d_lat=padder) - # bottom_right_off = bottom_right.spherical_offsets_by(d_lon=padder, d_lat=-padder) + + top_left = SkyCoord(cutout_args.ra_high * u.deg, cutout_args.dec_high * u.deg) + bottom_right = SkyCoord(cutout_args.ra_low * u.deg, cutout_args.dec_low * u.deg) + top_left_off = top_left.spherical_offsets_by(d_lon=-padder, d_lat=padder) + bottom_right_off = bottom_right.spherical_offsets_by(d_lon=padder, d_lat=-padder) x_left, y_bottom = skycoord_to_pixel(bottom_left_off, cube.wcs) x_right, y_top = skycoord_to_pixel(top_right_off, cube.wcs) + _x_left, _y_top = skycoord_to_pixel(top_left_off, cube.wcs) + _x_right, _y_bottom = skycoord_to_pixel(bottom_right_off, cube.wcs) - # Round for cutout - yp_lo_idx = int(np.floor(y_bottom)) - yp_hi_idx = int(np.ceil(y_top)) - xp_lo_idx = int(np.floor(x_right)) - xp_hi_idx = int(np.ceil(x_left)) + # Compare all points in case of insanity at the poles + yp_lo_idx = int(np.floor(min(y_bottom, _y_bottom, y_top, _y_top))) + yp_hi_idx = int(np.ceil(max(y_bottom, _y_bottom, y_top, _y_top))) + xp_lo_idx = int(np.floor(min(x_left, x_right, _x_left, _x_right))) + xp_hi_idx = int(np.ceil(max(x_left, x_right, _x_left, _x_right))) # Use subcube for header transformation cutout_cube = cube[:, yp_lo_idx:yp_hi_idx, xp_lo_idx:xp_hi_idx] @@ -572,7 +575,7 @@ def main(args: argparse.Namespace) -> None: args (argparse.Namespace): Command-line args verbose (bool, optional): Verbose output. Defaults to True. """ - cutout_islands( + cutout_islands.with_options(task_runner=ConcurrentTaskRunner)( field=args.field, directory=args.datadir, host=args.host, From 5a8e8cb590e26a63c2adcd36ef8fef54366e62d4 Mon Sep 17 00:00:00 2001 From: "Thomson, Alec (CASS, Kensington)" Date: Sun, 12 May 2024 17:55:48 +1000 Subject: [PATCH 3/6] Use celestial --- arrakis/cutout.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/arrakis/cutout.py b/arrakis/cutout.py index 02910db6..aab09151 100644 --- a/arrakis/cutout.py +++ b/arrakis/cutout.py @@ -160,10 +160,10 @@ def cutout_image( top_left_off = top_left.spherical_offsets_by(d_lon=-padder, d_lat=padder) bottom_right_off = bottom_right.spherical_offsets_by(d_lon=padder, d_lat=-padder) - x_left, y_bottom = skycoord_to_pixel(bottom_left_off, cube.wcs) - x_right, y_top = skycoord_to_pixel(top_right_off, cube.wcs) - _x_left, _y_top = skycoord_to_pixel(top_left_off, cube.wcs) - _x_right, _y_bottom = skycoord_to_pixel(bottom_right_off, cube.wcs) + x_left, y_bottom = skycoord_to_pixel(bottom_left_off, cube.wcs.celestial) + x_right, y_top = skycoord_to_pixel(top_right_off, cube.wcs.celestial) + _x_left, _y_top = skycoord_to_pixel(top_left_off, cube.wcs.celestial) + _x_right, _y_bottom = skycoord_to_pixel(bottom_right_off, cube.wcs.celestial) # Compare all points in case of insanity at the poles yp_lo_idx = int(np.floor(min(y_bottom, _y_bottom, y_top, _y_top))) From 65a5efd7adb4b9ba591a056ec100cb3597abefde Mon Sep 17 00:00:00 2001 From: "Thomson, Alec (CASS, Kensington)" Date: Mon, 13 May 2024 00:31:19 +1000 Subject: [PATCH 4/6] East vs West --- arrakis/cutout.py | 144 ++++++++++++++++++++------------------ arrakis/rmsynth_oncuts.py | 9 ++- 2 files changed, 81 insertions(+), 72 deletions(-) diff --git a/arrakis/cutout.py b/arrakis/cutout.py index aab09151..712b8d49 100644 --- a/arrakis/cutout.py +++ b/arrakis/cutout.py @@ -16,7 +16,7 @@ import numpy as np import pandas as pd import pymongo -from astropy.coordinates import Latitude, Longitude, SkyCoord +from astropy.coordinates import SkyCoord from astropy.io import fits from astropy.utils import iers from astropy.utils.exceptions import AstropyWarning @@ -56,9 +56,9 @@ class CutoutArgs(Struct): """Arguments for cutout function""" """Name of the source""" - ra_high: float + ra_left: float """Upper RA bound in degrees""" - ra_low: float + ra_right: float """Lower RA bound in degrees""" dec_high: float """Upper DEC bound in degrees""" @@ -149,16 +149,27 @@ def cutout_image( padder: float = cube.header["BMAJ"] * u.deg * pad - top_right = SkyCoord(cutout_args.ra_low * u.deg, cutout_args.dec_high * u.deg) - bottom_left = SkyCoord(cutout_args.ra_high * u.deg, cutout_args.dec_low * u.deg) + top_right = SkyCoord(cutout_args.ra_right * u.deg, cutout_args.dec_high * u.deg) + bottom_left = SkyCoord(cutout_args.ra_left * u.deg, cutout_args.dec_low * u.deg) + top_left = SkyCoord(cutout_args.ra_left * u.deg, cutout_args.dec_high * u.deg) + bottom_right = SkyCoord(cutout_args.ra_right * u.deg, cutout_args.dec_low * u.deg) - top_right_off = top_right.spherical_offsets_by(d_lon=padder, d_lat=padder) - bottom_left_off = bottom_left.spherical_offsets_by(d_lon=-padder, d_lat=-padder) - - top_left = SkyCoord(cutout_args.ra_high * u.deg, cutout_args.dec_high * u.deg) - bottom_right = SkyCoord(cutout_args.ra_low * u.deg, cutout_args.dec_low * u.deg) - top_left_off = top_left.spherical_offsets_by(d_lon=-padder, d_lat=padder) - bottom_right_off = bottom_right.spherical_offsets_by(d_lon=padder, d_lat=-padder) + top_right_off = top_right.directional_offset_by( + position_angle=-45 * u.deg, + separation=padder * np.sqrt(2), + ) + bottom_left_off = bottom_left.directional_offset_by( + position_angle=135 * u.deg, + separation=padder * np.sqrt(2), + ) + top_left_off = top_left.directional_offset_by( + position_angle=45 * u.deg, + separation=padder * np.sqrt(2), + ) + bottom_right_off = bottom_right.directional_offset_by( + position_angle=-135 * u.deg, + separation=padder * np.sqrt(2), + ) x_left, y_bottom = skycoord_to_pixel(bottom_left_off, cube.wcs.celestial) x_right, y_top = skycoord_to_pixel(top_right_off, cube.wcs.celestial) @@ -239,71 +250,64 @@ def get_args( ras: u.Quantity = comps.RA.values * u.deg decs: u.Quantity = comps.Dec.values * u.deg majs: List[float] = comps.Maj.values * u.arcsec + coords = SkyCoord(ras, decs, frame="icrs") + padder = np.max(majs) coords = SkyCoord(ras, decs) try: - ra_max = np.max(coords.ra) - ra_i_max = np.argmax(coords.ra) - ra_off = Longitude(majs[ra_i_max]) - - ra_min = np.min(coords.ra) - ra_i_min = np.argmin(coords.ra) - ra_off = Longitude(majs[ra_i_min]) - - dec_max = np.max(coords.dec) - dec_i_max = np.argmax(coords.dec) - dec_off = Latitude(majs[dec_i_max]) - - dec_min = np.min(coords.dec) - dec_i_min = np.argmin(coords.dec) - dec_off = Latitude(majs[dec_i_min]) + northern_boundary = SkyCoord(ras, np.ones_like(ras.value) * decs.max()) + for i, coord_1 in enumerate(northern_boundary): + for j, coord_2 in enumerate(northern_boundary): + sep = coord_1.separation(coord_2) + pa = coord_1.position_angle(coord_2) + if i == 0 and j == 0: + largest_sep = sep + largest_pa = pa + largest_i = i + largest_j = j + if sep > largest_sep: + largest_i = i + largest_j = j + largest_sep = sep + largest_pa = pa + + if largest_pa > 90 * u.deg: + left_point = northern_boundary[largest_i] + right_point = northern_boundary[largest_j] + else: + left_point = northern_boundary[largest_j] + right_point = northern_boundary[largest_i] # Use SkyCoords to account for poles and meridian - top_right = SkyCoord(ra_min, dec_max) - top_left = SkyCoord(ra_max, dec_max) - bottom_right = SkyCoord(ra_min, dec_min) - bottom_left = SkyCoord(ra_max, dec_min) - top_right_off = top_right.spherical_offsets_by(d_lon=ra_off, d_lat=dec_off) - top_left_off = top_left.spherical_offsets_by(d_lon=-ra_off, d_lat=dec_off) - bottom_right_off = bottom_right.spherical_offsets_by( - d_lon=ra_off, d_lat=-dec_off + top_right = SkyCoord(right_point.ra, decs.max()) + top_left = SkyCoord(left_point.ra, decs.max()) + bottom_right = SkyCoord(right_point.ra, decs.min()) + bottom_left = SkyCoord(left_point.ra, decs.min()) + + top_right_off = top_right.directional_offset_by( + position_angle=-45 * u.deg, + separation=padder * np.sqrt(2), ) - bottom_left_off = bottom_left.spherical_offsets_by( - d_lon=-ra_off, d_lat=-dec_off + bottom_left_off = bottom_left.directional_offset_by( + position_angle=135 * u.deg, + separation=padder * np.sqrt(2), ) - ra_high: float = np.max( - [ - top_right_off.ra.deg, - top_left_off.ra.deg, - bottom_right_off.ra.deg, - bottom_left_off.ra.deg, - ] + top_left_off = top_left.directional_offset_by( + position_angle=45 * u.deg, + separation=padder * np.sqrt(2), ) - ra_low: float = np.min( - [ - top_right_off.ra.deg, - top_left_off.ra.deg, - bottom_right_off.ra.deg, - bottom_left_off.ra.deg, - ] + bottom_right_off = bottom_right.directional_offset_by( + position_angle=-135 * u.deg, + separation=padder * np.sqrt(2), ) - dec_high: float = np.max( - [ - top_right_off.dec.deg, - top_left_off.dec.deg, - bottom_right_off.dec.deg, - bottom_left_off.dec.deg, - ] - ) - dec_low: float = np.min( - [ - top_right_off.dec.deg, - top_left_off.dec.deg, - bottom_right_off.dec.deg, - bottom_left_off.dec.deg, - ] + + dec_low = ( + np.min([bottom_left_off.dec.value, bottom_right_off.dec.value]) * u.deg ) + dec_high = np.max([top_left_off.dec.value, top_right_off.dec.value]) * u.deg + ra_right = top_right_off.ra + ra_left = top_left_off.ra except Exception as e: logger.debug(f"coords are {coords=}") @@ -311,10 +315,10 @@ def get_args( raise e return CutoutArgs( - ra_high=ra_high, - ra_low=ra_low, - dec_high=dec_high, - dec_low=dec_low, + ra_left=ra_left.to(u.deg).value, + ra_right=ra_right.to(u.deg).value, + dec_high=dec_high.to(u.deg).value, + dec_low=dec_low.to(u.deg).value, outdir=outdir, ) diff --git a/arrakis/rmsynth_oncuts.py b/arrakis/rmsynth_oncuts.py index 6a669b35..2162e6ab 100644 --- a/arrakis/rmsynth_oncuts.py +++ b/arrakis/rmsynth_oncuts.py @@ -26,6 +26,7 @@ from astropy.wcs import WCS from astropy.wcs.utils import proj_plane_pixel_scales from prefect import flow, task +from prefect.task_runners import ConcurrentTaskRunner from radio_beam import Beam from RMtools_1D import do_RMsynth_1D from RMtools_3D import do_RMsynth_3D @@ -296,7 +297,11 @@ def extract_single_spectrum( rms[np.isnan(rms)] = np.nanmedian(rms) wcs = WCS(header) x, y = np.array(wcs.celestial.world_to_pixel(coord)).round().astype(int) - spectrum_arr = np.array(data[:, y, x]) + try: + spectrum_arr = np.array(data[:, y, x]) + except Exception as e: + logger.error(f"Error extracting spectrum from {filename}") + raise e spectrum_arr[spectrum_arr == 0] = np.nan rms[~np.isfinite(spectrum_arr)] = np.nan bkg[~np.isfinite(spectrum_arr)] = np.nan @@ -1205,7 +1210,7 @@ def cli(): password=args.password, ) - main( + main.with_options(task_runner=ConcurrentTaskRunner)( field=args.field, outdir=Path(args.datadir), host=args.host, From ce1db3e9808cc97a3185a29f22d070c7b2bb601f Mon Sep 17 00:00:00 2001 From: "Thomson, Alec (CASS, Kensington)" Date: Mon, 13 May 2024 11:06:49 +1000 Subject: [PATCH 5/6] Improve boundary checker --- arrakis/cutout.py | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/arrakis/cutout.py b/arrakis/cutout.py index 712b8d49..5f9ec62c 100644 --- a/arrakis/cutout.py +++ b/arrakis/cutout.py @@ -256,9 +256,23 @@ def get_args( coords = SkyCoord(ras, decs) try: + # If in South - Northern boundary will have greatest RA range + # And vice versa northern_boundary = SkyCoord(ras, np.ones_like(ras.value) * decs.max()) - for i, coord_1 in enumerate(northern_boundary): - for j, coord_2 in enumerate(northern_boundary): + southern_boundary = SkyCoord(ras, np.ones_like(ras.value) * decs.min()) + + if abs(decs.max()) > abs(decs.min()): + # North - use Southern + longest_boundary = southern_boundary + else: + # South - use Northern + longest_boundary = northern_boundary + + # Find the coords with greatest separation + # That gives us our image size + # PA tells us East from West wrap + for i, coord_1 in enumerate(longest_boundary): + for j, coord_2 in enumerate(longest_boundary): sep = coord_1.separation(coord_2) pa = coord_1.position_angle(coord_2) if i == 0 and j == 0: @@ -272,12 +286,14 @@ def get_args( largest_sep = sep largest_pa = pa - if largest_pa > 90 * u.deg: - left_point = northern_boundary[largest_i] - right_point = northern_boundary[largest_j] + if largest_pa > 180 * u.deg: + # PA points West / left + left_point = longest_boundary[largest_i] + right_point = longest_boundary[largest_j] else: - left_point = northern_boundary[largest_j] - right_point = northern_boundary[largest_i] + # PA points East / right + left_point = longest_boundary[largest_j] + right_point = longest_boundary[largest_i] # Use SkyCoords to account for poles and meridian top_right = SkyCoord(right_point.ra, decs.max()) From d2ff4f7e566e3e4aa53472181787e565618e6014 Mon Sep 17 00:00:00 2001 From: "Thomson, Alec (CASS, Kensington)" Date: Mon, 13 May 2024 11:59:27 +1000 Subject: [PATCH 6/6] Use numpy --- arrakis/cutout.py | 61 +++++++++++++++++++---------------------------- 1 file changed, 24 insertions(+), 37 deletions(-) diff --git a/arrakis/cutout.py b/arrakis/cutout.py index 5f9ec62c..4f37e090 100644 --- a/arrakis/cutout.py +++ b/arrakis/cutout.py @@ -253,75 +253,62 @@ def get_args( coords = SkyCoord(ras, decs, frame="icrs") padder = np.max(majs) - coords = SkyCoord(ras, decs) - try: # If in South - Northern boundary will have greatest RA range # And vice versa northern_boundary = SkyCoord(ras, np.ones_like(ras.value) * decs.max()) southern_boundary = SkyCoord(ras, np.ones_like(ras.value) * decs.min()) - if abs(decs.max()) > abs(decs.min()): + if np.abs(decs.max()) > np.abs(decs.min()): # North - use Southern longest_boundary = southern_boundary else: # South - use Northern longest_boundary = northern_boundary - # Find the coords with greatest separation - # That gives us our image size - # PA tells us East from West wrap - for i, coord_1 in enumerate(longest_boundary): - for j, coord_2 in enumerate(longest_boundary): - sep = coord_1.separation(coord_2) - pa = coord_1.position_angle(coord_2) - if i == 0 and j == 0: - largest_sep = sep - largest_pa = pa - largest_i = i - largest_j = j - if sep > largest_sep: - largest_i = i - largest_j = j - largest_sep = sep - largest_pa = pa - - if largest_pa > 180 * u.deg: - # PA points West / left + # Compute separations and position angles between all pairs of coordinates + separations = longest_boundary[:, np.newaxis].separation(longest_boundary) + position_angles = longest_boundary[:, np.newaxis].position_angle( + longest_boundary + ) + + # Find the index of the pair with the largest separation + largest_i, largest_j = np.unravel_index( + np.argmax(separations), separations.shape + ) + + # Determine left and right points based on position angle + if position_angles[largest_i, largest_j] > 180 * u.deg: left_point = longest_boundary[largest_i] right_point = longest_boundary[largest_j] else: - # PA points East / right left_point = longest_boundary[largest_j] right_point = longest_boundary[largest_i] - # Use SkyCoords to account for poles and meridian + # Compute coordinates for top right, top left, bottom right, and bottom left points top_right = SkyCoord(right_point.ra, decs.max()) top_left = SkyCoord(left_point.ra, decs.max()) bottom_right = SkyCoord(right_point.ra, decs.min()) bottom_left = SkyCoord(left_point.ra, decs.min()) + # Compute offsets for top right, top left, bottom right, and bottom left points + offset = padder * np.sqrt(2) top_right_off = top_right.directional_offset_by( - position_angle=-45 * u.deg, - separation=padder * np.sqrt(2), + position_angle=-45 * u.deg, separation=offset ) bottom_left_off = bottom_left.directional_offset_by( - position_angle=135 * u.deg, - separation=padder * np.sqrt(2), + position_angle=135 * u.deg, separation=offset ) top_left_off = top_left.directional_offset_by( - position_angle=45 * u.deg, - separation=padder * np.sqrt(2), + position_angle=45 * u.deg, separation=offset ) bottom_right_off = bottom_right.directional_offset_by( - position_angle=-135 * u.deg, - separation=padder * np.sqrt(2), + position_angle=-135 * u.deg, separation=offset ) - dec_low = ( - np.min([bottom_left_off.dec.value, bottom_right_off.dec.value]) * u.deg - ) - dec_high = np.max([top_left_off.dec.value, top_right_off.dec.value]) * u.deg + # Compute dec_low, dec_high, ra_right, and ra_left + dec_low = min(bottom_left_off.dec.value, bottom_right_off.dec.value) * u.deg + dec_high = max(top_left_off.dec.value, top_right_off.dec.value) * u.deg ra_right = top_right_off.ra ra_left = top_left_off.ra