diff --git a/webbpsf/mast_wss.py b/webbpsf/mast_wss.py index 4874aba7..da6a903a 100644 --- a/webbpsf/mast_wss.py +++ b/webbpsf/mast_wss.py @@ -589,7 +589,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 @@ -599,7 +599,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'] } @@ -616,10 +615,32 @@ 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 4c7aea8e..60e19d5b 100644 --- a/webbpsf/match_data.py +++ b/webbpsf/match_data.py @@ -54,19 +54,21 @@ 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'] - 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 # The apername keyword is not always correct for cases with dual-channel coronagraphy # in some such cases, APERNAME != PPS_APER. Let's ensure we have the proper apername for this channel: 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 de2e4c05..870913a3 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1827,8 +1827,8 @@ def delta_wfe_around_time(datetime, plot=True, ax=None, vmax=0.05, return_filena #### Functions for image comparisons -def show_wfs_ta_img(visitid, ax=None, return_handles=False): - """ Retrieve and display a WFS target acq image""" +def show_nrc_ta_img(visitid, ax=None, return_handles=False): + """ Retrieve and display a NIRCam target acq image""" hdul = webbpsf.mast_wss.get_visit_nrc_ta_image(visitid) @@ -1841,24 +1841,22 @@ 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) - if return_handles: return hdul, ax, norm, cmap, bglevel -def nrc_ta_image_comparison(visitid): + +def nrc_ta_image_comparison(visitid, verbose=False, show_centroids=False): """ Retrieve a NIRCam target acq image and compare to a simulation Parameters: @@ -1872,30 +1870,95 @@ 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 # Plot + 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: + 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']) + 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(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') + axes[1].imshow(im_sim_scaled_aligned, norm=norm, cmap=cmap, origin='lower') axes[1].set_title(f"Simulated PSF in F212N\nusing {opdname}") @@ -1917,6 +1980,6 @@ def nrc_ta_image_comparison(visitid): plt.tight_layout() - outname = f'wfs_ta_comparison_{visitid}.pdf' + outname = f'nrc_ta_comparison_{visitid}.pdf' plt.savefig(outname) print(f" => {outname}")