From b176baafba513a6e2b1ee51b9f03ef464b703a67 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 28 Jun 2023 16:52:42 +0200 Subject: [PATCH 1/2] extract from imagehdu so that full fov is covered --- scopesim/optics/fov_utils.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/scopesim/optics/fov_utils.py b/scopesim/optics/fov_utils.py index 43684c6e..62ed362f 100644 --- a/scopesim/optics/fov_utils.py +++ b/scopesim/optics/fov_utils.py @@ -273,7 +273,7 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): Parameters ---------- imagehdu : fits.ImageHDU - The field ImageHDU, either an image of a wavelength [um] cube + The field ImageHDU, either an image or a cube with wavelength [um] fov_volume : dict Contains {"xs": [xmin, xmax], "ys": [ymin, ymax], "waves": [wave_min, wave_max], @@ -294,7 +294,11 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): y0s, y1s = max(min(y_hdu), min(y_fov)), min(max(y_hdu), max(y_fov)) xp, yp = imp_utils.val2pix(hdr, np.array([x0s, x1s]), np.array([y0s, y1s])) - (x0p, x1p), (y0p, y1p) = np.round(xp).astype(int), np.round(yp).astype(int) + x0p = np.floor(xp[0]).astype(int) + x1p = np.ceil(xp[1]).astype(int) + y0p = np.floor(yp[0]).astype(int) + y1p = np.ceil(yp[1]).astype(int) + # (x0p, x1p), (y0p, y1p) = np.round(xp).astype(int), np.round(yp).astype(int) if x0p == x1p: x1p += 1 if y0p == y1p: @@ -331,7 +335,7 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): i0p, i1p = np.where(mask)[0][0], np.where(mask)[0][-1] f0 = (abs(hdu_waves[i0p] - fov_waves[0] + 0.5 * wdel) % wdel) / wdel # blue edge f1 = (abs(hdu_waves[i1p] - fov_waves[1] - 0.5 * wdel) % wdel) / wdel # red edge - data = imagehdu.data[i0p:i1p+1, y0p:y1p, x0p:x1p] + data = imagehdu.data[i0p:i1p+1, y0p:y1p+1, x0p:x1p+1] data[0, :, :] *= f0 if i1p > i0p: data[-1, :, :] *= f1 @@ -353,7 +357,7 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): "BUNIT": hdr["BUNIT"]}) else: - data = imagehdu.data[y0p:y1p, x0p:x1p] + data = imagehdu.data[y0p:y1p+1, x0p:x1p+1] new_hdr["SPEC_REF"] = hdr.get("SPEC_REF") new_imagehdu = fits.ImageHDU(data=data) From 4566eea98d6c4a4eba3eb835c5cc00daed4e5a7a Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 28 Jun 2023 22:39:19 +0200 Subject: [PATCH 2/2] Ensure pixels in range --- scopesim/optics/fov_utils.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/scopesim/optics/fov_utils.py b/scopesim/optics/fov_utils.py index 62ed362f..8c554b49 100644 --- a/scopesim/optics/fov_utils.py +++ b/scopesim/optics/fov_utils.py @@ -286,7 +286,7 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): """ hdr = imagehdu.header new_hdr = {} - + naxis1, naxis2 = hdr["NAXIS1"], hdr["NAXIS2"] x_hdu, y_hdu = imp_utils.calc_footprint(imagehdu) # field edges in "deg" x_fov, y_fov = fov_volume["xs"], fov_volume["ys"] @@ -294,10 +294,10 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): y0s, y1s = max(min(y_hdu), min(y_fov)), min(max(y_hdu), max(y_fov)) xp, yp = imp_utils.val2pix(hdr, np.array([x0s, x1s]), np.array([y0s, y1s])) - x0p = np.floor(xp[0]).astype(int) - x1p = np.ceil(xp[1]).astype(int) - y0p = np.floor(yp[0]).astype(int) - y1p = np.ceil(yp[1]).astype(int) + x0p = max(0, np.floor(xp[0]).astype(int)) + x1p = min(naxis1, np.ceil(xp[1]).astype(int)) + y0p = max(0, np.floor(yp[0]).astype(int)) + y1p = min(naxis2, np.ceil(yp[1]).astype(int)) # (x0p, x1p), (y0p, y1p) = np.round(xp).astype(int), np.round(yp).astype(int) if x0p == x1p: x1p += 1 @@ -335,7 +335,7 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): i0p, i1p = np.where(mask)[0][0], np.where(mask)[0][-1] f0 = (abs(hdu_waves[i0p] - fov_waves[0] + 0.5 * wdel) % wdel) / wdel # blue edge f1 = (abs(hdu_waves[i1p] - fov_waves[1] - 0.5 * wdel) % wdel) / wdel # red edge - data = imagehdu.data[i0p:i1p+1, y0p:y1p+1, x0p:x1p+1] + data = imagehdu.data[i0p:i1p+1, y0p:y1p, x0p:x1p] data[0, :, :] *= f0 if i1p > i0p: data[-1, :, :] *= f1 @@ -357,7 +357,7 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): "BUNIT": hdr["BUNIT"]}) else: - data = imagehdu.data[y0p:y1p+1, x0p:x1p+1] + data = imagehdu.data[y0p:y1p, x0p:x1p] new_hdr["SPEC_REF"] = hdr.get("SPEC_REF") new_imagehdu = fits.ImageHDU(data=data)