From ddcbcce00b1bcb0ca749b33fd99fe2fce7935d08 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 31 Jan 2024 22:36:32 -0500 Subject: [PATCH 1/4] generalize nrc_ta comparison code to work on any NIRCam TA, not just WFS. Add more careful bad pix handling using DQ mask. --- webbpsf/mast_wss.py | 2 +- webbpsf/match_data.py | 5 +++-- webbpsf/trending.py | 30 ++++++++++++++++++------------ 3 files changed, 22 insertions(+), 15 deletions(-) diff --git a/webbpsf/mast_wss.py b/webbpsf/mast_wss.py index 11e36358..47e837c2 100644 --- a/webbpsf/mast_wss.py +++ b/webbpsf/mast_wss.py @@ -598,7 +598,7 @@ def get_visit_nrc_ta_image(visitid, verbose=True): from astroquery.mast import Mast keywords = { 'visit_id': [visitid[1:]], # note: drop the initial character 'V' - 'category': ['CAL'], + #'category': ['CAL'], 'exp_type': ['NRC_TACQ'] } diff --git a/webbpsf/match_data.py b/webbpsf/match_data.py index cbf68961..72290c9c 100644 --- a/webbpsf/match_data.py +++ b/webbpsf/match_data.py @@ -49,8 +49,9 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic if inst.name == 'NIRCam': if header['PUPIL'].startswith('MASK'): inst.pupil_mask = header['PUPIL'] - inst.image_mask = header['CORONMSK'].replace('MASKA', 'MASK') # note, have to modify the value slightly for - # consistency with the labels used in webbpsf + if 'CORONMSK' in header: + inst.image_mask = header['CORONMSK'].replace('MASKA', 'MASK') # note, have to modify the value slightly for + # consistency with the labels used in webbpsf inst.set_position_from_aperture_name(header['APERNAME']) # Redo this, in case the image_mask setting auto switched it to something else elif inst.name == 'MIRI': if inst.filter in ['F1065C', 'F1140C', 'F1550C']: diff --git a/webbpsf/trending.py b/webbpsf/trending.py index ea85b2a6..9cae611c 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1619,7 +1619,7 @@ def vprint(*args, **kwargs): #### Functions for image comparisons -def show_wfs_ta_img(visitid, ax=None, return_handles=False): +def show_nrc_ta_img(visitid, ax=None, return_handles=False): """ Retrieve and display a WFS target acq image""" hdul = webbpsf.mast_wss.get_visit_nrc_ta_image(visitid) @@ -1633,15 +1633,13 @@ def show_wfs_ta_img(visitid, ax=None, return_handles=False): cmap = matplotlib.cm.viridis.copy() cmap.set_bad('orange') - norm = matplotlib.colors.AsinhNorm(linear_width = vmax*0.003, vmax=vmax, #vmin=0) vmin=-1*rsig) - if ax is None: ax = plt.gca() ax.imshow(ta_img - bglevel, norm=norm, cmap=cmap, origin='lower') - ax.set_title(f"WFS TA on {visitid}\n{hdul[0].header['DATE-OBS']}") + ax.set_title(f"NIRCam TA on {visitid}\n{hdul[0].header['DATE-OBS']}") ax.set_ylabel("[Pixels]") ax.text(0.05, 0.9, hdul[0].header['TARGPROP'], color='white', transform=ax.transAxes) @@ -1650,7 +1648,7 @@ def show_wfs_ta_img(visitid, ax=None, return_handles=False): if return_handles: return hdul, ax, norm, cmap, bglevel -def nrc_ta_image_comparison(visitid): +def nrc_ta_image_comparison(visitid, verbose=False): """ Retrieve a NIRCam target acq image and compare to a simulation Parameters: @@ -1664,26 +1662,34 @@ def nrc_ta_image_comparison(visitid): fig, axes = plt.subplots(figsize=(10,5), ncols=3) # Get and plot the observed TA image - hdul, ax, norm, cmap, bglevel = show_wfs_ta_img(visitid, ax=axes[0], return_handles=True) - im_obs = hdul['SCI'].data + hdul, ax, norm, cmap, bglevel = show_nrc_ta_img(visitid, ax=axes[0], return_handles=True) + im_obs = hdul['SCI'].data.copy() im_obs_err = hdul['ERR'].data + im_obs_dq = hdul['DQ'].data + + im_obs_clean = im_obs.copy() + im_obs_clean[im_obs_dq & 1] = np.nan # Mask out any DO_NOT_USE pixels. + im_obs_clean = astropy.convolution.interpolate_replace_nans(im_obs_clean, kernel=np.ones((5,5))) # Make a matching sim nrc = webbpsf.setup_sim_to_match_file(hdul, verbose=False) opdname = nrc.pupilopd[0].header['CORR_ID'] + "-NRCA3_FP1-1.fits" - psf = nrc.calc_psf(fov_pixels=64) + if verbose: + print(f"Calculating PSF to match that TA image...") + psf = nrc.calc_psf(fov_pixels=im_obs.shape[0]) # Align and Shift: - im_sim = psf['DET_DIST'].data # Use the extension including distortion and IPC - im_obs_clean = astropy.convolution.interpolate_replace_nans(im_obs, kernel=np.ones((5,5))) - shift, _, _ = phase_cross_correlation(im_obs_clean, im_sim, upsample_factor=32) + if verbose: + print(f"Shift to register sim to data: {shift} pix") im_sim_shifted = scipy.ndimage.shift(im_sim, shift, order=5) # figure out the background level and scale factor - scalefactor = np.nanmax(im_obs) / im_sim.max() + scalefactor = np.nanmax(im_obs_clean) / im_sim.max() + if verbose: + print(f"Scale factor to match sim to data: {scalefactor}") im_sim_scaled_aligned = im_sim_shifted*scalefactor From 0d5aaaa3a65e85e190eb5cb7ee397cdac364db07 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 31 Jan 2024 23:32:46 -0500 Subject: [PATCH 2/4] add optional code to overplot centroid on the TA image, including the OSS on-board measurement --- webbpsf/trending.py | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/webbpsf/trending.py b/webbpsf/trending.py index 9cae611c..60c826e9 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1644,11 +1644,10 @@ def show_nrc_ta_img(visitid, ax=None, return_handles=False): ax.text(0.05, 0.9, hdul[0].header['TARGPROP'], color='white', transform=ax.transAxes) - if return_handles: return hdul, ax, norm, cmap, bglevel -def nrc_ta_image_comparison(visitid, verbose=False): +def nrc_ta_image_comparison(visitid, verbose=False, show_centroid=False): """ Retrieve a NIRCam target acq image and compare to a simulation Parameters: @@ -1694,6 +1693,29 @@ def nrc_ta_image_comparison(visitid, verbose=False): im_sim_scaled_aligned = im_sim_shifted*scalefactor # Plot + if show_centroid: + # First, see if we can retrieve the on-board TA centroid measurment from the OSS engineering DB in MAST + try: + import misc_jwst # Optional dependency, including engineering database access tools + # If we have that, retrieve the log for this visit, extract the OSS centroids, and convert to same + # coordinate frame as used here: + osslog = misc_jwst.engdb.get_ictm_event_log(hdul[0].header['VSTSTART'], hdul[0].header['VISITEND']) + oss_cen = misc_jwst.engdb.extract_oss_TA_centroids(osslog, 'V' + hdul[0].header['VISIT_ID']) + # Convert from full-frame (as used by OSS) to detector subarray coords: + oss_cen_sci = nrc._detector_geom_info.aperture.det_to_sci( *oss_cen) + oss_cen_sci_pythonic = np.asarray(oss_cen_sci) - 1 # convert from 1-based pixel indexing to 0-based + oss_centroid_text = f"\n OSS centroid: {oss_cen_sci_pythonic[0]:.2f}, {oss_cen_sci_pythonic[1]:.2f}" + #print(oss_cen_sci_pythonic) + except ImportError: + oss_centroid_text = "" + + cen = webbpsf.fwcentroid.fwcentroid(im_obs_clean) + axes[0].scatter(cen[1], cen[0], color='red', marker='+', s=50) + axes[0].text(0.95, 0.05, f' Centroid: {cen[1]:.2f}, {cen[0]:.2f}'+oss_centroid_text, + horizontalalignment='right', verticalalignment='bottom', + transform=axes[0].transAxes, + color='white') + axes[1].imshow(im_sim_scaled_aligned, norm=norm, cmap=cmap, origin='lower') axes[1].set_title(f"Simulated PSF in F212N\nusing {opdname}") From 749e37fabeb65d7e3142e3b4366bdedeb263c7c9 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Mon, 5 Feb 2024 22:24:43 -0500 Subject: [PATCH 3/4] update TA code: allow proprietary TA data access, refactor plotting to show WCS, OSS, and webbpsf centroids. And PEP8 cleanup. --- webbpsf/mast_wss.py | 27 ++++++++++++++++++--- webbpsf/match_data.py | 4 ++-- webbpsf/trending.py | 55 +++++++++++++++++++++++++++++++++++-------- 3 files changed, 71 insertions(+), 15 deletions(-) diff --git a/webbpsf/mast_wss.py b/webbpsf/mast_wss.py index 47e837c2..d640770a 100644 --- a/webbpsf/mast_wss.py +++ b/webbpsf/mast_wss.py @@ -588,7 +588,7 @@ def get_corrections(opdtable): @functools.lru_cache -def get_visit_nrc_ta_image(visitid, verbose=True): +def get_visit_nrc_ta_image(visitid, verbose=True, kind='cal'): """Retrieve from MAST the NIRCam target acq image for a given visit. This retrieves an image from MAST and returns it as a HDUList variable @@ -598,7 +598,6 @@ def get_visit_nrc_ta_image(visitid, verbose=True): from astroquery.mast import Mast keywords = { 'visit_id': [visitid[1:]], # note: drop the initial character 'V' - #'category': ['CAL'], 'exp_type': ['NRC_TACQ'] } @@ -615,9 +614,31 @@ def set_params(parameters): t = Mast.service_request(service, params) filename = t[0]['filename'] + # If user manually specifies rate or uncal, retrieve that instead + if kind == 'rate' or kind == 'uncal': + filename = filename.replace('_cal.fits', f'_{kind}.fits') + if verbose: print(f"TA filename: {filename}") mast_file_url = f"https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:JWST/product/{filename}" - ta_hdul = fits.open(mast_file_url) + import urllib + + try: + ta_hdul = fits.open(mast_file_url) + except urllib.error.HTTPError as err: + if err.code == 401: # Unauthorized + # Use MAST API to allow retrieval of exclusive access data, if relevant + import astroquery + import tempfile + mast_api_token = os.environ.get('MAST_API_TOKEN', None) + mast_obs = astroquery.mast.ObservationsClass(mast_api_token) + uri = f"mast:JWST/product/{filename}" + with tempfile.NamedTemporaryFile() as temp: + mast_obs.download_file(uri, local_path=temp.name, cache=False) + ta_hdul = astropy.io.fits.open(temp.name) + else: + raise # re-raise any errors other than 401 for permissions denied + + return ta_hdul diff --git a/webbpsf/match_data.py b/webbpsf/match_data.py index 72290c9c..6dfa959f 100644 --- a/webbpsf/match_data.py +++ b/webbpsf/match_data.py @@ -50,8 +50,8 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic if header['PUPIL'].startswith('MASK'): inst.pupil_mask = header['PUPIL'] if 'CORONMSK' in header: - inst.image_mask = header['CORONMSK'].replace('MASKA', 'MASK') # note, have to modify the value slightly for - # consistency with the labels used in webbpsf + inst.image_mask = header['CORONMSK'].replace('MASKA', 'MASK') # note, have to modify the value slightly for + # consistency with the labels used in webbpsf inst.set_position_from_aperture_name(header['APERNAME']) # Redo this, in case the image_mask setting auto switched it to something else elif inst.name == 'MIRI': if inst.filter in ['F1065C', 'F1140C', 'F1550C']: diff --git a/webbpsf/trending.py b/webbpsf/trending.py index 60c826e9..f903e7ce 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1620,7 +1620,7 @@ def vprint(*args, **kwargs): #### Functions for image comparisons def show_nrc_ta_img(visitid, ax=None, return_handles=False): - """ Retrieve and display a WFS target acq image""" + """ Retrieve and display a NIRCam target acq image""" hdul = webbpsf.mast_wss.get_visit_nrc_ta_image(visitid) @@ -1647,6 +1647,7 @@ def show_nrc_ta_img(visitid, ax=None, return_handles=False): if return_handles: return hdul, ax, norm, cmap, bglevel + def nrc_ta_image_comparison(visitid, verbose=False, show_centroid=False): """ Retrieve a NIRCam target acq image and compare to a simulation @@ -1668,7 +1669,7 @@ def nrc_ta_image_comparison(visitid, verbose=False, show_centroid=False): im_obs_clean = im_obs.copy() im_obs_clean[im_obs_dq & 1] = np.nan # Mask out any DO_NOT_USE pixels. - im_obs_clean = astropy.convolution.interpolate_replace_nans(im_obs_clean, kernel=np.ones((5,5))) + im_obs_clean = astropy.convolution.interpolate_replace_nans(im_obs_clean, kernel=np.ones((5, 5))) # Make a matching sim nrc = webbpsf.setup_sim_to_match_file(hdul, verbose=False) @@ -1694,24 +1695,58 @@ def nrc_ta_image_comparison(visitid, verbose=False, show_centroid=False): # Plot if show_centroid: + ### OSS CENTROIDS ### # First, see if we can retrieve the on-board TA centroid measurment from the OSS engineering DB in MAST try: import misc_jwst # Optional dependency, including engineering database access tools # If we have that, retrieve the log for this visit, extract the OSS centroids, and convert to same # coordinate frame as used here: osslog = misc_jwst.engdb.get_ictm_event_log(hdul[0].header['VSTSTART'], hdul[0].header['VISITEND']) - oss_cen = misc_jwst.engdb.extract_oss_TA_centroids(osslog, 'V' + hdul[0].header['VISIT_ID']) - # Convert from full-frame (as used by OSS) to detector subarray coords: - oss_cen_sci = nrc._detector_geom_info.aperture.det_to_sci( *oss_cen) - oss_cen_sci_pythonic = np.asarray(oss_cen_sci) - 1 # convert from 1-based pixel indexing to 0-based - oss_centroid_text = f"\n OSS centroid: {oss_cen_sci_pythonic[0]:.2f}, {oss_cen_sci_pythonic[1]:.2f}" - #print(oss_cen_sci_pythonic) + try: + oss_cen = misc_jwst.engdb.extract_oss_TA_centroids(osslog, 'V' + hdul[0].header['VISIT_ID']) + # Convert from full-frame (as used by OSS) to detector subarray coords: + oss_cen_sci = nrc._detector_geom_info.aperture.det_to_sci(*oss_cen) + oss_cen_sci_pythonic = np.asarray(oss_cen_sci) - 1 # convert from 1-based pixel indexing to 0-based + oss_centroid_text = f"\n OSS centroid: {oss_cen_sci_pythonic[0]:.2f}, {oss_cen_sci_pythonic[1]:.2f}" + axes[0].scatter(oss_cen_sci_pythonic[0], oss_cen_sci_pythonic[1], color='0.5', marker='+', s=50) + axes[0].text(oss_cen_sci_pythonic[0], oss_cen_sci_pythonic[1], 'OSS ', color='0.9', verticalalignment='center', horizontalalignment='right') + if verbose: + print(f"OSS centroid on board: {oss_cen} (full det coord frame, 1-based)") + print(f"OSS centroid converted: {oss_cen_sci_pythonic} (sci frame in {nrc._detector_geom_info.aperture.AperName}, 0-based)") + full_ap = nrc.siaf[nrc._detector_geom_info.aperture.AperName[0:5] + "_FULL"] + oss_cen_full_sci = np.asarray(full_ap.det_to_sci(*oss_cen)) - 1 + print(f"OSS centroid converted: {oss_cen_full_sci} (sci frame in {full_ap.AperName}, 0-based)") + + except RuntimeError: + if verbose: + print("Could not parse TA coordinates from log. TA may have failed?") + oss_cen_sci_pythonic = None + + ### WCS COORDINATES ### + import jwst.datamodels + model = jwst.datamodels.open(hdul) + targ_coords = astropy.coordinates.SkyCoord(model.meta.target.ra, model.meta.target.dec, frame='icrs', unit=u.deg) + targ_coords_pix = model.meta.wcs.world_to_pixel(targ_coords) # returns x, y + axes[0].scatter(targ_coords_pix[0], targ_coords_pix[1], color='magenta', marker='+', s=50) + axes[0].text(targ_coords_pix[0], targ_coords_pix[1]+2, 'WCS', color='magenta', verticalalignment='bottom', horizontalalignment='center') + axes[0].text(0.95, 0.18, f'From WCS: {targ_coords_pix[0]:.2f}, {targ_coords_pix[1]:.2f}', + horizontalalignment='right', verticalalignment='bottom', + transform=axes[0].transAxes, + color='white') + + if verbose: + print(f"Star coords from WCS: {targ_coords_pix}") + if oss_cen_sci_pythonic is not None: + print(f"WCS offset = {np.asarray(targ_coords_pix) - oss_cen_sci_pythonic} pix") + except ImportError: oss_centroid_text = "" + ### WEBBPSF CENTROIDS ### cen = webbpsf.fwcentroid.fwcentroid(im_obs_clean) axes[0].scatter(cen[1], cen[0], color='red', marker='+', s=50) - axes[0].text(0.95, 0.05, f' Centroid: {cen[1]:.2f}, {cen[0]:.2f}'+oss_centroid_text, + axes[0].text(cen[1], cen[0], ' webbpsf', color='red', verticalalignment='center') + axes[0].text(0.95, 0.05, f' webbpsf Centroid: {cen[1]:.2f}, {cen[0]:.2f}'+oss_centroid_text, horizontalalignment='right', verticalalignment='bottom', transform=axes[0].transAxes, color='white') @@ -1737,6 +1772,6 @@ def nrc_ta_image_comparison(visitid, verbose=False, show_centroid=False): plt.tight_layout() - outname = f'wfs_ta_comparison_{visitid}.pdf' + outname = f'nrc_ta_comparison_{visitid}.pdf' plt.savefig(outname) print(f" => {outname}") From ff4904ebe6ec4b10524f115c01957ccb5ad76c7c Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Tue, 2 Apr 2024 12:31:12 -0400 Subject: [PATCH 4/4] improve match data handling of NRC TA images, again --- webbpsf/match_data.py | 9 +++++---- webbpsf/trending.py | 4 ++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/webbpsf/match_data.py b/webbpsf/match_data.py index eddba31a..60e19d5b 100644 --- a/webbpsf/match_data.py +++ b/webbpsf/match_data.py @@ -54,7 +54,9 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic if inst.name == 'NIRCam': if header['PUPIL'].startswith('MASK'): if header['PUPIL'] == 'MASKBAR': - inst.pupil_mask = header['CORONMSK'].replace('MASKA', 'MASK') + # the FITS header is just 'BAR' but the value needed in webbpsf is either + # 'MASKLWB' or "MASKSWB' depending on channel. + inst.pupil_mask = 'MASKLWB' if header['CHANNEL']=='LONG' else 'MASKSWB' else: inst.pupil_mask = header['PUPIL'] if 'CORONMSK' in header: @@ -65,9 +67,8 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic apername = get_nrc_coron_apname(header) inst.set_position_from_aperture_name(apername) - elif header['PUPIL'].startswith('F'): - inst.filter = header['PUPIL'] - else: + elif header['PUPIL'] != 'CLEAR' and not header['PUPIL'].startswith('F'): # no action needed for these + # note that filters in the pupil wheel were handled already above inst.pupil_mask = header['PUPIL'] elif inst.name == 'MIRI': diff --git a/webbpsf/trending.py b/webbpsf/trending.py index 1e8c17f6..870913a3 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1856,7 +1856,7 @@ def show_nrc_ta_img(visitid, ax=None, return_handles=False): return hdul, ax, norm, cmap, bglevel -def nrc_ta_image_comparison(visitid, verbose=False, show_centroid=False): +def nrc_ta_image_comparison(visitid, verbose=False, show_centroids=False): """ Retrieve a NIRCam target acq image and compare to a simulation Parameters: @@ -1902,7 +1902,7 @@ def nrc_ta_image_comparison(visitid, verbose=False, show_centroid=False): im_sim_scaled_aligned = im_sim_shifted*scalefactor # Plot - if show_centroid: + if show_centroids: ### OSS CENTROIDS ### # First, see if we can retrieve the on-board TA centroid measurment from the OSS engineering DB in MAST try: