From 7b5ec5c77ccd7a209f95d0b6c08005663dde8c99 Mon Sep 17 00:00:00 2001 From: Sam Van Kooten Date: Fri, 25 Aug 2023 18:28:18 -0600 Subject: [PATCH] Use more points to find image bounds in moasics --- reproject/mosaicking/coadd.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py index efe2089c4..0014d1245 100644 --- a/reproject/mosaicking/coadd.py +++ b/reproject/mosaicking/coadd.py @@ -125,16 +125,29 @@ def reproject_and_coadd( # Since we might be reprojecting small images into a large mosaic we # want to make sure that for each image we reproject to an array with - # minimal footprint. We therefore find the pixel coordinates of corners - # in the initial image and transform this to pixel coordinates in the - # final image to figure out the final WCS and shape to reproject to for - # each tile. Note that in future if we are worried about significant - # distortions of the edges in the reprojection process we could simply - # add arbitrary numbers of midpoints to this list. + # minimal footprint. We therefore find the pixel coordinates of the + # edges of the initial image and transform this to pixel coordinates in + # the final image to figure out the final WCS and shape to reproject to + # for each tile. We strike a balance between transforming only the + # input-image corners, which is fast but can cause clipping in cases of + # significant distortion (when the edges of the input image become + # convex in the output projection), and transforming every edge pixel, + # which provides a lot of redundant information. ny, nx = array_in.shape - xc = np.array([-0.5, nx - 0.5, nx - 0.5, -0.5]) - yc = np.array([-0.5, -0.5, ny - 0.5, ny - 0.5]) - xc_out, yc_out = wcs_out.world_to_pixel(wcs_in.pixel_to_world(xc, yc)) + n_per_edge = 10 + xs = np.linspace(-0.5, nx - 0.5, n_per_edge) + ys = np.linspace(-0.5, ny - 0.5, n_per_edge) + xs = np.concatenate(( + xs, + np.full(n_per_edge, xs[-1]), + xs, + np.full(n_per_edge, xs[0]))) + ys = np.concatenate(( + np.full(n_per_edge, ys[0]), + ys, + np.full(n_per_edge, ys[-1]), + ys)) + xc_out, yc_out = wcs_out.world_to_pixel(wcs_in.pixel_to_world(xs, ys)) # Determine the cutout parameters