#!/usr/bin/env python3 """ Solar Orbiter / EUI Full Sun Imager (FSI) downloader + diagnostics. This script is designed for EUI/FSI Level-1 FITS images from the password-protected SIDC archive. It supports: - FSI 17.4 nm / 174 A files (also accepts "FSI171" as a user alias) - FSI 30.4 nm / 304 A files For every downloaded FITS image it creates: 1. Logarithmic image with the analysis-X interval marked. 2. Intensity histogram for the analysis region. 3. Intensity profile across the image centre in X. 4. Intensity profile across the image centre in Y. 5. X-row profile selected by the lowest summed intensity across X. 6. Y-column profile selected by the lowest summed intensity across Y. 7. X-row profile selected by the highest summed intensity across X. 8. Y-column profile selected by the highest summed intensity across Y. 9. TEMPINT, TEMP1DET, TEMP2DET temperatures plus DETGAINL, DETGAINH, GAINCOMB, READOUTM, DOWNLOAM and GAINTHRE from FITS headers. 10. Plain-text diagnostics report with statistics and selected rows/columns. Credentials ----------- Create pass_EUI.txt next to this script, for example: USERNAME = "abc" PASSWORD = "cde" The credentials file is parsed as text and is never executed. Existing non-empty FITS files are reused and are not downloaded again. FITS files are opened with memmap=False so Astropy can correctly handle BSCALE/BZERO/BLANK-scaled EUI Level-1 image data. 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 from astropy.io import fits # ============================================================ # USER SETTINGS # ============================================================ START_TIME = "2022-07-08 00:00:00" END_TIME = "2022-07-08 23:59:59" SERVER_ROOT = "https://www.sidc.be/EUI/data_internal/L1" # Accepted requested values: # "FSI171" -> matches FSI 17.4 nm / 174 A products # "FSI174" -> matches FSI 17.4 nm / 174 A products # "FSI304" -> matches FSI 30.4 nm / 304 A products REQUESTED_CHANNELS = ("FSI171", "FSI304") DOWNLOAD_DIR = Path("./eui_data") DIAGNOSTIC_DIR = Path("./eui_diagnostics") SKIP_EXISTING_FILES = True # Image display. USE_LOG_SCALE = True PERCENTILE_MIN = 0.5 PERCENTILE_MAX = 99.8 CMAP = "viridis" PLOT_DPI = 180 # ------------------------------------------------------------ # ANALYSIS REGION # ------------------------------------------------------------ # Quantitative diagnostics use only this X interval, inclusive. # The interval is clipped automatically to the actual image width. ANALYSIS_X_MIN = 0 ANALYSIS_X_MAX = 3000 # Zero/non-positive pixels are excluded from quantitative statistics. # This prevents padded/background zeros from forcing the minimum/median to 0. MIN_VALID_INTENSITY = 0.0 # Histogram. HIST_BINS = 150 HIST_LOG_X = True HIST_LOG_Y = True # Authentication. CREDENTIALS_FILE = Path("./pass_EUI.txt") READ_CREDENTIALS_FROM_FILE = True ASK_FOR_LOGIN_IF_FILE_MISSING = True USERNAME = None PASSWORD = None # ============================================================ # 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) # ============================================================ # AUTHENTICATION # ============================================================ def read_credentials_file(path): """Read USERNAME and PASSWORD from pass_EUI.txt without executing it.""" path = Path(path) 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") assignment_re = re.compile( r"^\s*(USERNAME|PASSWORD)\s*=\s*(['\"])(.*?)\2\s*(?:#.*)?$" ) values = {} 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 USERNAME = \"...\" and PASSWORD = \"...\"." ) return username, password, path def get_credentials(): username = USERNAME password = PASSWORD if READ_CREDENTIALS_FROM_FILE: try: username, password, path = read_credentials_file(CREDENTIALS_FILE) print(f"EUI credentials loaded from: {path}") print(f"Username: {username}") # Password is intentionally never printed. 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("EUI archive authentication") print("--------------------------") username = input("Username: ").strip() password = getpass.getpass("Password: ") return username, password def create_session(username=None, password=None): session = requests.Session() session.headers.update({ "User-Agent": ( "Mozilla/5.0 " "Solar-Orbiter-EUI-FSI-Diagnostics/1.0" ) }) if username and password: session.auth = (username, password) return session # ============================================================ # ARCHIVE SEARCH # ============================================================ def directory_url(day: datetime): return f"{SERVER_ROOT}/{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 EUI username/password." ) if response.status_code == 403: raise RuntimeError( "HTTP 403 Forbidden. The account may not have permission." ) 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", ".fit", ".fit.gz") ): urls.append(urljoin(url, href)) return sorted(set(urls)) def normalise_requested_channels(requested_channels): result = set() for channel in requested_channels: ch = str(channel).upper().replace(" ", "") if ch in {"FSI171", "FSI174", "171", "174"}: result.add("FSI174") elif ch in {"FSI304", "304"}: result.add("FSI304") else: raise ValueError( f"Unsupported channel {channel!r}. " "Use FSI171, FSI174, or FSI304." ) return result def filename_channel(filename): """Best-effort FSI channel classification from filename.""" name = filename.lower() # Common EUI descriptor-style strings. if "fsi304" in name or "fsi-304" in name: return "FSI304" if ( "fsi174" in name or "fsi-174" in name or "fsi171" in name or "fsi-171" in name ): return "FSI174" # More permissive fallback, but only if FSI occurs. if "fsi" in name: if "304" in name: return "FSI304" if "174" in name or "171" in name: return "FSI174" return None def extract_observation_time_from_filename(filename): """Best-effort UTC time extraction from EUI filenames.""" match = re.search( r"_(\d{8}T\d{6}(?:\d{3,6})?)", filename, ) if not match: return None stamp = match.group(1) for fmt in ( "%Y%m%dT%H%M%S%f", "%Y%m%dT%H%M%S", ): try: return datetime.strptime(stamp, fmt) except ValueError: pass return None def find_eui_files( start_time, end_time, session, requested_channels=REQUESTED_CHANNELS, ): start = parse_time(start_time) end = parse_time(end_time) if end < start: raise ValueError("END_TIME must be later than START_TIME.") wanted = normalise_requested_channels(requested_channels) matches = [] for day in iter_days(start, end): url = directory_url(day) try: files = get_directory_files(url, 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 channel = filename_channel(filename) if channel not in wanted: continue obs_time = extract_observation_time_from_filename(filename) # If the filename includes a time, enforce the exact requested # time range. If not, the day-directory match is retained. if obs_time is not None and not (start <= obs_time <= end): continue matches.append(file_url) return sorted(set(matches)) # ============================================================ # DOWNLOAD # ============================================================ def download_file( url, output_directory, session, skip_existing=True, ): """Download one archive file, or reuse an existing local copy. If a non-empty file with exactly the same archive filename already exists in ``output_directory``, it is returned immediately and no HTTP download request is made. A zero-byte local file is treated as incomplete and is downloaded again. A stale ``*.part`` file is removed before a new download starts. """ output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) destination = output_directory / Path(url).name temporary = destination.with_suffix(destination.suffix + ".part") # -------------------------------------------------------- # REUSE EXISTING DATA # -------------------------------------------------------- if destination.is_file(): size = destination.stat().st_size if skip_existing and size > 0: print( f"Data already exist -- reusing: {destination.name} " f"({size / 1024 / 1024:.2f} MB)" ) return destination if size == 0: print( f"Existing file is empty: {destination.name}; " "it will be downloaded again." ) # Remove a leftover incomplete transfer before retrying. if temporary.exists(): print(f"Removing incomplete download: {temporary.name}") temporary.unlink() 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 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() # Basic transfer-integrity check when Content-Length is available. if total and received != total: raise IOError( f"Incomplete download for {destination.name}: " f"received {received} bytes, expected {total}." ) temporary.replace(destination) except Exception: if temporary.exists(): temporary.unlink() raise print(f"Saved: {destination}") return destination def local_file_for_url(url, output_directory=DOWNLOAD_DIR): """Return the expected local Path for an archive URL.""" return Path(output_directory) / Path(url).name def local_data_exist(url, output_directory=DOWNLOAD_DIR): """True only when the matching local file exists and is non-empty.""" path = local_file_for_url(url, output_directory) return path.is_file() and path.stat().st_size > 0 def download_eui_range( start_time, end_time, requested_channels=REQUESTED_CHANNELS, output_directory=DOWNLOAD_DIR, username=None, password=None, ): """Find requested archive products and download only missing files.""" output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) session = create_session( username=username, password=password, ) urls = find_eui_files( start_time=start_time, end_time=end_time, session=session, requested_channels=requested_channels, ) print() print(f"Found {len(urls)} matching EUI/FSI archive file(s).") already_local = [ url for url in urls if local_data_exist(url, output_directory) ] missing = [ url for url in urls if not local_data_exist(url, output_directory) ] print(f"Already present locally : {len(already_local)}") print(f"Need downloading : {len(missing)}") paths = [] # First add files that are already present. They are still processed below. for url in already_local: path = local_file_for_url(url, output_directory) print(f"Reusing local data: {path.name}") paths.append(path) # Download only genuinely missing/empty files. for url in missing: 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}") # Stable archive order and no duplicates. return sorted(set(paths), key=lambda p: p.name) # ============================================================ # FITS IMAGE EXTRACTION # ============================================================ def reduce_to_2d_image(data): arr = np.asarray(data, dtype=np.float64) arr = np.squeeze(arr) while arr.ndim > 2: arr = np.nanmean(arr, axis=0) if arr.ndim != 2: raise RuntimeError( f"Could not reduce FITS data to 2-D; shape={arr.shape}" ) return arr def extract_best_eui_image(hdul): """Select the largest image-like HDU. This works for ordinary primary-image FITS products and for FITS files where Astropy exposes a compressed image extension as an image HDU. """ candidates = [] for index, hdu in enumerate(hdul): data = getattr(hdu, "data", None) if data is None or not isinstance(data, np.ndarray): continue if data.ndim < 2: continue try: image = reduce_to_2d_image(data) except Exception: continue if image.size == 0: continue score = int(image.shape[0] * image.shape[1]) candidates.append( (score, index, hdu, image) ) if not candidates: raise RuntimeError( "No usable 2-D image HDU was found in the FITS file." ) candidates.sort( key=lambda item: item[0], reverse=True, ) _, index, hdu, image = candidates[0] return index, hdu, image # ============================================================ # IMAGE DIAGNOSTICS # ============================================================ def finite_image_statistics(image): """Statistics from finite pixels strictly above MIN_VALID_INTENSITY.""" image = np.asarray(image, dtype=float) finite_mask = np.isfinite(image) valid_mask = finite_mask & (image > MIN_VALID_INTENSITY) valid = image[valid_mask] if valid.size == 0: raise ValueError( "No valid pixels above " f"MIN_VALID_INTENSITY={MIN_VALID_INTENSITY}." ) finite_count = int(np.count_nonzero(finite_mask)) nonfinite_count = int(image.size - finite_count) excluded_count = int(np.count_nonzero(finite_mask & ~valid_mask)) return { "finite_count": finite_count, "nonfinite_count": nonfinite_count, "nonfinite_fraction": float(nonfinite_count / image.size), "valid_count": int(valid.size), "excluded_nonpositive_count": excluded_count, "excluded_nonpositive_fraction": float(excluded_count / image.size), "minimum": float(np.min(valid)), "maximum": float(np.max(valid)), "mean": float(np.mean(valid)), "median": float(np.median(valid)), "std": float(np.std(valid)), } def analysis_x_bounds(nx): """Return inclusive analysis-X bounds clipped to the image width.""" x0 = max(0, int(ANALYSIS_X_MIN)) x1 = min(nx - 1, int(ANALYSIS_X_MAX)) if x0 > x1: raise ValueError( f"Invalid analysis X interval {ANALYSIS_X_MIN}.." f"{ANALYSIS_X_MAX} for image width nx={nx}." ) return x0, x1 def axis_sum_extrema(image, x0, x1): """Find low/high summed-intensity rows/columns using valid positive pixels.""" science = np.asarray(image[:, x0:x1 + 1], dtype=float) valid = np.isfinite(science) & (science > MIN_VALID_INTENSITY) row_valid_count = valid.sum(axis=1) col_valid_count = valid.sum(axis=0) row_sum = np.sum(np.where(valid, science, 0.0), axis=1, dtype=float) col_sum_local = np.sum(np.where(valid, science, 0.0), axis=0, dtype=float) row_sum = np.asarray(row_sum, dtype=float) col_sum_local = np.asarray(col_sum_local, dtype=float) # Entirely invalid/zero rows and columns are excluded. row_sum[row_valid_count == 0] = np.nan col_sum_local[col_valid_count == 0] = np.nan if not np.any(np.isfinite(row_sum)): raise ValueError("No valid Y-row sums in the analysis region.") if not np.any(np.isfinite(col_sum_local)): raise ValueError("No valid X-column sums in the analysis region.") y_min = int(np.nanargmin(row_sum)) y_max = int(np.nanargmax(row_sum)) x_min_local = int(np.nanargmin(col_sum_local)) x_max_local = int(np.nanargmax(col_sum_local)) x_min = int(x0 + x_min_local) x_max = int(x0 + x_max_local) return { "x0": int(x0), "x1": int(x1), "row_sum": row_sum, "col_sum": col_sum_local, "row_valid_count": row_valid_count, "col_valid_count": col_valid_count, "y_min": y_min, "y_max": y_max, "x_min": x_min, "x_max": x_max, "y_min_sum": float(row_sum[y_min]), "y_max_sum": float(row_sum[y_max]), "x_min_sum": float(col_sum_local[x_min_local]), "x_max_sum": float(col_sum_local[x_max_local]), } def center_profiles(image, x0, x1): """Return center profiles using the center of the analysis-X interval.""" ny, _ = image.shape x_center = int((x0 + x1) // 2) y_center = int(ny // 2) return { "x_center": x_center, "y_center": y_center, "x_coordinates": np.arange(x0, x1 + 1), "x_profile": image[y_center, x0:x1 + 1].astype( float, copy=True, ), "y_coordinates": np.arange(ny), "y_profile": image[:, x_center].astype( float, copy=True, ), } def profiles_from_axis_sum_extrema(image, extrema): """Return profiles selected from row/column summed-intensity extrema.""" x0 = extrema["x0"] x1 = extrema["x1"] y_min = extrema["y_min"] y_max = extrema["y_max"] x_min = extrema["x_min"] x_max = extrema["x_max"] return { "x_coordinates": np.arange(x0, x1 + 1), "y_coordinates": np.arange(image.shape[0]), # X-row profiles chosen from the minimum/maximum ROW SUMS. "x_row_low_sum": image[y_min, x0:x1 + 1].astype( float, copy=True, ), "x_row_high_sum": image[y_max, x0:x1 + 1].astype( float, copy=True, ), # Y-column profiles chosen from the minimum/maximum COLUMN SUMS. "y_col_low_sum": image[:, x_min].astype( float, copy=True, ), "y_col_high_sum": image[:, x_max].astype( float, copy=True, ), } def log_display_normalization(image): positive = image[ np.isfinite(image) & (image > 0) ] if positive.size == 0: raise ValueError( "Image contains no positive finite values for LogNorm." ) vmin, vmax = np.nanpercentile( positive, [PERCENTILE_MIN, PERCENTILE_MAX], ) 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, ) norm = LogNorm( vmin=float(vmin), vmax=float(vmax), clip=True, ) return masked, norm def infer_channel(filename, header): text = filename.lower() for key in ( "FILENAME", "FILTER", "WAVELNTH", "WAVELENGTH", "DETECTOR", "INSTRUME", "TELESCOP", ): value = header.get(key) if value is not None: text += " " + str(value).lower() if "304" in text or "30.4" in text: return "FSI 304" if ( "174" in text or "171" in text or "17.4" in text ): return "FSI 174" return "FSI" def native_unit(header): for key in ("BUNIT", "UNIT", "DATUNIT"): value = str( header.get(key, "") ).strip() if value: return value return "Intensity" def first_header_value(header, keys): """Return the first non-empty FITS header value among candidate keys.""" for key in keys: if key in header: value = header.get(key) if value not in (None, ""): return value, key return None, None def extract_temperature_and_gain(header): """Extract EUI detector temperatures and commanded gain/readout settings. Detector temperatures: TEMPINT - internal APS detector temperature [K] TEMP1DET - last measured APS detector temperature before date-average [K] TEMP2DET - earliest measured APS detector temperature after date-average [K] Gain/readout settings: DETGAINL - commanded low gain value DETGAINH - commanded high-gain value GAINCOMB - commanded low/high gain combination READOUTM - commanded FEE readout mode DOWNLOAM - commanded FEE download mode GAINTHRE - commanded threshold value for H/L gain """ temperature_definitions = ( ("TEMPINT", "Internal APS detector temperature"), ("TEMP1DET", "Last measured APS detector temperature before date-average"), ("TEMP2DET", "Earliest measured APS detector temperature after date-average"), ) temperatures = [] for key, description in temperature_definitions: temperatures.append({ "key": key, "value": header.get(key), "unit": "K", "description": description, }) gain_definitions = ( ("DETGAINL", "Commanded low gain value"), ("DETGAINH", "Commanded high-gain value"), ("GAINCOMB", "Commanded low/high gain combination"), ("READOUTM", "Commanded FEE readout mode"), ("DOWNLOAM", "Commanded FEE download mode"), ("GAINTHRE", "Commanded threshold value for H/L gain"), ) gain_parameters = [] for key, description in gain_definitions: gain_parameters.append({ "key": key, "value": header.get(key), "description": description, }) return { "temperatures": temperatures, "gain_parameters": gain_parameters, } def format_temperature_entry(entry): """Format one EUI detector-temperature entry for plot/report output.""" value = entry["value"] if value is None: value_text = "N/A" else: try: value_text = f"{float(value):.3f} K" except (TypeError, ValueError): value_text = f"{value} K" return ( f"{entry['key']}={value_text} " f"({entry['description']})" ) def format_gain_entry(entry): """Format one EUI/FEE gain or readout parameter.""" value = entry["value"] value_text = "N/A" if value is None else str(value) return ( f"{entry['key']}={value_text} " f"({entry['description']})" ) def format_header_quantity(value, key, default_label): if value is None: return f"{default_label}=N/A" return ( f"{default_label}={value}" + (f" [{key}]" if key else "") ) def analyze_eui_file(filename): filename = Path(filename) # EUI L1 images may use BSCALE/BZERO/BLANK. memmap=False is required # so Astropy can correctly apply FITS image scaling. with fits.open( filename, memmap=False, do_not_scale_image_data=False, uint=True, ) as hdul: hdu_index, hdu, image = extract_best_eui_image(hdul) header = hdu.header.copy() ny, nx = image.shape x0, x1 = analysis_x_bounds(nx) # All quantitative image statistics are restricted to X=0..3000 # (clipped to the actual image width). science_image = image[:, x0:x1 + 1] stats = finite_image_statistics(science_image) extrema = axis_sum_extrema(image, x0, x1) centers = center_profiles(image, x0, x1) selected_profiles = profiles_from_axis_sum_extrema( image, extrema, ) housekeeping = extract_temperature_and_gain(header) return { "filename": filename, "image": image, "science_image": science_image, "header": header, "hdu_index": int(hdu_index), "stats": stats, "extrema": extrema, "centers": centers, "selected_profiles": selected_profiles, "channel": infer_channel(filename.name, header), "unit": native_unit(header), "housekeeping": housekeeping, "analysis_x_min": int(x0), "analysis_x_max": int(x1), } # ============================================================ # DIAGNOSTIC PLOT # ============================================================ def plot_eui_diagnostics( result, output_directory=DIAGNOSTIC_DIR, ): output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) filename = result["filename"] image = result["image"] science_image = result["science_image"] stats = result["stats"] extrema = result["extrema"] centers = result["centers"] profiles = result["selected_profiles"] channel = result["channel"] unit = result["unit"] housekeeping = result["housekeeping"] x0 = result["analysis_x_min"] x1 = result["analysis_x_max"] fig = plt.figure(figsize=(18, 20), constrained_layout=True) gs = fig.add_gridspec( 5, 2, height_ratios=[0.72, 1.25, 0.9, 0.9, 0.9], ) ax_table = fig.add_subplot(gs[0, :]) ax_image = fig.add_subplot(gs[1, 0]) ax_hist = fig.add_subplot(gs[1, 1]) ax_center_x = fig.add_subplot(gs[2, 0]) ax_center_y = fig.add_subplot(gs[2, 1]) ax_low_x = fig.add_subplot(gs[3, 0]) ax_low_y = fig.add_subplot(gs[3, 1]) ax_high_x = fig.add_subplot(gs[4, 0]) ax_high_y = fig.add_subplot(gs[4, 1]) # ======================================================== # TECHNICAL PARAMETERS TABLE # ======================================================== ax_table.axis("off") temperature_map = { entry["key"]: entry["value"] for entry in housekeeping["temperatures"] } gain_map = { entry["key"]: entry["value"] for entry in housekeeping["gain_parameters"] } def fmt_number(value, digits=6): if value is None: return "N/A" try: return f"{float(value):.{digits}g}" except (TypeError, ValueError): return str(value) def fmt_temperature(value): if value is None: return "N/A" try: return f"{float(value):.3f} K" except (TypeError, ValueError): return f"{value} K" technical_rows = [ ["Analysis X range", f"{x0} .. {x1} px"], ["Mean intensity", fmt_number(stats["mean"])], ["Standard deviation", fmt_number(stats["std"])], ["Minimum intensity", fmt_number(stats["minimum"])], ["Maximum intensity", fmt_number(stats["maximum"])], ["TEMPINT", fmt_temperature(temperature_map.get("TEMPINT"))], ["TEMP1DET", fmt_temperature(temperature_map.get("TEMP1DET"))], ["TEMP2DET", fmt_temperature(temperature_map.get("TEMP2DET"))], ["DETGAINL", fmt_number(gain_map.get("DETGAINL"))], ["DETGAINH", fmt_number(gain_map.get("DETGAINH"))], ["GAINCOMB", "N/A" if gain_map.get("GAINCOMB") is None else str(gain_map.get("GAINCOMB"))], ["READOUTM", "N/A" if gain_map.get("READOUTM") is None else str(gain_map.get("READOUTM"))], ["DOWNLOAM", "N/A" if gain_map.get("DOWNLOAM") is None else str(gain_map.get("DOWNLOAM"))], ["GAINTHRE", "N/A" if gain_map.get("GAINTHRE") is None else str(gain_map.get("GAINTHRE"))], ] left_rows = technical_rows[:8] right_rows = technical_rows[8:] n_rows = max(len(left_rows), len(right_rows)) while len(left_rows) < n_rows: left_rows.append(["", ""]) while len(right_rows) < n_rows: right_rows.append(["", ""]) table_rows = [left_rows[i] + right_rows[i] for i in range(n_rows)] table = ax_table.table( cellText=table_rows, colLabels=["Parameter", "Value", "Parameter", "Value"], cellLoc="right", colLoc="right", loc="center", colWidths=[0.24, 0.22, 0.24, 0.22], ) table.auto_set_font_size(False) table.set_fontsize(9) table.scale(1.0, 1.22) for cell in table.get_celld().values(): cell.get_text().set_ha("right") cell.get_text().set_va("center") ax_table.set_title("Technical parameters", loc="right", fontsize=11, pad=3) # ======================================================== # LOG-SCALED IMAGE # ======================================================== if USE_LOG_SCALE: try: image_for_plot, norm = log_display_normalization(image) im = ax_image.imshow(image_for_plot, origin="lower", aspect="equal", cmap=CMAP, norm=norm, interpolation="nearest") except ValueError: im = ax_image.imshow(image, origin="lower", aspect="equal", cmap=CMAP, interpolation="nearest") else: im = ax_image.imshow(image, origin="lower", aspect="equal", cmap=CMAP, interpolation="nearest") fig.colorbar(im, ax=ax_image, pad=0.02, label=unit) ax_image.axvline(x0, linestyle="--", linewidth=2.4, color="white", alpha=0.8) ax_image.axvline(x1, linestyle="--", linewidth=2.4, color="white", alpha=0.8) ax_image.axhline(extrema["y_min"], linestyle=":", linewidth=3.3, color="cyan", label=f"lowest row sum: y={extrema['y_min']}") ax_image.axhline(extrema["y_max"], linestyle=":", linewidth=3.3, color="magenta", label=f"highest row sum: y={extrema['y_max']}") ax_image.axvline(extrema["x_min"], linestyle=":", linewidth=3.3, color="cyan", label=f"lowest column sum: x={extrema['x_min']}") ax_image.axvline(extrema["x_max"], linestyle=":", linewidth=3.3, color="magenta", label=f"highest column sum: x={extrema['x_max']}") ax_image.set_title(f"{channel} image — logarithmic display\nanalysis region: {x0} <= X <= {x1}") ax_image.set_xlabel("X [pixel]") ax_image.set_ylabel("Y [pixel]") ax_image.legend(fontsize=7, loc="upper right") finite = science_image[np.isfinite(science_image)] hist_values = finite[finite > MIN_VALID_INTENSITY] if hist_values.size: if HIST_LOG_X and np.nanmin(hist_values) > 0 and np.nanmax(hist_values) > np.nanmin(hist_values): bins = np.geomspace(np.nanmin(hist_values), np.nanmax(hist_values), HIST_BINS) else: bins = HIST_BINS ax_hist.hist(hist_values, bins=bins) if HIST_LOG_X and hist_values.size and np.any(hist_values > 0): ax_hist.set_xscale("log") if HIST_LOG_Y: ax_hist.set_yscale("log") ax_hist.axvline(stats["mean"], linestyle="--", linewidth=3, label=f"mean={stats['mean']:.5g}") ax_hist.axvline(stats["median"], linestyle=":", linewidth=3.6, label=f"median={stats['median']:.5g}") ax_hist.set_title(f"Intensity histogram, X={x0}..{x1}") ax_hist.set_xlabel(unit) ax_hist.set_ylabel("Number of pixels") ax_hist.legend(fontsize=8) ax_center_x.plot(centers["x_coordinates"], centers["x_profile"], linewidth=3) ax_center_x.set_title(f"Centre X-row profile (y={centers['y_center']})") ax_center_x.set_xlabel("X [pixel]") ax_center_x.set_ylabel(unit) ax_center_x.grid(alpha=0.2) ax_center_y.plot(centers["y_coordinates"], centers["y_profile"], linewidth=3) ax_center_y.set_title(f"Centre Y-column profile (x={centers['x_center']})") ax_center_y.set_xlabel("Y [pixel]") ax_center_y.set_ylabel(unit) ax_center_y.grid(alpha=0.2) ax_low_x.plot(profiles["x_coordinates"], profiles["x_row_low_sum"], linewidth=3) ax_low_x.set_title(f"X-row with LOWEST summed intensity (y={extrema['y_min']}, sum={extrema['y_min_sum']:.6g})") ax_low_x.set_xlabel("X [pixel]") ax_low_x.set_ylabel(unit) ax_low_x.grid(alpha=0.2) ax_low_y.plot(profiles["y_coordinates"], profiles["y_col_low_sum"], linewidth=3) ax_low_y.set_title(f"Y-column with LOWEST summed intensity (x={extrema['x_min']}, sum={extrema['x_min_sum']:.6g})") ax_low_y.set_xlabel("Y [pixel]") ax_low_y.set_ylabel(unit) ax_low_y.grid(alpha=0.2) ax_high_x.plot(profiles["x_coordinates"], profiles["x_row_high_sum"], linewidth=3) ax_high_x.set_title(f"X-row with HIGHEST summed intensity (y={extrema['y_max']}, sum={extrema['y_max_sum']:.6g})") ax_high_x.set_xlabel("X [pixel]") ax_high_x.set_ylabel(unit) ax_high_x.grid(alpha=0.2) ax_high_y.plot(profiles["y_coordinates"], profiles["y_col_high_sum"], linewidth=3) ax_high_y.set_title(f"Y-column with HIGHEST summed intensity (x={extrema['x_max']}, sum={extrema['x_max_sum']:.6g})") ax_high_y.set_xlabel("Y [pixel]") ax_high_y.set_ylabel(unit) ax_high_y.grid(alpha=0.2) fig.suptitle(f"{filename.name}\n{channel} diagnostics", fontsize=13) safe_channel = re.sub(r"[^A-Za-z0-9_.-]+", "_", channel).strip("_") output_file = output_directory / f"{filename.stem}_{safe_channel}_diagnostics.png" fig.savefig(output_file, dpi=PLOT_DPI, bbox_inches="tight") plt.close(fig) return output_file # ============================================================ # TEXT REPORT # ============================================================ def write_eui_diagnostics_report( result, output_directory=DIAGNOSTIC_DIR, ): output_directory = Path(output_directory) output_directory.mkdir(parents=True, exist_ok=True) filename = result["filename"] image = result["image"] stats = result["stats"] extrema = result["extrema"] centers = result["centers"] channel = result["channel"] unit = result["unit"] housekeeping = result["housekeeping"] x0 = result["analysis_x_min"] x1 = result["analysis_x_max"] lines = [] lines.append("=" * 78) lines.append( "SOLAR ORBITER EUI / FSI IMAGE DIAGNOSTICS REPORT" ) lines.append("=" * 78) lines.append(f"File: {filename.name}") lines.append(f"Channel: {channel}") lines.append(f"Image HDU index: {result['hdu_index']}") lines.append(f"Image shape [y, x]: {list(image.shape)}") lines.append(f"Intensity unit: {unit}") lines.append(f"Analysis X interval: {x0}..{x1}") lines.append("") lines.append("INSTRUMENT / HOUSEKEEPING") lines.append("-" * 78) lines.append("APS DETECTOR TEMPERATURES") for entry in housekeeping["temperatures"]: lines.append(format_temperature_entry(entry)) lines.append("") lines.append("GAIN / FEE SETTINGS") for entry in housekeeping["gain_parameters"]: lines.append(format_gain_entry(entry)) lines.append("") lines.append("IMAGE STATISTICS — ANALYSIS REGION ONLY") lines.append("-" * 78) lines.append(f"Mean: {stats['mean']:.10g}") lines.append(f"Standard deviation: {stats['std']:.10g}") lines.append(f"Minimum: {stats['minimum']:.10g}") lines.append(f"Maximum: {stats['maximum']:.10g}") lines.append("") lines.append("SUMMED-INTENSITY ROW / COLUMN SELECTIONS") lines.append("-" * 78) lines.append( f"Lowest X-row sum: y={extrema['y_min']}, " f"sum={extrema['y_min_sum']:.10g}" ) lines.append( f"Highest X-row sum: y={extrema['y_max']}, " f"sum={extrema['y_max_sum']:.10g}" ) lines.append( f"Lowest Y-column sum: x={extrema['x_min']}, " f"sum={extrema['x_min_sum']:.10g}" ) lines.append( f"Highest Y-column sum: x={extrema['x_max']}, " f"sum={extrema['x_max_sum']:.10g}" ) lines.append("") lines.append("CENTER PROFILE POSITIONS") lines.append("-" * 78) lines.append(f"Centre X-row: y={centers['y_center']}") lines.append(f"Centre Y-column: x={centers['x_center']}") lines.append("") output_file = output_directory / ( f"{filename.stem}_diagnostics.txt" ) output_file.write_text( "\n".join(lines) + "\n", encoding="utf-8", ) return output_file # ============================================================ # PROCESS ONE FILE # ============================================================ def process_eui_file( filename, output_directory=DIAGNOSTIC_DIR, ): filename = Path(filename) print() print("-" * 78) print(f"Processing: {filename.name}") result = analyze_eui_file(filename) plot_path = plot_eui_diagnostics( result, output_directory=output_directory, ) report_path = write_eui_diagnostics_report( result, output_directory=output_directory, ) print(f"Diagnostic plot: {plot_path}") print(f"Text report : {report_path}") return { "plot": plot_path, "text_report": report_path, } # ============================================================ # MAIN # ============================================================ def main(): DOWNLOAD_DIR.mkdir( parents=True, exist_ok=True, ) DIAGNOSTIC_DIR.mkdir( parents=True, exist_ok=True, ) print("=" * 78) print( "Solar Orbiter / EUI FSI downloader + diagnostics" ) print("=" * 78) print(f"Start : {START_TIME}") print(f"End : {END_TIME}") print(f"Server : {SERVER_ROOT}") print(f"Channels : {REQUESTED_CHANNELS}") print(f"Log display: {USE_LOG_SCALE}") username, password = get_credentials() files = download_eui_range( START_TIME, END_TIME, requested_channels=REQUESTED_CHANNELS, output_directory=DOWNLOAD_DIR, username=username, password=password, ) if not files: print() print( "No matching EUI/FSI FITS files were found." ) return products = [] for file in files: try: products.append( process_eui_file( file, output_directory=DIAGNOSTIC_DIR, ) ) except Exception as exc: print( f"ERROR processing {file.name}: {exc}" ) print() print("=" * 78) print("Finished") print("=" * 78) print( f"Matched/downloaded files: {len(files)}" ) print( f"Processed files : {len(products)}" ) print( f"FITS directory : " f"{DOWNLOAD_DIR.resolve()}" ) print( f"Diagnostics directory : " f"{DIAGNOSTIC_DIR.resolve()}" ) if __name__ == "__main__": main()