#!/usr/bin/env python3 """ Solar Orbiter / SPICE downloader and data quality check tool. For every matching SPICE L2 FITS file, this script creates TWO figures: 1) SIMPLE SPECTRAL SUM S(x, y) = sum_lambda I_lambda(x, y, lambda) 2) PHYSICAL WAVELENGTH INTEGRAL I(x, y) = sum_lambda I_lambda(x, y, lambda) * Delta_lambda Both figures use logarithmic display scaling by default. The QC extension also creates, for every observational spectral window, a multi-panel diagnostic figure and a plain-text report. Quantitative QC is restricted by default to 120 <= y <= 700. The QC figure includes a dedicated panel containing: 1) average spectrum over the whole science region, 2) spectrum at the strongest integrated-intensity spatial pixel, 3) spectrum at the weakest integrated-intensity spatial pixel. Archive credentials can be read from pass.txt containing USERNAME and PASSWORD assignments. Dependencies ------------ pip install numpy matplotlib astropy requests beautifulsoup4 """ from __future__ import annotations import getpass import re from datetime import datetime, timedelta from pathlib import Path from urllib.parse import urljoin import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LogNorm import requests from bs4 import BeautifulSoup import astropy.units as u from astropy.io import fits # ============================================================ # USER SETTINGS # ============================================================ START_TIME = "2022-01-05 09:00:00" END_TIME = "2022-01-05 10:00:00" LEVEL = 2 PRODUCT = "spice-n-ras" SERVER_ROOT = "https://astro-sdc-db.uio.no/vol/spice/fits" DOWNLOAD_DIR = Path("./spice_data") PLOT_DIR = Path("./spice_plots") QC_DIR = Path("./spice_qc") # ============================================================ # AUTHENTICATION # ============================================================ # Credentials are read from a small local text file by default. # # Example pass.txt: # USERNAME = "abc" # PASSWORD = "cde" # # For safety, keep pass.txt outside version control. CREDENTIALS_FILE = Path("./pass.txt") READ_CREDENTIALS_FROM_FILE = True ASK_FOR_LOGIN_IF_FILE_MISSING = True # Optional hard-coded fallback. Normally keep these as None. USERNAME = None PASSWORD = None # ============================================================ # DOWNLOAD SETTINGS # ============================================================ # Existing non-empty files are reused. SKIP_EXISTING_FILES = True # ============================================================ # PLOT SETTINGS # ============================================================ CMAP = "viridis" # Logarithmic display scaling. USE_LOG_SCALE = True # Percentiles used separately for every spectral window. PERCENTILE_MIN = 1.0 PERCENTILE_MAX = 99.8 PLOT_DPI = 180 # Number of panels per row. NCOLS = 3 # ============================================================ # QUALITY-CONTROL SETTINGS # ============================================================ # Restrict quantitative diagnostics to the useful part of the slit. # Rows outside this interval are still shown in the diagnostic image, # but are excluded from min/max searches, averages, saturation tests, # NaN statistics used for QC status, and row/column outlier tests. USE_Y_SCIENCE_MASK = True Y_VALID_MIN = 120 Y_VALID_MAX = 700 # A row/column must have at least this fraction of finite samples to be # eligible for the min/max intensity selection. MIN_PROFILE_VALID_FRACTION = 0.95 # Robust row/column anomaly threshold, expressed in median absolute # deviations (MAD). OUTLIER_MAD_THRESHOLD = 6.0 # Conservative heuristic saturation detector. A spatial spectrum is # flagged only if it is among the brightest spectra and contains a # flat top of nearly identical adjacent spectral samples. These are # CANDIDATES, not instrument-confirmed saturation flags. SATURATION_BRIGHT_PERCENTILE = 99.9 SATURATION_MIN_PLATEAU_BINS = 2 SATURATION_RELATIVE_TOLERANCE = 1e-6 SATURATION_ABSOLUTE_TOLERANCE = 0.0 # QC status thresholds for non-finite samples inside the science region. QC_WARN_NONFINITE_FRACTION = 1e-3 QC_FAIL_NONFINITE_FRACTION = 5e-2 # Create the new diagnostic products in addition to the two existing # integrated-intensity figures. RUN_QUALITY_CONTROL = True # ============================================================ # TIME UTILITIES # ============================================================ def parse_time(value: str | datetime) -> datetime: if isinstance(value, datetime): return value value = value.replace("T", " ") for fmt in ( "%Y-%m-%d %H:%M:%S", "%Y-%m-%d %H:%M", "%Y-%m-%d", ): try: return datetime.strptime(value, fmt) except ValueError: pass raise ValueError(f"Cannot understand time: {value}") def iter_days(start: datetime, end: datetime): day = datetime(start.year, start.month, start.day) last = datetime(end.year, end.month, end.day) while day <= last: yield day day += timedelta(days=1) # ============================================================ # HTTP / AUTHENTICATION # ============================================================ def create_session(username=None, password=None): session = requests.Session() session.headers.update({ "User-Agent": ( "Mozilla/5.0 " "Solar-Orbiter-SPICE-Python-Downloader/2.0" ) }) if username and password: session.auth = (username, password) return session def read_credentials_file(path): """Read USERNAME and PASSWORD assignments from a simple text file. Expected format, for example: USERNAME = "abc" PASSWORD = "cde" The file is parsed as text; it is NOT executed. """ path = Path(path) # Prefer a pass.txt placed next to this script. If the configured path is # absolute, use it directly. if not path.is_absolute(): script_candidate = Path(__file__).resolve().parent / path cwd_candidate = Path.cwd() / path if script_candidate.exists(): path = script_candidate elif cwd_candidate.exists(): path = cwd_candidate else: path = script_candidate if not path.exists(): raise FileNotFoundError(f"Credentials file not found: {path}") content = path.read_text(encoding="utf-8") values = {} assignment_re = re.compile( r"^\s*(USERNAME|PASSWORD)\s*=\s*(['\"])(.*?)\2\s*(?:#.*)?$" ) for line in content.splitlines(): match = assignment_re.match(line) if match: values[match.group(1)] = match.group(3) username = values.get("USERNAME") password = values.get("PASSWORD") if not username or not password: raise ValueError( f"{path} must contain both USERNAME = \"...\" and PASSWORD = \"...\"." ) return username, password, path def get_credentials(): """Obtain SPICE archive credentials, preferring pass.txt.""" username = USERNAME password = PASSWORD if READ_CREDENTIALS_FROM_FILE: try: username, password, credentials_path = read_credentials_file( CREDENTIALS_FILE ) print(f"SPICE credentials loaded from: {credentials_path}") print(f"Username: {username}") # Never print the password. return username, password except Exception as exc: if not ASK_FOR_LOGIN_IF_FILE_MISSING: raise print() print(f"WARNING: Could not read credentials file: {exc}") print("Falling back to interactive login.") if username and password: return username, password print() print("SPICE archive authentication") print("----------------------------") username = input("Username: ").strip() password = getpass.getpass("Password: ") return username, password # ============================================================ # ARCHIVE SEARCH # ============================================================ SPICE_TIME_RE = re.compile(r"_(\d{8}T\d{6})_") def observation_time_from_filename(filename: str): match = SPICE_TIME_RE.search(filename) if not match: return None return datetime.strptime(match.group(1), "%Y%m%dT%H%M%S") def directory_url(day: datetime, level: int = 2): return ( f"{SERVER_ROOT}/" f"level{level}/" f"{day:%Y/%m/%d}/" ) def get_directory_files(url, session, timeout=60): print(f"Scanning: {url}") response = session.get(url, timeout=timeout) if response.status_code == 401: raise RuntimeError( "HTTP 401 Unauthorized. Check SPICE username/password." ) if response.status_code == 403: raise RuntimeError( "HTTP 403 Forbidden. The account may not have permission " "to access this directory." ) response.raise_for_status() soup = BeautifulSoup(response.text, "html.parser") urls = [] for link in soup.find_all("a", href=True): href = link["href"] if href.lower().endswith((".fits", ".fits.gz")): urls.append(urljoin(url, href)) return urls def find_spice_files( start_time, end_time, session, level=2, product="spice-n-ras", ): start = parse_time(start_time) end = parse_time(end_time) if end < start: raise ValueError("END_TIME must be later than START_TIME.") matches = [] for day in iter_days(start, end): url = directory_url(day, level=level) try: files = get_directory_files(url, session=session) except Exception as exc: print(f"WARNING: Could not read {url}") print(f" {exc}") continue for file_url in files: filename = Path(file_url).name if f"solo_L{level}_".lower() not in filename.lower(): continue if product and product.lower() not in filename.lower(): continue obs_time = observation_time_from_filename(filename) if obs_time is None: continue if start <= obs_time <= end: matches.append(file_url) return sorted(set(matches)) # ============================================================ # DOWNLOAD # ============================================================ def download_file( url, output_directory, session, skip_existing=True, ): output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) destination = output_directory / Path(url).name # Reuse an existing non-empty file. if destination.exists(): size = destination.stat().st_size if skip_existing and size > 0: print( f"Already downloaded: {destination.name} " f"({size / 1024 / 1024:.2f} MB) -- reusing it." ) return destination if size == 0: print( f"Existing file is empty: {destination.name}; " "downloading it again." ) print(f"Downloading: {url}") with session.get(url, stream=True, timeout=180) as response: if response.status_code == 401: raise RuntimeError( "HTTP 401 Unauthorized while downloading." ) if response.status_code == 403: raise RuntimeError( "HTTP 403 Forbidden while downloading." ) response.raise_for_status() total = int(response.headers.get("content-length", 0)) received = 0 # Download first to a temporary file so a failed transfer does # not leave a file that looks complete on the next run. temporary = destination.with_suffix(destination.suffix + ".part") try: with open(temporary, "wb") as f: for chunk in response.iter_content( chunk_size=1024 * 1024 ): if not chunk: continue f.write(chunk) received += len(chunk) if total: pct = 100.0 * received / total print( f"\r {pct:6.2f} %", end="", flush=True, ) if total: print() temporary.replace(destination) except Exception: if temporary.exists(): temporary.unlink() raise print(f"Saved: {destination}") return destination def download_spice_range( start_time, end_time, level=2, product="spice-n-ras", output_directory="./spice_data", username=None, password=None, ): session = create_session( username=username, password=password, ) urls = find_spice_files( start_time=start_time, end_time=end_time, session=session, level=level, product=product, ) print() print(f"Found {len(urls)} matching SPICE file(s).") paths = [] for url in urls: try: path = download_file( url, output_directory, session=session, skip_existing=SKIP_EXISTING_FILES, ) paths.append(path) except Exception as exc: print(f"ERROR downloading {url}") print(f" {exc}") return paths # ============================================================ # FITS / SPICE AXES # ============================================================ def fits_axis_to_numpy_axis(fits_axis: int, ndim: int): """ FITS axis numbering is reversed relative to NumPy indexing. """ return ndim - fits_axis def find_spectral_axis(header, ndim): """ Identify the wavelength/dispersion axis from FITS WCS metadata. Returns ------- numpy_axis, fits_axis """ for fits_axis in range(1, ndim + 1): ctype = str( header.get(f"CTYPE{fits_axis}", "") ).upper() cname = str( header.get(f"CNAME{fits_axis}", "") ).upper() text = f"{ctype} {cname}" if any( token in text for token in ( "WAVE", "AWAV", "FREQ", "SPECT", "DISP", ) ): return ( fits_axis_to_numpy_axis(fits_axis, ndim), fits_axis, ) # Standard SPICE fallback: # FITS axis 3 = dispersion/wavelength. if ndim >= 3: fits_axis = 3 return fits_axis_to_numpy_axis(fits_axis, ndim), fits_axis raise RuntimeError("Could not identify spectral axis.") def spectral_pixel_width_nm(header, fits_spectral_axis): """ Return absolute wavelength step in nm for a linear WCS. """ cdelt_key = f"CDELT{fits_spectral_axis}" if cdelt_key in header: delta = abs(float(header[cdelt_key])) else: cd_key = ( f"CD{fits_spectral_axis}_" f"{fits_spectral_axis}" ) if cd_key not in header: return None delta = abs(float(header[cd_key])) unit_string = header.get( f"CUNIT{fits_spectral_axis}", "", ) if not unit_string: return None try: wavelength_unit = u.Unit(unit_string) return float( (delta * wavelength_unit).to_value(u.nm) ) except Exception: return None def is_observational_hdu(hdu): if hdu.data is None: return False if not isinstance(hdu.data, np.ndarray): return False if hdu.data.ndim < 3: return False header = hdu.header if "WIN_TYPE" in header: return True btype = str(header.get("BTYPE", "")).lower() if "radiance" in btype or "intensity" in btype: return True # Keep this fallback because many SPICE science windows have EXTNAME. if str(header.get("EXTNAME", "")).strip(): return True return False # ============================================================ # CALCULATE BOTH MAP TYPES # ============================================================ def calculate_intensity_maps(hdu): """ Calculate both requested versions from one spectral window. Returns ------- simple_sum : 2-D ndarray Sum over the spectral pixels: sum_lambda I_lambda physical_integral : 2-D ndarray or None Wavelength integral: sum_lambda I_lambda * Delta_lambda None if Delta_lambda cannot be determined. delta_nm : float or None Spectral pixel width in nm. numpy_spectral_axis : int NumPy axis that was integrated. fits_spectral_axis : int Corresponding FITS axis number. """ data = np.asarray(hdu.data, dtype=np.float64) data[~np.isfinite(data)] = np.nan numpy_spectral_axis, fits_spectral_axis = ( find_spectral_axis( hdu.header, data.ndim, ) ) # -------------------------------------------------------- # VERSION 1: SIMPLE SUM # -------------------------------------------------------- simple_sum = np.nansum( data, axis=numpy_spectral_axis, ) simple_sum = np.squeeze(simple_sum) # Standard SPICE n-ras usually has a singleton time dimension. # If a non-singleton extra dimension remains, average it. while simple_sum.ndim > 2: simple_sum = np.nanmean( simple_sum, axis=0, ) # -------------------------------------------------------- # VERSION 2: PHYSICAL WAVELENGTH INTEGRAL # -------------------------------------------------------- delta_nm = spectral_pixel_width_nm( hdu.header, fits_spectral_axis, ) if delta_nm is None: physical_integral = None else: physical_integral = simple_sum * delta_nm return ( simple_sum, physical_integral, delta_nm, numpy_spectral_axis, fits_spectral_axis, ) # ============================================================ # TITLES / UNITS # ============================================================ def clean_window_title(header, hdu_name): candidates = [ header.get("EXTNAME"), header.get("LINE_ID"), header.get("WINTAB"), hdu_name, ] for value in candidates: if value: title = str(value).strip() title = re.sub( r"^\[[^\]]+\]\s*", "", title, ) return title return "SPICE spectral window" def physical_output_unit(header): """ Best-effort colorbar unit after multiplying spectral radiance by wavelength in nm. """ bunit = str(header.get("BUNIT", "")).strip() if not bunit: return "" result = bunit for pattern in ( r"/\s*nm\b", r"\bnm\s*\^\s*-1\b", r"\bnm\s*\*\*\s*-1\b", r"\bnm-1\b", ): result = re.sub( pattern, "", result, flags=re.IGNORECASE, ) return " ".join(result.split()) def simple_sum_output_unit(header): bunit = str(header.get("BUNIT", "")).strip() if not bunit: return "sum over spectral pixels" # A simple numeric sum has no d-lambda multiplication, so this is # intentionally described as a spectral-pixel sum rather than as # a physical wavelength integral. return f"spectral-pixel sum [{bunit}]" # ============================================================ # DISPLAY NORMALIZATION # ============================================================ def display_normalization( image, use_log=True, pmin=1.0, pmax=99.8, ): """ Return image_for_plot, norm, vmin, vmax. For logarithmic display, non-positive values are masked. """ finite = image[np.isfinite(image)] if finite.size == 0: raise ValueError("Image contains no finite values.") if use_log: positive = image[ np.isfinite(image) & (image > 0) ] if positive.size == 0: raise ValueError( "Image contains no positive values for LogNorm." ) vmin, vmax = np.nanpercentile( positive, [pmin, pmax], ) if not np.isfinite(vmin) or vmin <= 0: vmin = np.nanmin(positive) if not np.isfinite(vmax) or vmax <= vmin: vmax = np.nanmax(positive) masked = np.ma.masked_where( (~np.isfinite(image)) | (image <= 0), image, ) return ( masked, LogNorm( vmin=float(vmin), vmax=float(vmax), clip=True, ), float(vmin), float(vmax), ) vmin, vmax = np.nanpercentile( finite, [pmin, pmax], ) return image, None, float(vmin), float(vmax) # ============================================================ # GENERIC MULTI-PANEL PLOTTER # ============================================================ def save_multi_panel_plot( filename, windows, output_directory, suffix, figure_title, use_log_scale=True, cmap="viridis", pmin=1.0, pmax=99.8, ): if not windows: print( f"WARNING: no windows available for {figure_title}" ) return None n_windows = len(windows) nrows = int(np.ceil(n_windows / NCOLS)) fig, axes = plt.subplots( nrows, NCOLS, figsize=(5.4 * NCOLS, 3.8 * nrows), squeeze=False, constrained_layout=True, ) axes = axes.ravel() for ax, window in zip(axes, windows): image = window["image"] try: image_for_plot, norm, vmin, vmax = ( display_normalization( image, use_log=use_log_scale, pmin=pmin, pmax=pmax, ) ) if norm is not None: im = ax.imshow( image_for_plot, origin="lower", aspect="auto", cmap=cmap, norm=norm, interpolation="nearest", ) else: im = ax.imshow( image_for_plot, origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax, interpolation="nearest", ) except ValueError as exc: print( f"WARNING: {window['title']}: {exc}; " "using unscaled linear display." ) im = ax.imshow( image, origin="lower", aspect="auto", cmap=cmap, interpolation="nearest", ) ax.set_title( window["title"], fontsize=11, ) ax.set_xlabel("Raster step [pixel]") ax.set_ylabel("Slit position [pixel]") cbar = fig.colorbar( im, ax=ax, pad=0.02, ) if window.get("unit"): cbar.set_label( window["unit"], fontsize=8, ) for ax in axes[n_windows:]: ax.axis("off") scale_text = ( "logarithmic display" if use_log_scale else "linear display" ) fig.suptitle( f"{filename.name}\n" f"{figure_title} — {scale_text}", fontsize=13, ) output_directory = Path(output_directory) output_directory.mkdir( parents=True, exist_ok=True, ) output_file = ( output_directory / f"{filename.stem}_{suffix}.png" ) fig.savefig( output_file, dpi=PLOT_DPI, bbox_inches="tight", ) plt.close(fig) print(f"Saved plot: {output_file}") return output_file # ============================================================ # PROCESS ONE FITS FILE AND CREATE TWO FIGURES # ============================================================ def plot_both_intensity_versions( filename, output_directory="./spice_plots", use_log_scale=True, cmap="viridis", pmin=1.0, pmax=99.8, ): """ Create two PNG files for one SPICE FITS observation: *_simple_sum_log.png sum(I_lambda) over the spectral pixels. *_wavelength_integrated_log.png sum(I_lambda) * Delta_lambda. """ filename = Path(filename) simple_windows = [] physical_windows = [] print() print(f"Opening FITS: {filename}") with fits.open( filename, memmap=True, ) as hdul: print(f"Number of FITS HDUs: {len(hdul)}") for hdu_index, hdu in enumerate(hdul): if not is_observational_hdu(hdu): continue print( f" HDU {hdu_index}: {hdu.name}, " f"shape={getattr(hdu.data, 'shape', None)}" ) try: ( simple_sum, physical_integral, delta_nm, numpy_spectral_axis, fits_spectral_axis, ) = calculate_intensity_maps(hdu) except Exception as exc: print( f" Skipping HDU {hdu_index}: {exc}" ) continue if simple_sum.ndim != 2: print( f" Skipping HDU {hdu_index}: " f"result has shape {simple_sum.shape}" ) continue title = clean_window_title( hdu.header, hdu.name, ) print( f" spectral axis: NumPy " f"{numpy_spectral_axis}, FITS {fits_spectral_axis}" ) if delta_nm is not None: print( f" Delta lambda: " f"{delta_nm:.6g} nm" ) else: print( " Delta lambda: unavailable" ) # ---------------- SIMPLE SUM ---------------- simple_windows.append({ "image": simple_sum, "title": title, "unit": simple_sum_output_unit( hdu.header ), }) # ------------- PHYSICAL INTEGRAL ------------ if physical_integral is not None: physical_windows.append({ "image": physical_integral, "title": title, "unit": physical_output_unit( hdu.header ), }) else: print( f" WARNING: cannot calculate physical " f"wavelength integral for {title} because " "the wavelength step could not be determined." ) # -------------------------------------------------------- # Save simple-sum figure # -------------------------------------------------------- scale_suffix = "log" if use_log_scale else "linear" simple_plot = save_multi_panel_plot( filename=filename, windows=simple_windows, output_directory=output_directory, suffix=f"simple_sum_{scale_suffix}", figure_title=( "Simple spectral-pixel sum: " "Σ Iλ" ), use_log_scale=use_log_scale, cmap=cmap, pmin=pmin, pmax=pmax, ) # -------------------------------------------------------- # Save wavelength-integrated figure # -------------------------------------------------------- physical_plot = save_multi_panel_plot( filename=filename, windows=physical_windows, output_directory=output_directory, suffix=f"wavelength_integrated_{scale_suffix}", figure_title=( "Wavelength-integrated intensity: " "Σ Iλ Δλ" ), use_log_scale=use_log_scale, cmap=cmap, pmin=pmin, pmax=pmax, ) return simple_plot, physical_plot # ============================================================ # SPICE QUALITY-CONTROL / DIAGNOSTIC TOOLS # ============================================================ def canonical_spectral_cube(hdu, replace_nonfinite=False): """Return a SPICE window as a canonical [lambda, y, x] cube. The spectral axis is identified from the FITS/WCS metadata. Singleton dimensions are removed. If an additional non-singleton dimension remains (for example time), it is averaged until only [lambda, y, x] remains. Parameters ---------- hdu : astropy.io.fits.ImageHDU replace_nonfinite : bool If True, +/-Inf are converted to NaN after the cube has been put in canonical order. """ data = np.asarray(hdu.data, dtype=np.float64) spectral_axis, fits_spectral_axis = find_spectral_axis( hdu.header, data.ndim ) cube = np.moveaxis(data, spectral_axis, 0) # Remove singleton non-spectral dimensions without touching lambda. squeeze_axes = tuple( ax for ax in range(1, cube.ndim) if cube.shape[ax] == 1 ) if squeeze_axes: cube = np.squeeze(cube, axis=squeeze_axes) # Preserve the last two dimensions as the spatial image. If another # dimension remains, average it. This matches the behavior of the # existing integrated-map routine for standard SPICE n-ras products. while cube.ndim > 3: cube = np.nanmean(cube, axis=1) if cube.ndim != 3: raise RuntimeError( f"Could not obtain [lambda, y, x] cube; got shape {cube.shape}." ) if replace_nonfinite: cube = cube.copy() cube[~np.isfinite(cube)] = np.nan return cube, fits_spectral_axis def spectral_coordinates_nm(header, fits_spectral_axis, n_lambda): """Best-effort wavelength coordinate in nm; otherwise pixel coordinate.""" crval = header.get(f"CRVAL{fits_spectral_axis}") crpix = header.get(f"CRPIX{fits_spectral_axis}") cdelt = header.get(f"CDELT{fits_spectral_axis}") if cdelt is None: cdelt = header.get( f"CD{fits_spectral_axis}_{fits_spectral_axis}" ) unit_string = header.get(f"CUNIT{fits_spectral_axis}", "") if crval is None or crpix is None or cdelt is None or not unit_string: return np.arange(n_lambda, dtype=float), "Spectral pixel" try: pix_fits = np.arange(n_lambda, dtype=float) + 1.0 values = float(crval) + (pix_fits - float(crpix)) * float(cdelt) unit = u.Unit(unit_string) values_nm = (values * unit).to_value(u.nm) return np.asarray(values_nm, dtype=float), "Wavelength [nm]" except Exception: return np.arange(n_lambda, dtype=float), "Spectral pixel" def science_y_bounds(ny, y_min=Y_VALID_MIN, y_max=Y_VALID_MAX): """Return inclusive science-region bounds clipped to the actual image.""" if not USE_Y_SCIENCE_MASK: return 0, ny - 1 lo = max(0, int(y_min)) hi = min(ny - 1, int(y_max)) if lo > hi: raise ValueError( f"Invalid science Y interval {y_min}..{y_max} for ny={ny}." ) return lo, hi def integrated_map_preserve_missing(cube): """Spectral sum with fully invalid spectra kept as NaN, not zero.""" finite = np.isfinite(cube) valid_count = finite.sum(axis=0) image = np.nansum(np.where(finite, cube, np.nan), axis=0) image = np.asarray(image, dtype=float) image[valid_count == 0] = np.nan return image, valid_count def safe_nanargextreme(values, eligible, mode="min"): """Index of min/max value among eligible finite positions, else None.""" values = np.asarray(values, dtype=float) eligible = np.asarray(eligible, dtype=bool) & np.isfinite(values) if not np.any(eligible): return None candidate_indices = np.flatnonzero(eligible) candidate_values = values[candidate_indices] if mode == "min": return int(candidate_indices[np.argmin(candidate_values)]) if mode == "max": return int(candidate_indices[np.argmax(candidate_values)]) raise ValueError("mode must be 'min' or 'max'") def robust_outlier_mask(values, valid_mask=None, threshold=OUTLIER_MAD_THRESHOLD): """Two-sided robust outlier mask using median and MAD.""" values = np.asarray(values, dtype=float) if valid_mask is None: valid_mask = np.isfinite(values) else: valid_mask = np.asarray(valid_mask, dtype=bool) & np.isfinite(values) result = np.zeros(values.shape, dtype=bool) sample = values[valid_mask] if sample.size < 5: return result med = np.nanmedian(sample) mad = np.nanmedian(np.abs(sample - med)) if not np.isfinite(mad) or mad == 0: return result robust_sigma = 1.4826 * mad result[valid_mask] = ( np.abs(values[valid_mask] - med) > threshold * robust_sigma ) return result def saturation_candidate_mask( cube_science, bright_percentile=SATURATION_BRIGHT_PERCENTILE, min_plateau_bins=SATURATION_MIN_PLATEAU_BINS, rtol=SATURATION_RELATIVE_TOLERANCE, atol=SATURATION_ABSOLUTE_TOLERANCE, ): """Return a conservative 2-D mask of possible saturated spectra. This is deliberately a heuristic. A spectrum is flagged only when its peak is in the very bright tail of all finite science-region samples and at least ``min_plateau_bins`` adjacent spectral bins are nearly identical to that peak. This is intended to locate flat-topped/clipped spectra. """ finite_values = cube_science[np.isfinite(cube_science)] spatial_shape = cube_science.shape[1:] if finite_values.size == 0: return np.zeros(spatial_shape, dtype=bool), np.nan bright_threshold = float( np.nanpercentile(finite_values, bright_percentile) ) with np.errstate(all="ignore"): peak = np.nanmax(cube_science, axis=0) bright = np.isfinite(peak) & (peak >= bright_threshold) tolerance = np.maximum(np.abs(peak) * rtol, atol) near_peak = ( np.isfinite(cube_science) & np.isfinite(peak)[None, :, :] & (np.abs(cube_science - peak[None, :, :]) <= tolerance[None, :, :]) ) run = np.zeros(spatial_shape, dtype=np.int16) max_run = np.zeros(spatial_shape, dtype=np.int16) for k in range(near_peak.shape[0]): run = np.where(near_peak[k], run + 1, 0) max_run = np.maximum(max_run, run) mask = bright & (max_run >= int(min_plateau_bins)) return mask, bright_threshold def qc_status(nonfinite_fraction, saturation_count, anomalous_rows, anomalous_cols): """Simple configurable PASS/WARNING/FAIL summary.""" if not np.isfinite(nonfinite_fraction): return "FAIL" if nonfinite_fraction >= QC_FAIL_NONFINITE_FRACTION: return "FAIL" if ( nonfinite_fraction >= QC_WARN_NONFINITE_FRACTION or saturation_count > 0 or anomalous_rows > 0 or anomalous_cols > 0 ): return "WARNING" return "PASS" def analyze_spice_window_qc(hdu, hdu_index): """Analyze one observational SPICE HDU / spectral window independently. All spatial extrema and spectra in this function belong only to this HDU; they are never calculated globally across the full multi-window FITS file. """ raw_cube, fits_spectral_axis = canonical_spectral_cube( hdu, replace_nonfinite=False ) cube = raw_cube.copy() cube[~np.isfinite(cube)] = np.nan n_lambda, ny, nx = cube.shape y0, y1 = science_y_bounds(ny) ys = slice(y0, y1 + 1) title = clean_window_title(hdu.header, hdu.name) wavelength, wavelength_label = spectral_coordinates_nm( hdu.header, fits_spectral_axis, n_lambda ) image, valid_count_map = integrated_map_preserve_missing(cube) nonfinite_fraction_map = 1.0 - valid_count_map / float(n_lambda) cube_science_raw = raw_cube[:, ys, :] cube_science = cube[:, ys, :] image_science = image[ys, :] finite_science = np.isfinite(cube_science_raw) total_science_samples = int(cube_science_raw.size) finite_count = int(finite_science.sum()) nan_count = int(np.isnan(cube_science_raw).sum()) posinf_count = int(np.isposinf(cube_science_raw).sum()) neginf_count = int(np.isneginf(cube_science_raw).sum()) nonfinite_count = total_science_samples - finite_count nonfinite_fraction = ( nonfinite_count / total_science_samples if total_science_samples else np.nan ) finite_values = cube_science_raw[finite_science] zero_count = int(np.count_nonzero(finite_values == 0)) negative_count = int(np.count_nonzero(finite_values < 0)) science_valid_count = np.isfinite(cube_science).sum(axis=0) fully_invalid_spectra = int(np.count_nonzero(science_valid_count == 0)) partially_invalid_spectra = int( np.count_nonzero( (science_valid_count > 0) & (science_valid_count < n_lambda) ) ) # X profile: total intensity over valid science Y rows. x_total = np.nansum(image_science, axis=0) x_finite_fraction = np.mean(np.isfinite(image_science), axis=0) x_eligible = x_finite_fraction >= MIN_PROFILE_VALID_FRACTION x_min = safe_nanargextreme(x_total, x_eligible, "min") x_max = safe_nanargextreme(x_total, x_eligible, "max") # Y profile: total intensity over X, restricted to science rows. y_total_science = np.nansum(image_science, axis=1) y_finite_fraction = np.mean(np.isfinite(image_science), axis=1) y_eligible_science = y_finite_fraction >= MIN_PROFILE_VALID_FRACTION y_min_local = safe_nanargextreme( y_total_science, y_eligible_science, "min" ) y_max_local = safe_nanargextreme( y_total_science, y_eligible_science, "max" ) y_min = None if y_min_local is None else int(y0 + y_min_local) y_max = None if y_max_local is None else int(y0 + y_max_local) # Average spectra requested by the user. full_spectrum = np.nanmean(cube_science, axis=(1, 2)) x_min_spectrum = ( None if x_min is None else np.nanmean(cube_science[:, :, x_min], axis=1) ) x_max_spectrum = ( None if x_max is None else np.nanmean(cube_science[:, :, x_max], axis=1) ) y_min_spectrum = ( None if y_min is None else np.nanmean(cube[:, y_min, :], axis=1) ) y_max_spectrum = ( None if y_max is None else np.nanmean(cube[:, y_max, :], axis=1) ) # -------------------------------------------------------- # Point-by-point extrema for THIS spectral window only. # # IMPORTANT: # These extrema are recalculated independently for every observational # HDU / spectral-range map. Nothing is shared between spectral windows. # # Search only inside the accepted science-Y interval. A spatial pixel # is eligible whenever its integrated-map value is finite and at least # one spectral sample is finite. We intentionally DO NOT apply the # stricter MIN_PROFILE_VALID_FRACTION here; that threshold is suitable # for row/column averages, but in some SPICE windows it can reject every # individual spatial pixel and therefore suppress the requested point # minimum/maximum. # -------------------------------------------------------- point_valid_fraction = science_valid_count / float(n_lambda) point_eligible = ( np.isfinite(image_science) & (science_valid_count > 0) ) point_min_x = point_min_y = None point_max_x = point_max_y = None point_min_intensity = None point_max_intensity = None point_min_valid_fraction = None point_max_valid_fraction = None point_min_spectrum = None point_max_spectrum = None if np.any(point_eligible): # Use this HDU's own integrated-intensity map only. eligible_values = np.where(point_eligible, image_science, np.nan) flat_min = int(np.nanargmin(eligible_values)) flat_max = int(np.nanargmax(eligible_values)) point_min_y_local, point_min_x = np.unravel_index( flat_min, image_science.shape ) point_max_y_local, point_max_x = np.unravel_index( flat_max, image_science.shape ) point_min_y = int(y0 + point_min_y_local) point_max_y = int(y0 + point_max_y_local) point_min_x = int(point_min_x) point_max_x = int(point_max_x) point_min_intensity = float( image_science[point_min_y_local, point_min_x] ) point_max_intensity = float( image_science[point_max_y_local, point_max_x] ) point_min_valid_fraction = float( point_valid_fraction[point_min_y_local, point_min_x] ) point_max_valid_fraction = float( point_valid_fraction[point_max_y_local, point_max_x] ) # These are the spectra of the exact spatial pixels selected from # this same spectral-range integrated map. point_min_spectrum = cube[:, point_min_y, point_min_x].copy() point_max_spectrum = cube[:, point_max_y, point_max_x].copy() saturation_science, bright_threshold = saturation_candidate_mask( cube_science ) saturation_full = np.zeros((ny, nx), dtype=bool) saturation_full[ys, :] = saturation_science saturation_count = int(np.count_nonzero(saturation_science)) x_outliers = robust_outlier_mask(x_total, x_eligible) y_outliers_science = robust_outlier_mask( y_total_science, y_eligible_science ) anomalous_cols = int(np.count_nonzero(x_outliers)) anomalous_rows = int(np.count_nonzero(y_outliers_science)) def finite_stat(func): if finite_values.size == 0: return None value = float(func(finite_values)) return value if np.isfinite(value) else None summary = { "hdu_index": int(hdu_index), "hdu_name": str(hdu.name), "title": title, "cube_shape_lambda_y_x": [int(v) for v in cube.shape], "science_y_min": int(y0), "science_y_max": int(y1), "sample_statistics": { "total_samples_science_region": total_science_samples, "finite_samples": finite_count, "nan_samples": nan_count, "positive_inf_samples": posinf_count, "negative_inf_samples": neginf_count, "nonfinite_samples": nonfinite_count, "nonfinite_fraction": float(nonfinite_fraction), "zero_samples": zero_count, "negative_samples": negative_count, "minimum": finite_stat(np.nanmin), "maximum": finite_stat(np.nanmax), "median": finite_stat(np.nanmedian), "mean": finite_stat(np.nanmean), "std": finite_stat(np.nanstd), }, "spectrum_completeness": { "spatial_spectra_science_region": int(science_valid_count.size), "fully_invalid_spectra": fully_invalid_spectra, "partially_invalid_spectra": partially_invalid_spectra, }, "extrema": { "x_min": x_min, "x_max": x_max, "y_min": y_min, "y_max": y_max, "minimum_profile_valid_fraction": float( MIN_PROFILE_VALID_FRACTION ), }, "point_extrema": { "minimum": { "x": point_min_x, "y": point_min_y, "integrated_intensity": point_min_intensity, "valid_spectral_fraction": point_min_valid_fraction, }, "maximum": { "x": point_max_x, "y": point_max_y, "integrated_intensity": point_max_intensity, "valid_spectral_fraction": point_max_valid_fraction, }, "selection": ( "Minimum and maximum finite spatial pixels are calculated " "independently from THIS spectral-window integrated map, " "inside the science Y interval. Fully invalid spectra are " "excluded." ), }, "saturation_candidates": { "method": ( "bright-tail + adjacent near-identical flat-top heuristic" ), "bright_percentile": float(SATURATION_BRIGHT_PERCENTILE), "bright_threshold": float(bright_threshold) if np.isfinite(bright_threshold) else None, "minimum_plateau_bins": int(SATURATION_MIN_PLATEAU_BINS), "candidate_spatial_pixels": saturation_count, "note": ( "Candidates are heuristic and are not instrument-confirmed " "saturation flags." ), }, "row_column_anomalies": { "mad_threshold": float(OUTLIER_MAD_THRESHOLD), "anomalous_x_columns": anomalous_cols, "anomalous_y_rows": anomalous_rows, "x_indices": [int(v) for v in np.flatnonzero(x_outliers)], "y_indices": [ int(y0 + v) for v in np.flatnonzero(y_outliers_science) ], }, } arrays = { "image": image, "nonfinite_fraction_map": nonfinite_fraction_map, "saturation_mask": saturation_full, "wavelength": wavelength, "wavelength_label": wavelength_label, "x_total": x_total, "x_eligible": x_eligible, "x_outliers": x_outliers, "y_total_science": y_total_science, "y_eligible_science": y_eligible_science, "y_outliers_science": y_outliers_science, "full_spectrum": full_spectrum, "x_min_spectrum": x_min_spectrum, "x_max_spectrum": x_max_spectrum, "y_min_spectrum": y_min_spectrum, "y_max_spectrum": y_max_spectrum, "point_min_spectrum": point_min_spectrum, "point_max_spectrum": point_max_spectrum, "point_min_xy": ( None if point_min_x is None else (point_min_x, point_min_y) ), "point_max_xy": ( None if point_max_x is None else (point_max_x, point_max_y) ), "y0": y0, "y1": y1, } return summary, arrays def plot_spice_window_qc(filename, summary, arrays, output_directory): """Create one QC figure for a spectral window. Panels ------ 1. Integrated intensity map with QC selections and point extrema. 2. NaN / Inf fraction map. 3. Total intensity versus raster step. 4. Total intensity versus slit position. 5. Average spectra for the full field and min/max X/Y cuts. 6. NEW: full science-region average spectrum + spectra at the spatial pixels with minimum and maximum integrated intensity. """ output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) image = arrays["image"] nonfinite_fraction_map = arrays["nonfinite_fraction_map"] saturation_mask = arrays["saturation_mask"] wavelength = arrays["wavelength"] wavelength_label = arrays["wavelength_label"] y0 = arrays["y0"] y1 = arrays["y1"] # Four rows: two image panels, two 1-D profiles, old spectrum panel, # and the new point-extrema spectrum panel. fig = plt.figure(figsize=(16, 16), constrained_layout=True) gs = fig.add_gridspec( 4, 2, height_ratios=[1.25, 0.85, 1.0, 1.0], ) ax_map = fig.add_subplot(gs[0, 0]) ax_nan = fig.add_subplot(gs[0, 1]) ax_x = fig.add_subplot(gs[1, 0]) ax_y = fig.add_subplot(gs[1, 1]) ax_spec = fig.add_subplot(gs[2, :]) ax_point_spec = fig.add_subplot(gs[3, :]) # Main map. try: image_for_plot, norm, vmin, vmax = display_normalization( image, use_log=USE_LOG_SCALE, pmin=PERCENTILE_MIN, pmax=PERCENTILE_MAX, ) if norm is not None: im = ax_map.imshow( image_for_plot, origin="lower", aspect="auto", cmap=CMAP, norm=norm, interpolation="nearest", ) else: im = ax_map.imshow( image_for_plot, origin="lower", aspect="auto", cmap=CMAP, vmin=vmin, vmax=vmax, interpolation="nearest", ) except Exception: im = ax_map.imshow( image, origin="lower", aspect="auto", cmap=CMAP, interpolation="nearest", ) fig.colorbar(im, ax=ax_map, pad=0.02, label="Spectral-pixel sum") ax_map.set_title("Integrated intensity + QC selections") ax_map.set_xlabel("Raster step x [pixel]") ax_map.set_ylabel("Slit position y [pixel]") # Shade non-science slit edges but keep them visible. if y0 > 0: ax_map.axhspan(-0.5, y0 - 0.5, alpha=0.18, color="gray") if y1 < image.shape[0] - 1: ax_map.axhspan( y1 + 0.5, image.shape[0] - 0.5, alpha=0.18, color="gray", ) ax_map.axhline(y0, linestyle="--", linewidth=0.8, color="black") ax_map.axhline(y1, linestyle="--", linewidth=0.8, color="black") ext = summary["extrema"] if ext["x_min"] is not None: ax_map.axvline( ext["x_min"], linestyle=":", linewidth=1.0, color="cyan", label="x min", ) if ext["x_max"] is not None: ax_map.axvline( ext["x_max"], linestyle=":", linewidth=1.0, color="magenta", label="x max", ) if ext["y_min"] is not None: ax_map.axhline( ext["y_min"], linestyle=":", linewidth=1.0, color="cyan", label="y min", ) if ext["y_max"] is not None: ax_map.axhline( ext["y_max"], linestyle=":", linewidth=1.0, color="magenta", label="y max", ) # Mark the actual spatial pixels having minimum/maximum integrated # intensity inside the science region. pext = summary["point_extrema"] pmin = pext["minimum"] pmax = pext["maximum"] if pmin["x"] is not None and pmin["y"] is not None: ax_map.scatter( [pmin["x"]], [pmin["y"]], marker="v", s=85, facecolors="none", edgecolors="cyan", linewidths=1.6, label=f"point min ({pmin['x']},{pmin['y']})", zorder=8, ) if pmax["x"] is not None and pmax["y"] is not None: ax_map.scatter( [pmax["x"]], [pmax["y"]], marker="^", s=85, facecolors="none", edgecolors="magenta", linewidths=1.6, label=f"point max ({pmax['x']},{pmax['y']})", zorder=8, ) if np.any(saturation_mask): ax_map.contour( saturation_mask.astype(float), levels=[0.5], colors="red", linewidths=1.0, origin="lower", ) handles, labels = ax_map.get_legend_handles_labels() if handles: unique = dict(zip(labels, handles)) ax_map.legend( unique.values(), unique.keys(), fontsize=7, loc="upper right", ) # Non-finite map. nan_im = ax_nan.imshow( nonfinite_fraction_map, origin="lower", aspect="auto", vmin=0.0, vmax=1.0, interpolation="nearest", ) fig.colorbar( nan_im, ax=ax_nan, pad=0.02, label="Non-finite fraction per spectrum", ) ax_nan.set_title("NaN / Inf fraction map") ax_nan.set_xlabel("Raster step x [pixel]") ax_nan.set_ylabel("Slit position y [pixel]") if y0 > 0: ax_nan.axhspan(-0.5, y0 - 0.5, alpha=0.18, color="gray") if y1 < image.shape[0] - 1: ax_nan.axhspan( y1 + 0.5, image.shape[0] - 0.5, alpha=0.18, color="gray", ) # X total profile. x = np.arange(arrays["x_total"].size) ax_x.plot(x, arrays["x_total"], linewidth=1.0) ax_x.set_title("Total intensity vs raster step (science Y only)") ax_x.set_xlabel("Raster step x [pixel]") ax_x.set_ylabel("Summed intensity") if ext["x_min"] is not None: ax_x.axvline( ext["x_min"], linestyle="--", linewidth=0.9, color="cyan", ) if ext["x_max"] is not None: ax_x.axvline( ext["x_max"], linestyle="--", linewidth=0.9, color="magenta", ) out_x = np.flatnonzero(arrays["x_outliers"]) if out_x.size: ax_x.scatter( out_x, arrays["x_total"][out_x], marker="x", color="red", ) # Y total profile. y_values = np.arange(y0, y1 + 1) ax_y.plot(y_values, arrays["y_total_science"], linewidth=1.0) ax_y.set_title("Total intensity vs slit position (science Y only)") ax_y.set_xlabel("Slit position y [pixel]") ax_y.set_ylabel("Summed intensity") if ext["y_min"] is not None: ax_y.axvline( ext["y_min"], linestyle="--", linewidth=0.9, color="cyan", ) if ext["y_max"] is not None: ax_y.axvline( ext["y_max"], linestyle="--", linewidth=0.9, color="magenta", ) out_y_local = np.flatnonzero(arrays["y_outliers_science"]) if out_y_local.size: out_y = y0 + out_y_local ax_y.scatter( out_y, arrays["y_total_science"][out_y_local], marker="x", color="red", ) # Existing spectral diagnostic. ax_spec.plot( wavelength, arrays["full_spectrum"], linewidth=2.0, label="Full science-region average", ) spectral_series = [ ("x_min_spectrum", ext["x_min"], "x min"), ("x_max_spectrum", ext["x_max"], "x max"), ("y_min_spectrum", ext["y_min"], "y min"), ("y_max_spectrum", ext["y_max"], "y max"), ] for key, coordinate, label in spectral_series: spectrum = arrays[key] if spectrum is not None: ax_spec.plot( wavelength, spectrum, linewidth=1.0, label=f"{label} = {coordinate}", ) ax_spec.set_title( "Average spectra for full science region and min/max X/Y cuts" ) ax_spec.set_xlabel(wavelength_label) ax_spec.set_ylabel(str(hdu_unit_from_summary(summary))) ax_spec.legend(ncol=3, fontsize=8) ax_spec.grid(alpha=0.2) # -------------------------------------------------------- # NEW PANEL: # 1) whole science-region average spectrum # 2) spectrum at maximum integrated-intensity spatial pixel # 3) spectrum at minimum integrated-intensity spatial pixel # -------------------------------------------------------- ax_point_spec.plot( wavelength, arrays["full_spectrum"], linewidth=2.2, label="Whole science-region average", ) if arrays["point_max_spectrum"] is not None: ax_point_spec.plot( wavelength, arrays["point_max_spectrum"], linewidth=1.4, label=( "Strongest integrated-intensity pixel " f"(x={pmax['x']}, y={pmax['y']})" ), ) if arrays["point_min_spectrum"] is not None: ax_point_spec.plot( wavelength, arrays["point_min_spectrum"], linewidth=1.4, label=( "Lowest integrated-intensity pixel " f"(x={pmin['x']}, y={pmin['y']})" ), ) ax_point_spec.set_title( "This spectral window: average + spectra at lowest/highest " "integrated-intensity pixels" ) ax_point_spec.set_xlabel(wavelength_label) ax_point_spec.set_ylabel(str(hdu_unit_from_summary(summary))) ax_point_spec.legend(fontsize=8) ax_point_spec.grid(alpha=0.2) sat = summary["saturation_candidates"]["candidate_spatial_pixels"] nf = summary["sample_statistics"]["nonfinite_fraction"] fig.suptitle( f"{filename.name}\n" f"{summary['title']} — diagnostics\n" f"science y={y0}..{y1}; " f"non-finite={100 * nf:.4f}%; " f"saturation candidates={sat}", fontsize=13, ) safe_title = re.sub( r"[^A-Za-z0-9_.-]+", "_", summary["title"], ).strip("_") output_file = output_directory / ( f"{filename.stem}_HDU{summary['hdu_index']:02d}_" f"{safe_title}_QC.png" ) fig.savefig( output_file, dpi=PLOT_DPI, bbox_inches="tight", ) plt.close(fig) return output_file def hdu_unit_from_summary(summary): """Small plotting helper; spectra keep the native HDU radiance unit.""" # The unit itself is attached later by run_spice_qc when available. return summary.get("native_unit", "Radiance / intensity") def write_spice_qc_report(filename, summaries, output_directory): """Write one plain-text human-readable QC report.""" output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) txt_path = output_directory / f"{filename.stem}_QC.txt" lines = [] lines.append("=" * 78) lines.append("SPICE DATA DIAGNOSTICS REPORT") lines.append("=" * 78) lines.append(f"File: {filename.name}") lines.append("") lines.append("Science-region policy") lines.append("-" * 78) lines.append( f"Quantitative QC uses slit rows {Y_VALID_MIN} <= y <= " f"{Y_VALID_MAX}, clipped to each spectral window." ) lines.append( "Rows outside this interval remain visible in plots but are excluded " "from quantitative QC." ) lines.append("") for s in summaries: st = s["sample_statistics"] comp = s["spectrum_completeness"] ext = s["extrema"] pext = s["point_extrema"] sat = s["saturation_candidates"] anom = s["row_column_anomalies"] lines.append("=" * 78) lines.append(f"{s['title']}") lines.append("=" * 78) lines.append(f"HDU index: {s['hdu_index']}") lines.append(f"HDU name: {s['hdu_name']}") lines.append( "Cube shape [lambda, y, x]: " f"{s['cube_shape_lambda_y_x']}" ) lines.append( "Science Y interval actually used: " f"{s['science_y_min']}..{s['science_y_max']}" ) lines.append("") lines.append("DATA COMPLETENESS") lines.append("-" * 78) lines.append( "Non-finite samples: " f"{st['nonfinite_samples']} / " f"{st['total_samples_science_region']} " f"({100.0 * st['nonfinite_fraction']:.6f}%)" ) lines.append( f"NaN: {st['nan_samples']} " f"+Inf: {st['positive_inf_samples']} " f"-Inf: {st['negative_inf_samples']}" ) lines.append( f"Zero samples: {st['zero_samples']} " f"Negative samples: {st['negative_samples']}" ) lines.append( "Fully invalid spatial spectra: " f"{comp['fully_invalid_spectra']}" ) lines.append( "Partially invalid spatial spectra: " f"{comp['partially_invalid_spectra']}" ) lines.append("") lines.append("ROW / COLUMN INTEGRATED-INTENSITY EXTREMA") lines.append("-" * 78) lines.append( f"x_min={ext['x_min']} x_max={ext['x_max']} " f"y_min={ext['y_min']} y_max={ext['y_max']}" ) lines.append("") lines.append("SPATIAL-PIXEL INTEGRATED-INTENSITY EXTREMA") lines.append("-" * 78) pmin = pext["minimum"] pmax = pext["maximum"] if pmin["x"] is None: lines.append("Minimum point: unavailable") else: lines.append( "Minimum point: " f"x={pmin['x']}, y={pmin['y']}, " f"integrated intensity={pmin['integrated_intensity']:.8g}, " f"valid spectral fraction=" f"{100.0 * pmin['valid_spectral_fraction']:.2f}%" ) if pmax["x"] is None: lines.append("Maximum point: unavailable") else: lines.append( "Maximum point: " f"x={pmax['x']}, y={pmax['y']}, " f"integrated intensity={pmax['integrated_intensity']:.8g}, " f"valid spectral fraction=" f"{100.0 * pmax['valid_spectral_fraction']:.2f}%" ) lines.append( "The QC figure includes spectra at these two spatial pixels " "together with the average spectrum of the whole science region." ) lines.append("") lines.append("SATURATION") lines.append("-" * 78) lines.append( "Possible saturation spatial pixels: " f"{sat['candidate_spatial_pixels']}" ) lines.append( "Method: heuristic bright-tail + flat-top test; candidates are " "not instrument-confirmed saturation flags." ) lines.append("") lines.append("ROBUST ROW / COLUMN ANOMALIES") lines.append("-" * 78) lines.append( f"Anomalous X columns: {anom['anomalous_x_columns']}" ) lines.append( f"Anomalous Y rows: {anom['anomalous_y_rows']}" ) if anom["x_indices"]: lines.append(f"Anomalous X indices: {anom['x_indices']}") if anom["y_indices"]: lines.append(f"Anomalous Y indices: {anom['y_indices']}") lines.append("") txt_path.write_text( "\n".join(lines) + "\n", encoding="utf-8", ) return txt_path def run_spice_qc(filename, output_directory=QC_DIR): """Run all QC tools for one SPICE FITS file. Products -------- * one diagnostic PNG per observational spectral window * one plain-text QC report per FITS file """ filename = Path(filename) output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) summaries = [] plots = [] print() print(f"Diagnostics opening FITS: {filename}") with fits.open(filename, memmap=True) as hdul: for hdu_index, hdu in enumerate(hdul): if not is_observational_hdu(hdu): continue try: summary, arrays = analyze_spice_window_qc(hdu, hdu_index) summary["native_unit"] = str( hdu.header.get("BUNIT", "Radiance / intensity") ).strip() or "Radiance / intensity" plot_path = plot_spice_window_qc( filename, summary, arrays, output_directory ) summaries.append(summary) plots.append(plot_path) print( f" Diagnostics: {summary['title']} -> " f"{plot_path.name}" ) except Exception as exc: print( f" DIAGNOSTICS ERROR HDU {hdu_index} ({hdu.name}): {exc}" ) if not summaries: raise RuntimeError("No observational HDUs could be analyzed for diagnostics.") txt_path = write_spice_qc_report( filename, summaries, output_directory ) print(f"Diagnostics text report: {txt_path}") return { "plots": plots, "text_report": txt_path, } # ============================================================ # MAIN # ============================================================ def main(): DOWNLOAD_DIR.mkdir( parents=True, exist_ok=True, ) PLOT_DIR.mkdir( parents=True, exist_ok=True, ) QC_DIR.mkdir( parents=True, exist_ok=True, ) print("=" * 70) print( "Solar Orbiter / SPICE dual intensity pipeline" ) print("=" * 70) print(f"Start : {START_TIME}") print(f"End : {END_TIME}") print(f"Level : L{LEVEL}") print(f"Product : {PRODUCT}") print(f"Server : {SERVER_ROOT}") print(f"Log display : {USE_LOG_SCALE}") print( f"Percentiles : " f"{PERCENTILE_MIN} - {PERCENTILE_MAX}" ) print( f"Skip existing downloads: " f"{SKIP_EXISTING_FILES}" ) username, password = get_credentials() print() print("=" * 70) print("Searching / downloading SPICE data") print("=" * 70) files = download_spice_range( START_TIME, END_TIME, level=LEVEL, product=PRODUCT, output_directory=DOWNLOAD_DIR, username=username, password=password, ) if not files: print() print( "No files found in the requested interval." ) return print() print("=" * 70) print("Creating BOTH intensity versions") print("=" * 70) created_plots = [] created_qc = [] for file in files: print() print("-" * 70) print(f"Processing: {file.name}") try: simple_plot, physical_plot = ( plot_both_intensity_versions( file, output_directory=PLOT_DIR, use_log_scale=USE_LOG_SCALE, cmap=CMAP, pmin=PERCENTILE_MIN, pmax=PERCENTILE_MAX, ) ) if simple_plot is not None: created_plots.append(simple_plot) if physical_plot is not None: created_plots.append(physical_plot) except Exception as exc: print( f"ERROR processing {file.name}: {exc}" ) if RUN_QUALITY_CONTROL: try: qc_result = run_spice_qc( file, output_directory=QC_DIR, ) created_qc.append(qc_result) except Exception as exc: print( f"ERROR quality-control analysis for {file.name}: {exc}" ) print() print("=" * 70) print("Finished") print("=" * 70) print(f"Matching FITS files : {len(files)}") print(f"Plots created : {len(created_plots)}") print(f"FITS directory : {DOWNLOAD_DIR.resolve()}") print(f"Plot directory : {PLOT_DIR.resolve()}") print(f"QC directory : {QC_DIR.resolve()}") print(f"QC analyses created : {len(created_qc)}") print() print("For each FITS file you should normally get:") print( " *_simple_sum_log.png" ) print( " *_wavelength_integrated_log.png" ) if RUN_QUALITY_CONTROL: print(" *_HDU##_..._QC.png") print(" *_QC.txt") if __name__ == "__main__": main()