From fcdb538407a826287b5a2e5a78f753f6b55d2bd7 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Mon, 19 Feb 2024 10:15:58 -0500 Subject: [PATCH 1/5] add show_wfs_around_obs function --- webbpsf/mast_wss.py | 95 +++++++++++++++++++++++++++++++ webbpsf/trending.py | 132 ++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 221 insertions(+), 6 deletions(-) diff --git a/webbpsf/mast_wss.py b/webbpsf/mast_wss.py index 11e36358..6355f5c1 100644 --- a/webbpsf/mast_wss.py +++ b/webbpsf/mast_wss.py @@ -221,6 +221,7 @@ def get_opd_at_time(date, choice='closest', verbose=False, output_path=None): return mast_retrieve_opd(closest_fn, output_path = output_path) + ### Functions for format conversion of OPDs def import_wss_opd(filename, npix_out=1024, verbose=False): @@ -621,3 +622,97 @@ def set_params(parameters): ta_hdul = fits.open(mast_file_url) return ta_hdul + + +# Functions for retrieving metadata about science observations + +def _query_program_visit_times_by_inst(program, instrument, verbose=False): + """ Get the start and end times of all completed visits in a program, per instrument. + Not intended for general use; this is mostly a helper to query_program_visit_times. + + Getting the vststart_mjd and visitend_mjd fields requires using the instrument keywords + interface, so one has to specify which instrument ahead of time. + + Parameters + ---------- + program : int or str + Program ID + instrument : str + instrument name + verbose : bool + be more verbose in output? + + returns list of (visitid, start, end) tuples. + + """ + + from astroquery.mast import Mast + svc_table = {'MIRI':'Mast.Jwst.Filtered.Miri', + 'NIRCAM': 'Mast.Jwst.Filtered.NIRCam', + 'NIRSPEC': 'Mast.Jwst.Filtered.NIRSpec', + 'NIRISS': 'Mast.Jwst.Filtered.NIRISS', + } + + service = svc_table[instrument.upper()] + + collist = 'filename, program, observtn, visit_id, vststart_mjd, visitend_mjd, bstrtime' + all_columns = False + + def set_params(parameters): + return [{"paramName" : p, "values" : v} for p, v in parameters.items()] + + + keywords = {'program': [str(program),]} + parameters = {'columns': '*' if all_columns else collist, + 'filters': set_params(keywords)} + + if verbose: + print("MAST query parameters:") + print(parameters) + + responsetable = Mast.service_request(service, parameters) + responsetable.sort(keys='bstrtime') + + + visit_times = [] + + for row in responsetable: + visit_times.append( ('V'+row['visit_id'], row['vststart_mjd'], row['visitend_mjd'])) + + visit_times= set(visit_times) + return list(visit_times) + +def query_program_visit_times(program, verbose=False): + """ Get the start and end times of all completed visits in a program. + + Parameters + ---------- + program : int or str + Program ID + verbose : bool + be more verbose in output? + + Returns astropy Table with columns for visit ID and start and end times. + """ + + from astroquery.mast import Observations + obs = Observations.query_criteria(obs_collection=["JWST"], proposal_id=[program]) + # Annoyingly, that query interface doesn't return start/end times + instruments = [val.split('/')[0] for val in set(obs['instrument_name'])] + + visit_times = [] + for inst in instruments: + if verbose: + print(f"querying for visits using {inst}") + visit_times += _query_program_visit_times_by_inst(program, inst) + + + vids = [v[0] for v in visit_times] + starts =astropy.time.Time([float(v[1]) for v in visit_times], format='mjd') + ends = astropy.time.Time([float(v[2]) for v in visit_times], format='mjd') + + #visit_times = np.asarray(visit_times) + return astropy.table.Table([vids, starts, ends], + names=('visit_id', 'start_mjd', 'end_mjd')) + + diff --git a/webbpsf/trending.py b/webbpsf/trending.py index 0815427c..38a7143a 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -952,12 +952,11 @@ def get_month_start_end(year, month): return start_date, end_date -def filter_opdtable_for_month(year, mon, opdtable): - """Filter existing opdtable for a given month - This includes the last measurement in the prior month too (if applicable), so we can compute a delta +def filter_opdtable_for_daterange(start_date, end_date, opdtable): + """Filter existing opdtable for a given time range + This includes the last measurement in the prior time range too (if applicable), so we can compute a delta to the first one """ - start_date, end_date = get_month_start_end(year, mon) # Start a little early, such that we are going to have at least 1 WFS before the start date pre_start_date = astropy.time.Time(start_date) - astropy.time.TimeDelta(4 * u.day) opdtable = webbpsf.mast_wss.filter_opd_table( @@ -975,6 +974,30 @@ def filter_opdtable_for_month(year, mon, opdtable): return opdtable +def filter_opdtable_for_month(year, mon, opdtable): + """Filter existing opdtable for a given month + This includes the last measurement in the prior month too (if applicable), so we can compute a delta + to the first one + """ + start_date, end_date = get_month_start_end(year, mon) + return filter_opdtable_for_daterange(start_date, end_date, opdtable) + + + +def get_opdtable_for_daterange(start_date, end_date): + """Return table of OPD measurements for date range. + + This includes the last measurement preceding this date range, too, so we + can compute the first delta at the start of this range. + """ + # Retrieve full OPD table, then trim to the selected time period + opdtable0 = webbpsf.mast_wss.retrieve_mast_opd_table() + opdtable0 = webbpsf.mast_wss.deduplicate_opd_table(opdtable0) + + opdtable = filter_opdtable_for_daterange(start_date, end_date, opdtable0) + return opdtable + + def get_opdtable_for_month(year, mon): """Return table of OPD measurements for a given month. retrieve the opd table and filter according to month/year designation @@ -1043,8 +1066,6 @@ def get_dates_for_pid(pid, project='jwst'): return - - def monthly_trending_plot(year, month, verbose=True, instrument='NIRCam', filter='F200W', vmax=200, pid=None, opdtable=None): """Make monthly trending plot showing OPDs, mirror moves, RMS WFE, and the resulting PSF EEs @@ -1416,6 +1437,7 @@ def plot_phase_retrieval_crosscheck(fn, vmax_fraction=1.0): return fig + def plot_wfs_obs_delta(fn1, fn2, vmax_fraction=1.0): """ Display comparison of two weak lens observations @@ -1516,6 +1538,8 @@ def show_wfs_around_obs(filename, verbose='True'): """Make a helpful plot showing available WFS before and after some given science observation. This can be used to help inform how much WFE variability there was around that time. + See also show_wfs_during_program. + Parameters ---------- filename : str @@ -1624,6 +1648,102 @@ def vprint(*args, **kwargs): ax4.set_title("Delta WFE\nAfter-Before", color='C1', fontweight='bold') +def show_wfs_during_program(program, verbose=False, ax = None, ref_wavefront_date=None): + """ Show WFS data for the entire time interval in which a given program was observed. + + Plots time series of the WFS measured RMS WFE as seen in NIRCam, the start and end times + of all the observations in that program, and the delta wavefront RMS relative to either + the median wavefront during that time period, or to a specifed date. + + See also show_wfs_around_obs. + + Parameters + ---------- + program : str or int + Program ID number + verbose : str + be more verbose in output? + ax : matplotlib.Axes instance, or None + an existing Axes to plot into, or else a new one will be created. + ref_wavefront_date : date-like str or None + Optional, to specify which date's wavefront sensing should be used as + the reference wavefront for computation of the plotted delta WFE RMS. + The closest wavefront in time to the specified date will be used. + If not set, the median wavefront over the entire time period will be used. + """ + + # Query mast for when the observations took place + science_visit_table = webbpsf.mast_wss.query_program_visit_times(program, verbose=verbose) + science_visit_table.sort(keys=['start_mjd']) + + # Figure out reasonable start and end dates for the time interval to display + sci_duration = science_visit_table['start_mjd'].max() - science_visit_table['start_mjd'].min() + plot_padding_time_range = max(sci_duration.to(u.day)*0.2, 4*u.day) + start_date = science_visit_table['start_mjd'].min() - plot_padding_time_range + end_date = science_visit_table['end_mjd'].max() + plot_padding_time_range + + # Look up wavefront sensing and mirror move corrections for that range + opdtable = get_opdtable_for_daterange(start_date, end_date) + corrections_table = webbpsf.mast_wss.get_corrections(opdtable) + + # Iterate over the WFS measurements to retrieve the OPDs and RMS WFE + wfs_dates = [] + rms_obs = [] + opds = [] + for row_index in range(len(opdtable)): + opd_fn = opdtable[row_index]['fileName'] + opd, opd_hdul = webbpsf.trending._read_opd(opd_fn) + + opds.append(opd) + wfs_dates.append(opdtable[row_index]['date']) + rms_obs.append(opd_hdul[1].header['RMS_WFE']*1000) + wfs_dates_array = astropy.time.Time(wfs_dates, format='isot') + opd_array = np.asarray(opds) + + # Compute delta WFE, either relative to median or a specifed date + if ref_wavefront_date: + refdate = astropy.time.Time(ref_wavefront_date) + delta_times = wfs_dates_array - refdate + closest = np.argmin(np.abs(delta_times)) + reference_opd = opds[closest] + ref_label = "WFS near "+ref_wavefront_date + if verbose: + print(f"for date {ref_wavefront_date}, using opd {closest}") + else: + median_opd = np.median(opd_array, axis=0) + reference_opd = median_opd + ref_label = 'median wavefront' + delta_opds = opd_array - reference_opd + mask = opd_array.sum(axis=0)!=0 + delta_rmses = [webbpsf.utils.rms(d, mask=mask)*1000 for d in delta_opds] + + # Plot! + if ax is None: + fig, ax = plt.subplots(figsize=(8,4), ncols=1, nrows=1) + + ax.plot_date(wfs_dates_array.plot_date, rms_obs, '+', + color='C1', ls='-', label='Measured RMS Wavefront Error at NIRCam NRCA3') + ax.plot_date(wfs_dates_array.plot_date, delta_rmses,'none', + color='C0', ls='--', label=f'RMS Delta Wavefront relative to {ref_label}') + + plot_sci_y = 25 + ax.scatter(science_visit_table['start_mjd'].plot_date, + [plot_sci_y]*len(science_visit_table), s=150, marker="*", color='black' , + label=f'Program {program} observations') + for row in science_visit_table: + ax.fill_betweenx([0,120], row['start_mjd'].plot_date, (row['end_mjd']).plot_date, + color='gray', alpha=0.2) + ax.text(row['start_mjd'].plot_date, plot_sci_y+4, row['visit_id'], rotation=45, fontsize='small') + + ax.set_ylim(0,120) + ax.legend(framealpha=0.99) + ax.xaxis.set_major_formatter( + matplotlib.dates.ConciseDateFormatter(ax.xaxis.get_major_locator())) + ax.set_ylabel("Measured Wavefront Error RMS [nm]") + ax.set_xlim(start_date.plot_date, end_date.plot_date) + + + #### Functions for image comparisons def show_wfs_ta_img(visitid, ax=None, return_handles=False): """ Retrieve and display a WFS target acq image""" From a63a4c5be6dbf6c08ecf813e24c6c58e1ddeab54 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Thu, 21 Mar 2024 08:59:21 -0400 Subject: [PATCH 2/5] fix edge case bug that might lead to unnecessary repeat queries of the same instrument --- webbpsf/mast_wss.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/webbpsf/mast_wss.py b/webbpsf/mast_wss.py index 6355f5c1..4874aba7 100644 --- a/webbpsf/mast_wss.py +++ b/webbpsf/mast_wss.py @@ -698,7 +698,7 @@ def query_program_visit_times(program, verbose=False): from astroquery.mast import Observations obs = Observations.query_criteria(obs_collection=["JWST"], proposal_id=[program]) # Annoyingly, that query interface doesn't return start/end times - instruments = [val.split('/')[0] for val in set(obs['instrument_name'])] + instruments = set([val.split('/')[0] for val in set(obs['instrument_name'])]) visit_times = [] for inst in instruments: From a2918d06ab48207916cdd9caec52ea5629be7d12 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Thu, 21 Mar 2024 08:59:58 -0400 Subject: [PATCH 3/5] improve options and flexibility for show_wfs_during_program --- webbpsf/trending.py | 45 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/webbpsf/trending.py b/webbpsf/trending.py index 38a7143a..c7d1cbe4 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1648,7 +1648,9 @@ def vprint(*args, **kwargs): ax4.set_title("Delta WFE\nAfter-Before", color='C1', fontweight='bold') -def show_wfs_during_program(program, verbose=False, ax = None, ref_wavefront_date=None): +def show_wfs_during_program(program, verbose=False, ax = None, ref_wavefront_date=None, + ref_wavefront_visit=None, + start_date=None, end_date=None): """ Show WFS data for the entire time interval in which a given program was observed. Plots time series of the WFS measured RMS WFE as seen in NIRCam, the start and end times @@ -1666,21 +1668,40 @@ def show_wfs_during_program(program, verbose=False, ax = None, ref_wavefront_dat ax : matplotlib.Axes instance, or None an existing Axes to plot into, or else a new one will be created. ref_wavefront_date : date-like str or None - Optional, to specify which date's wavefront sensing should be used as - the reference wavefront for computation of the plotted delta WFE RMS. - The closest wavefront in time to the specified date will be used. - If not set, the median wavefront over the entire time period will be used. + Optional, to specify which date's wavefront sensing should be used as + the reference wavefront for computation of the plotted delta WFE RMS. + The closest wavefront in time to the specified date will be used. + If not set, the median wavefront over the entire time period will be used. + ref_wavefront_visit : str + Another way to specify which date's wavefront sensing should be used as the + reference. Provide a science visit ID from the program (e.g. "V01234002001") + and the closest wavefront to the observation start time of that visit will be used. + start_date, end_date : strings or astropy.time.Time + Start and/or end dates for the time period to display. If not set, + reasonable default values will be computed that include all observations + for that science program plus some padding time on either side. """ # Query mast for when the observations took place science_visit_table = webbpsf.mast_wss.query_program_visit_times(program, verbose=verbose) science_visit_table.sort(keys=['start_mjd']) + if verbose: + for row in science_visit_table: + print(f" Found {row['visit_id']} starting at {row['start_mjd'].iso}") + - # Figure out reasonable start and end dates for the time interval to display + # Figure out reasonable start and end dates for the time interval to display, + # or use values provided by the user sci_duration = science_visit_table['start_mjd'].max() - science_visit_table['start_mjd'].min() plot_padding_time_range = max(sci_duration.to(u.day)*0.2, 4*u.day) - start_date = science_visit_table['start_mjd'].min() - plot_padding_time_range - end_date = science_visit_table['end_mjd'].max() + plot_padding_time_range + if start_date is None: + start_date = science_visit_table['start_mjd'].min() - plot_padding_time_range + else: + start_date = astropy.time.Time(start_date) + if end_date is None: + end_date = science_visit_table['end_mjd'].max() + plot_padding_time_range + else: + end_date = astropy.time.Time(end_date) # Look up wavefront sensing and mirror move corrections for that range opdtable = get_opdtable_for_daterange(start_date, end_date) @@ -1701,6 +1722,14 @@ def show_wfs_during_program(program, verbose=False, ax = None, ref_wavefront_dat opd_array = np.asarray(opds) # Compute delta WFE, either relative to median or a specifed date + if ref_wavefront_visit is not None: + for row in science_visit_table: + if row['visit_id']==ref_wavefront_visit: + ref_wavefront_date = row['start_mjd'].iso + if verbose: + print(f" Will compare wavefronts relative to {row['visit_id']} starting at {row['start_mjd'].iso}") + + if ref_wavefront_date: refdate = astropy.time.Time(ref_wavefront_date) delta_times = wfs_dates_array - refdate From 59b0d0855a12749d0271d4c9277ac4368b6fdce6 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Thu, 21 Mar 2024 10:01:59 -0400 Subject: [PATCH 4/5] improve verbose message outputs --- webbpsf/trending.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/webbpsf/trending.py b/webbpsf/trending.py index c7d1cbe4..510b3a0b 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1737,11 +1737,14 @@ def show_wfs_during_program(program, verbose=False, ax = None, ref_wavefront_dat reference_opd = opds[closest] ref_label = "WFS near "+ref_wavefront_date if verbose: - print(f"for date {ref_wavefront_date}, using opd {closest}") + print(f"Computing delta OPDs relative to date {ref_wavefront_date}, using opd {opdtable[closest]['fileName']}") else: median_opd = np.median(opd_array, axis=0) reference_opd = median_opd ref_label = 'median wavefront' + if verbose: + print("Computing delta OPDs relative to median OPD over that time period.") + delta_opds = opd_array - reference_opd mask = opd_array.sum(axis=0)!=0 delta_rmses = [webbpsf.utils.rms(d, mask=mask)*1000 for d in delta_opds] From 17bd67601adbd93e536260b7a7de1a95c8803a55 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Mon, 1 Apr 2024 14:17:45 -0400 Subject: [PATCH 5/5] add auto download of OPD files in show_wfs_during_program if they were not previously retrieved --- webbpsf/trending.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/webbpsf/trending.py b/webbpsf/trending.py index 510b3a0b..9af92b60 100644 --- a/webbpsf/trending.py +++ b/webbpsf/trending.py @@ -1713,7 +1713,11 @@ def show_wfs_during_program(program, verbose=False, ax = None, ref_wavefront_dat opds = [] for row_index in range(len(opdtable)): opd_fn = opdtable[row_index]['fileName'] - opd, opd_hdul = webbpsf.trending._read_opd(opd_fn) + try: + opd, opd_hdul = webbpsf.trending._read_opd(opd_fn) + except FileNotFoundError: + webbpsf.mast_wss.mast_retrieve_opd(opd_fn, verbose=verbose) + opd, opd_hdul = webbpsf.trending._read_opd(opd_fn) opds.append(opd) wfs_dates.append(opdtable[row_index]['date'])