From 2792c6c6d84577a6a320628bb06b085477d387ac Mon Sep 17 00:00:00 2001 From: anand_raichoor Date: Wed, 7 Aug 2024 11:10:34 -0700 Subject: [PATCH 1/7] add create_ctedet_rowbyrow_pdf() and few related utilities --- py/desispec/night_qa.py | 238 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 226 insertions(+), 12 deletions(-) diff --git a/py/desispec/night_qa.py b/py/desispec/night_qa.py index 9f29c5de4..763178c1e 100644 --- a/py/desispec/night_qa.py +++ b/py/desispec/night_qa.py @@ -26,10 +26,11 @@ from desitarget.geomask import match_to # AR desispec from desispec.fiberbitmasking import get_skysub_fiberbitmask_val -from desispec.io import findfile +from desispec.io import findfile, read_image from desispec.io.util import replace_prefix, get_tempfilename from desispec.calibfinder import CalibFinder from desispec.scripts import preproc +from desispec.correct_cte import get_rowbyrow_image_model from desispec.tile_qa_plot import get_tilecov # AR matplotlib import matplotlib @@ -45,6 +46,9 @@ cameras = ["b", "r", "z"] petals = np.arange(10, dtype=int) +ctedet_mydicts = None + + def get_nightqa_outfns(outdir, night): """ Utility function to get nightqa file names. @@ -819,9 +823,9 @@ def create_badcol_png(outpng, night, prod, n_previous_nights=10): os.rename(tmp_outpng, outpng) -def _read_ctedet(night, prod, ctedet_expid, petal, camera): +def _read_ctedet_campet(night, prod, ctedet_expid, petal, camera): """ - Internal function used by create_ctedet_pdf(), reading the ctedet_expid preproc info. + Internal function used by _read_ctedet reading the ctedet_expid preproc info for one {camera, petal}. Args: night: night (int) @@ -839,9 +843,10 @@ def _read_ctedet(night, prod, ctedet_expid, petal, camera): if os.path.isfile(fn): mydict = {} mydict["fn"] = fn - # AR read - with fitsio.FITS(fn) as fx: - mydict["img"] = fx["IMAGE"].read() + # AR read (use read_image(), so that it can be used by get_rowbyrow_image_model) + #with fitsio.FITS(fn) as fx: + # mydict["img"] = fx["IMAGE"].read() + mydict["img"] = read_image(fn) # AR check if we re displaying a 1s FLAT hdr = fitsio.read_header(fn, "IMAGE") if (hdr["OBSTYPE"] == "FLAT") & (hdr["REQTIME"] == 1): @@ -860,6 +865,41 @@ def _read_ctedet(night, prod, ctedet_expid, petal, camera): return None +def _read_ctedet(night, prod, ctedet_expid, nproc): + """ + Internal function used by create_ctedet_pdf() and create_ctedet_rowbyrow_pdf(), + reading the ctedet_expid preproc info. + + Args: + night: night (int) + prod: full path to prod folder, e.g. /global/cfs/cdirs/desi/spectro/redux/blanc (string) + ctedet_expid: EXPID for the CTE diagnosis (1s FLAT, or darker science exposure) (int) + nproc: number of processes running at the same time (int) + + Returns: + mydicts: a list of dictionaries with the IMAGE data, plus various infos, + assuming a preproc file is available, otherwise ``None``. + the list is for [b0, r0, z0, b1, r1, z1, ..., b9, r9, z9] + """ + myargs = [] + for petal in petals: + for camera in cameras: + myargs.append( + [ + night, + prod, + ctedet_expid, + petal, + camera, + ] + ) + # AR launching pool + pool = multiprocessing.Pool(processes=nproc) + with pool: + mydicts = pool.starmap(_read_ctedet_campet, myargs) + return mydicts + + def create_ctedet_pdf(outpdf, night, prod, ctedet_expid, nproc, nrow=21, xmin=None, xmax=None, ylim=(-5, 10), yoffcols={"A" : -2.5, "B" : -3.5, "C" : -2.5, "D" : -3.5}, ): @@ -887,6 +927,12 @@ def create_ctedet_pdf(outpdf, night, prod, ctedet_expid, nproc, nrow=21, xmin=No # AR temporary output file tmp_outpdf = get_tempfilename(outpdf) + # AR read ctedet + global ctedet_mydicts + if ctedet_mydicts is None: + ctedet_mydicts = _read_ctedet(night, prod, ctedet_expid, nproc) + + # AR myargs = [] for petal in petals: for camera in cameras: @@ -899,16 +945,13 @@ def create_ctedet_pdf(outpdf, night, prod, ctedet_expid, nproc, nrow=21, xmin=No camera, ] ) - # AR launching pool - pool = multiprocessing.Pool(processes=nproc) - with pool: - mydicts = pool.starmap(_read_ctedet, myargs) + # AR plotting clim = (-5, 5) tmpcols = plt.rcParams["axes.prop_cycle"].by_key()["color"] colors = {"A" : tmpcols[1], "B" : tmpcols[1], "C" : tmpcols[0], "D" : tmpcols[0]} with PdfPages(tmp_outpdf) as pdf: - for mydict in mydicts: + for mydict in ctedet_mydicts: petcam_xmin, petcam_xmax = xmin, xmax fig = plt.figure(figsize=(30, 5)) gs = gridspec.GridSpec(2, 1, wspace=0.1, height_ratios = [1, 4]) @@ -928,7 +971,7 @@ def create_ctedet_pdf(outpdf, night, prod, ctedet_expid, nproc, nrow=21, xmin=No color="k", fontsize=20, fontweight="bold", ha="center", va="center", transform=ax1d.transAxes, ) - img = mydict["img"] + img = mydict["img"].pix ny, nx = img.shape if petcam_xmin is None: petcam_xmin = 0 @@ -972,6 +1015,177 @@ def create_ctedet_pdf(outpdf, night, prod, ctedet_expid, nproc, nrow=21, xmin=No os.rename(tmp_outpdf, outpdf) +def _compute_rowbyrow(ims, nproc): + """ + Internal function used by create_ctedet_rowbyrow_pdf(), + executing get_rowbyrow_image_model() on a list of read_image(ctedet_fn). + + Args: + ims: list of im objects, outputs of read_image(ctedet_fn); can contain some None values (list of objects) + nproc: number of processes running at the same time (int) + + Returns: + mods: list of mod objects, outputs of get_rowbyrow_image_model(im) for each im of ims (list of objects) + + Note: + mods[i]=None if ims[i]=None + """ + # AR indexes with ims != None + ii = [i for i in range(len(ims)) if ims[i] is not None] + + if len(ii) == 0: + return [None for _ in range(len(ims))] + + else: + tmp_ims = [ims[i] for i in ii] + pool = multiprocessing.Pool(nproc) + with pool: + tmp_mods = pool.map(get_rowbyrow_image_model, tmp_ims) + + # AR mods with mods[i]=None if ims[i]=None + mods = [None for _ in range(len(ims))] + for i, tmp_mod in zip(ii, tmp_mods): + mods[i] = tmp_mod + # + return mods + + +def create_ctedet_rowbyrow_pdf(outpdf, night, prod, ctedet_expid, nproc, dpi=100): + """ + For a given night, create a pdf with a CTE diagnosis (from preproc files). + + Args: + outpdf: output pdf file (string) + night: night (int) + prod: full path to prod folder, e.g. /global/cfs/cdirs/desi/spectro/redux/blanc (string) + ctedet_expid: EXPID for the CTE diagnosis (1s FLAT, or darker science exposure) (int) + nproc: number of processes running at the same time (int) + dpi (optional, default to 100): dpi value in plt.savefig() (int) + + Note: + * Credits to E. Schlafly. + """ + + + # AR temporary output file + tmp_outpdf = get_tempfilename(outpdf) + + campets = [] + for petal in petals: + for camera in cameras: + campets.append("{}{}".format(camera, petal)) + + # AR read ctedet images + global ctedet_mydicts + if ctedet_mydicts is None: + ctedet_mydicts = _read_ctedet(night, prod, ctedet_expid, nproc) + + # AR compute the row-by-row model + ims = [ctedet_mydict["img"] for ctedet_mydict in ctedet_mydicts] + mods = _compute_rowbyrow(ims, nproc) + + # AR plotting + # AR remarks: + # AR - the (x,y) conversions for the side panels + # AR are a bit counter-intuitive, as ax.imshow() + # AR reverses the displayed axes... + # AR - the panels positioning is not very elegant, + # AR probably could be coded in a nicer way! + clim = (-2, 2) + cmap = matplotlib.cm.coolwarm + tmpcols = plt.rcParams["axes.prop_cycle"].by_key()["color"] + colors = {"A" : tmpcols[1], "B" : tmpcols[1], "C" : tmpcols[0], "D" : tmpcols[0]} + with PdfPages(tmp_outpdf) as pdf: + + for campet, mydict, mod in zip(campets, ctedet_mydicts, mods): + camera, petal = campet[0], campet[1] + title = "EXPID={} {}{}".format(ctedet_expid, camera, petal) + + if camera == cameras[0]: + fig = plt.figure(figsize=(20 * len(cameras) / 3, 10)) + gs = gridspec.GridSpec(1, len(cameras), wspace=0.1) + ic = 0 + + ax = fig.add_subplot(gs[0, ic]) + ax.set_title(title) + ax.set_xlabel(r"Fiber direction $\Longrightarrow$") + if ic == 0: + ax.set_ylabel(r"Wavelength direction $\Longrightarrow$") + if mydict is not None: + assert(mydict["img"].meta["CAMERA"].lower() == campet.lower()) + ress = mydict["img"].pix - mod + median = np.nanmedian(ress) + ax.text(0.05, 0.95, "res_median = {:.2f}".format(median), transform=ax.transAxes) + im = ax.imshow(ress - median, cmap=cmap, vmin=clim[0], vmax=clim[1]) + # AR amp labels + cte correction infos + # AR ampsec: + # AR - in the header they are 1-index + # AR - we change those to 0-index (consistent with OFFCOLS and CTECOLS) + tmpdy = -0.035 * ress.shape[0] + for amp in ["a", "b", "c", "d"]: + key = "ampsec{}".format(amp) + # 2-amp readout only has 2 amps, so check + if key in mydict["img"].meta: + ampsec = mydict["img"].meta[key] + ampsecx, ampsecy = ampsec.replace("[", "").replace("]", "").split(",") + ampsecx = "{}:{}".format( + int(ampsecx.split(":")[0]) - 1, + int(ampsecx.split(":")[1]) - 1, + ) + ampsecy = "{}:{}".format( + int(ampsecy.split(":")[0]) - 1, + int(ampsecy.split(":")[1]) - 1, + ) + ampsec_old = ampsec + ampsec = "[{}, {}]".format(ampsecx, ampsecy) + tmpx = np.mean([int(_) for _ in ampsecx.split(":")]) + tmpy = np.mean([int(_) for _ in ampsecy.split(":")]) + ax.text(tmpx, tmpy, amp.upper(), fontweight="bold", ha="center", va="center") + tmpy += tmpdy + ax.text(tmpx, tmpy, "[{}, {}]".format(ampsecx, ampsecy), fontweight="bold", ha="center", va="center") + tmpy += tmpdy + #ax.text(ampsecx, ampsecy, "{}\n{}".format(amp.upper(), ampsec), fontweight="bold", ha="center", va="center") + # AR known columns with offset? + # AR stored in format like "12:24,1700:1900" + colmins, colmaxs, colors = [], [], [] + for key, color in zip( + ["OFFCOLS{}".format(amp.upper()) , "CTECOLS{}".format(amp.upper())], + ["g", "darkgreen"], + ): + if key in mydict and mydict[key] is not None: + ax.text(tmpx, tmpy, "{}: {}".format(key, mydict[key]), color=color, fontweight="bold", ha="center", va="center") + tmpy += tmpdy + for colrange in mydict[key].split(","): + colmins.append(int(colrange.split(":")[0])) + colmaxs.append(int(colrange.split(":")[1])) + colors.append(color) + for colmin, colmax, color in zip(colmins, colmaxs, colors): + ax.plot([colmin, colmax], [tmpy, tmpy], color=color, lw=3, alpha=0.75, zorder=1) + tmpy += tmpdy + # AR flip the y-axis to have y coords increasing towards up + # AR (e.g. see Guy+2023 Fig.4) + ax.set_ylim(ax.get_ylim()[::-1]) + if (camera == cameras[-1]): + p = ax.get_position().get_points().flatten() + cax = fig.add_axes([ + p[0] + 1.05 * (p[2] - p[0]), + p[1], + 0.05 * (p[2] - p[0]), + 1.0 * (p[3] - p[1]) + ]) + cbar = plt.colorbar(im, ax=ax, cax=cax, orientation="vertical", ticklocation="right", pad=0, extend="both") + clab = "preproc - rowbyrow_model - res_median [Units : ?]".format(median) + cbar.set_label(clab) + cbar.mappable.set_clim(clim) + ic += 1 + if camera == cameras[-1]: + pdf.savefig(fig, bbox_inches="tight", dpi=dpi) + plt.close() + + # AR move to final location + os.rename(tmp_outpdf, outpdf) + + def _read_sframesky(night, prod, expid): """ Internal function called by create_sframesky_pdf, which generates From 6395fe277fcf6680ad597bbea114be5f97bec556 Mon Sep 17 00:00:00 2001 From: anand_raichoor Date: Wed, 7 Aug 2024 11:13:23 -0700 Subject: [PATCH 2/7] add ctedetrowbyrow in get_nightqa_outfns() --- py/desispec/night_qa.py | 1 + 1 file changed, 1 insertion(+) diff --git a/py/desispec/night_qa.py b/py/desispec/night_qa.py index 763178c1e..eef4719cb 100644 --- a/py/desispec/night_qa.py +++ b/py/desispec/night_qa.py @@ -66,6 +66,7 @@ def get_nightqa_outfns(outdir, night): "morningdark" : os.path.join(outdir, "morningdark-{}.pdf".format(night)), "badcol" : os.path.join(outdir, "badcol-{}.png".format(night)), "ctedet" : os.path.join(outdir, "ctedet-{}.pdf".format(night)), + "ctedetrowbyrow" : os.path.join(outdir, "ctedetrowbyrow-{}.pdf".format(night)), "sframesky" : os.path.join(outdir, "sframesky-{}.pdf".format(night)), "tileqa" : os.path.join(outdir, "tileqa-{}.pdf".format(night)), "skyzfiber" : os.path.join(outdir, "skyzfiber-{}.png".format(night)), From 7be9cdd57da0c1f1acc967443d2f4fa9dbc0ac61 Mon Sep 17 00:00:00 2001 From: anand_raichoor Date: Wed, 7 Aug 2024 11:44:29 -0700 Subject: [PATCH 3/7] update write_nightqa_html() for ctedetrowbyrow --- py/desispec/night_qa.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/py/desispec/night_qa.py b/py/desispec/night_qa.py index eef4719cb..91bc0f745 100644 --- a/py/desispec/night_qa.py +++ b/py/desispec/night_qa.py @@ -2315,19 +2315,21 @@ def write_nightqa_html(outfns, night, prod, css, expids, tileids, surveys): # AR - morningdark # AR - badcol # AR - ctedet + # AR - ctedetrowbyrow # AR - sframesky # AR - tileqa # AR - skyzfiber # AR - petalnz for case, caselab, width, text in zip( - ["dark", "morningdark", "badcol", "ctedet", "sframesky", "tileqa", "skyzfiber", "petalnz"], - ["DARK", "Morning DARK", "bad columns", "CTE detector", "sframesky", "Tile QA", "SKY Z vs. FIBER", "Per-petal n(z)"], - ["100%", "100%", "35%", "100%", "75%", "90%", "90%", "100%"], + ["dark", "morningdark", "badcol", "ctedet", "ctedetrowbyrow", "sframesky", "tileqa", "skyzfiber", "petalnz"], + ["DARK", "Morning DARK", "bad columns", "CTE detector (amps boundaries)", "CTE detector (row-by-row)", "sframesky", "Tile QA", "SKY Z vs. FIBER", "Per-petal n(z)"], + ["100%", "100%", "35%", "100%", "100%", "75%", "90%", "90%", "100%"], [ "This pdf displays the 300s (binned) DARK (one page per spectrograph; non-valid pixels are displayed in red)\nThe side panels report the median profiles for each pair of amps along each direction.\nWatch it and report unsual features (easy to say!)", "This pdf displays the 1200s (binned) morning DARK (one page per spectrograph; non-valid pixels are displayed in red)\nThe side panels report the median profiles for each pair of amps along each direction.\nWatch it and report unsual features (easy to say!)", "This plot displays the histograms of the bad columns.\nWatch it and report unsual features (easy to say!)", "This pdf displays a small diagnosis to detect CTE anormal behaviour (one petal-camera per page)\nWatch it and report unusual features (typically if the lower enveloppe of the blue or orange curve is systematically lower than the other one).", + "This pdf displays another diagnosis to detect CTE anormal behaviour (one petal per page)\nWatch it and report unusual features (typically if some significant 'block pattern' along some fibers is not already identified by a OFFCOLS or CTECOLS correction).", "This pdf displays the sframe image for the sky fibers for each Main exposure (one exposure per page).\nPixels with IVAR=0 are displayed in yellow.\nWatch it and report unsual features (easy to say!)", "This pdf displays the tile-qa-TILEID-thru{}.png files for the Main tiles (one tile per page).\nWatch it, in particular the Z vs. FIBER plot, and report unsual features (easy to say!)".format(night), "This plot displays all the SKY fibers for the {} night.\nWatch it and report unsual features (easy to say!)".format(night), From c980704664db6c84a54ba5038212e25295013ada Mon Sep 17 00:00:00 2001 From: anand_raichoor Date: Wed, 7 Aug 2024 11:45:39 -0700 Subject: [PATCH 4/7] updates for newly introduced ctedetrowbyrow --- bin/desi_night_qa | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/bin/desi_night_qa b/bin/desi_night_qa index c66b50eed..ea000da68 100755 --- a/bin/desi_night_qa +++ b/bin/desi_night_qa @@ -28,6 +28,7 @@ from desispec.night_qa import ( create_dark_pdf, create_badcol_png, create_ctedet_pdf, + create_ctedet_rowbyrow_pdf, create_sframesky_pdf, create_tileqa_pdf, create_skyzfiber_png, @@ -58,7 +59,7 @@ def parse(options=None): parser.add_argument("--steps", type = str, default = ",".join(steps_all), required = False, help = "comma-separated list of steps to execute (default={})".format(",".join(steps_all))) parser.add_argument("--nproc", type = int, default = 1, required = False, - help="number of parallel processes for create_dark_pdf(), create_ctedet_pdf() and create_sframesky_pdf() (default=1)") + help="number of parallel processes for create_dark_pdf(), create_ctedet_pdf(), create_ctedet_rowbyrow_pdf() and create_sframesky_pdf() (default=1)") parser.add_argument("--dark-bkgsub-science-cameras", type = str, default = "b", required = False, help="for the dark/morningdark, comma-separated list of the cameras to be additionally processed with the --bkgsub-for-science argument (default=b)") @@ -169,6 +170,7 @@ def main(): if "ctedet" in steps_tbd: if ctedet_expid is not None: create_ctedet_pdf(outfns["ctedet"], args.night, args.prod, ctedet_expid, args.nproc) + create_ctedet_rowbyrow_pdf(outfns["ctedetrowbyrow"], args.night, args.prod, ctedet_expid, args.nproc) # AR sframesky if "sframesky" in steps_tbd: From 094c3be2f5c788d72ae93efd3b83442912f121e7 Mon Sep 17 00:00:00 2001 From: anand_raichoor Date: Wed, 7 Aug 2024 13:53:16 -0700 Subject: [PATCH 5/7] fix in _read_ctedet() docstr --- py/desispec/night_qa.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/py/desispec/night_qa.py b/py/desispec/night_qa.py index 91bc0f745..88f821516 100644 --- a/py/desispec/night_qa.py +++ b/py/desispec/night_qa.py @@ -879,8 +879,8 @@ def _read_ctedet(night, prod, ctedet_expid, nproc): Returns: mydicts: a list of dictionaries with the IMAGE data, plus various infos, - assuming a preproc file is available, otherwise ``None``. - the list is for [b0, r0, z0, b1, r1, z1, ..., b9, r9, z9] + assuming a preproc file is available, otherwise ``None``. + the list is for [b0, r0, z0, b1, r1, z1, ..., b9, r9, z9] """ myargs = [] for petal in petals: From fbdba5686063d268f6a684364bb0d7626f7fb674 Mon Sep 17 00:00:00 2001 From: anand_raichoor Date: Wed, 7 Aug 2024 15:19:35 -0700 Subject: [PATCH 6/7] bugfix: handling of ctedetrowbyrow --- bin/desi_night_qa | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bin/desi_night_qa b/bin/desi_night_qa index ea000da68..35c4d0f68 100755 --- a/bin/desi_night_qa +++ b/bin/desi_night_qa @@ -148,7 +148,7 @@ def main(): morningdark_expid = get_morning_dark_night_expid(args.night, args.prod) # AR CTE detector expid - if "ctedet" in steps_tbd: + if np.in1d(["ctedet", "ctedetrowbyrow"], steps_tbd).sum() > 0: ctedet_expid = get_ctedet_night_expid(args.night, args.prod) # AR dark @@ -167,7 +167,7 @@ def main(): create_badcol_png(outfns["badcol"], args.night, args.prod) # AR CTE detector - if "ctedet" in steps_tbd: + if np.in1d(["ctedet", "ctedetrowbyrow"], steps_tbd).sum() > 0: if ctedet_expid is not None: create_ctedet_pdf(outfns["ctedet"], args.night, args.prod, ctedet_expid, args.nproc) create_ctedet_rowbyrow_pdf(outfns["ctedetrowbyrow"], args.night, args.prod, ctedet_expid, args.nproc) From 6af0fd676d5d84bc300a7ac02c427316ccda241e Mon Sep 17 00:00:00 2001 From: anand_raichoor Date: Wed, 7 Aug 2024 15:20:54 -0700 Subject: [PATCH 7/7] bugfix: handling of missing campets in create_ctedet_rowbyrow_pdf() --- py/desispec/night_qa.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/desispec/night_qa.py b/py/desispec/night_qa.py index 88f821516..36975e37e 100644 --- a/py/desispec/night_qa.py +++ b/py/desispec/night_qa.py @@ -1082,7 +1082,7 @@ def create_ctedet_rowbyrow_pdf(outpdf, night, prod, ctedet_expid, nproc, dpi=100 ctedet_mydicts = _read_ctedet(night, prod, ctedet_expid, nproc) # AR compute the row-by-row model - ims = [ctedet_mydict["img"] for ctedet_mydict in ctedet_mydicts] + ims = [ctedet_mydict["img"] if ctedet_mydict is not None else None for ctedet_mydict in ctedet_mydicts] mods = _compute_rowbyrow(ims, nproc) # AR plotting