diff --git a/.gitignore b/.gitignore index b71134a..4233169 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,10 @@ __pycache__/ *.pyc *.swp +*.ipynb +**/.ipynb_checkpoints/ +.ipynb_checkpoints/ data/ output/galleries/ scripts/*/output/ diff --git a/README.md b/README.md index e3feeb5..0cd360b 100644 --- a/README.md +++ b/README.md @@ -1,28 +1,125 @@ # HamSCI Amateur Radio Communications (HARC) Database Plotting Toolkit +Tools for loading and plotting HF amateur radio propagation spot data — from +RBN, WSPRNet, and PSKReporter — as distributed via the MIT Haystack CEDAR +[Madrigal database](http://cedar.openmadrigal.org/). + ## Requirements -* Ubuntu 18.04 LTS -* Anaconda Python 3.7 -## Install Instructions -``` -conda install -c conda-forge cartopy +* Linux or macOS (tested on Ubuntu 22.04+). Windows is not supported directly; + WSL2 works. +* Python 3.10 or newer. +* Conda recommended (miniconda or Anaconda) for managing the scientific Python + stack. + +Python dependencies are listed in `requirements.txt`. The main ones: +`numpy`, `pandas`, `xarray`, `netCDF4`, `tables` (pytables), `h5py`, +`matplotlib`, `cartopy`, `madrigalWeb`. + +## Install + +```bash +# Create and activate a conda environment +conda create -n harc_plot -c conda-forge python=3.11 cartopy pytables netCDF4 xarray pandas matplotlib tqdm +conda activate harc_plot + +# Clone and install git clone https://github.com/HamSCI/harc_plot.git cd harc_plot pip install -e . ``` -## Test Installation -From the harc_plot directory: -``` +`cartopy` is installed via conda because its binary dependencies (PROJ, GEOS) +are awkward to get from pip. + +## Data input + +The toolkit reads spot data in two formats: + +1. **Madrigal HDF5 (preferred)** — files named `rsdYYYY-MM-DD.01.hdf5`, placed + in `data/madrigal/`. Downloadable from the Madrigal database via + `madrigalWeb`. A sample (`data/madrigal/rsd2015-06-10.01.hdf5`, ~39 MB, + covering all three networks) ships with the repository for testing. +2. **Legacy CSV** — files named `YYYY-MM-DD.csv.bz2` placed in `data/spot_csvs/` + (pre-Madrigal HamSCI pipeline). + +`harc_plot.gen_lib.load_spots_csv()` auto-detects which format is present. + +### A note about Madrigal timestamps + +A subset of Madrigal ham-radio files (empirically: most of 2023 and 2024, +other years clean) have their `ut1_unix`/`ut2_unix` columns shifted by a +local-timezone offset. The HDF5 loader (`load_madrigal_hdf5()`) detects this +per file by comparing `ut1_unix` against the authoritative `y/m/d/h/m/s` byte +columns and corrects or falls back as needed. You will see a `RuntimeWarning` +on any affected file. + +## Test installation + +From the `harc_plot` directory: + +```bash cd scripts/basic_histogram ./histograms.py ``` -The following plots will be generated: + +This loads the bundled sample (`data/madrigal/rsd2015-06-10.01.hdf5`), bins it +into 2D histograms, and produces four plots: + ``` -output/galleries/histograms/World/ut_hrs/20170701.0000UT-20170702.0000UT.data.nc.ut_hrs.spot_density.png -output/galleries/histograms/World/slt_mid/20170701.0000UT-20170702.0000UT.data.nc.slt_mid.spot_density.png -output/galleries/histograms/World/ut_hrs/dailies/20170701.0000UT-20170702.0000UT.data.nc.ut_hrs.spot_density.png -output/galleries/histograms/World/slt_mid/dailies/20170701.0000UT-20170702.0000UT.data.nc.slt_mid.spot_density.png +output/galleries/histograms/World/ut_hrs/20150610.0000UT-20150611.0000UT.data.nc.ut_hrs.spot_density.png +output/galleries/histograms/World/slt_mid/20150610.0000UT-20150611.0000UT.data.nc.slt_mid.spot_density.png +output/galleries/histograms/World/ut_hrs/dailies/20150610.0000UT-20150611.0000UT.data.nc.ut_hrs.spot_density.png +output/galleries/histograms/World/slt_mid/dailies/20150610.0000UT-20150611.0000UT.data.nc.slt_mid.spot_density.png ``` -![Example Plot](example_plots/20170701.0000UT-20170702.0000UT.data.nc.ut_hrs.spot_density.png) + +![Example Plot](example_plots/20150610.0000UT-20150611.0000UT.data.nc.ut_hrs.spot_density.png) + +## Data sources and acknowledgements + +### Solar F10.7 cm radio flux + +Solar F10.7 cm radio flux values bundled in `data/omni/` (and re-downloaded by +`data/omni/update_omni_data.py`) are obtained from the NASA OMNIWeb service +() as part of the OMNI2 hourly compiled dataset, +which ingests F10.7 observations made at the Dominion Radio Astrophysical +Observatory (DRAO), Penticton, British Columbia, Canada, by Natural Resources +Canada. + +**Suggested acknowledgement for publications using this toolkit's F10.7 data:** + +> Solar F10.7 cm radio flux data were obtained from the NASA OMNIWeb service +> (https://omniweb.gsfc.nasa.gov) as part of the OMNI2 hourly compiled dataset +> (Papitashvili & King, 2020), which ingests F10.7 observations made at the +> Dominion Radio Astrophysical Observatory, Penticton, British Columbia, +> Canada, by Natural Resources Canada. + +Note that F10.7 is determined **once per day**, so it cannot resolve solar +flares. Use GOES X-ray irradiance for flare timing (see below). + +### GOES X-ray irradiance + +`harc_plot.goes.read_goes_ncei()` retrieves 1-minute GOES XRS data from the NOAA +National Centers for Environmental Information (NCEI) over HTTPS, caching files in +`data/goes/`. Two archives are used and normalized to a common set of column names +(`A_AVG` for 0.05–0.4 nm, `B_AVG` for 0.1–0.8 nm): + +* GOES 8–15 — GOES Space Environment Monitor archive, one netCDF per month +* GOES 16–19 — GOES-R series L2 `xrsf-l2-avg1m` product, one netCDF per day + +The older `read_goes()` targets `ftp://satdat.ngdc.noaa.gov`, which NOAA retired +along with anonymous FTP; it now always fails and returns `None`. Prefer +`read_goes_ncei()`. + +**Flare magnitudes are not comparable across satellite generations.** The pre-R +and GOES-R XRS instruments report on different flux scales — for the 2017-09-06 +flare, GOES-15 gives X9.3 (the accepted classification), GOES-13 X10.4, and +GOES-16 X5.2. `goes_xrs_candidates()` therefore prefers the *primary operational* +satellite for each epoch, which reproduces the standard classifications. + +**Suggested acknowledgement for publications using this toolkit's GOES data:** + +> GOES X-ray irradiance data were obtained from the NOAA National Centers for +> Environmental Information (https://www.ncei.noaa.gov), from the GOES Space +> Environment Monitor archive (GOES 8–15) and the GOES-R series Level 2 X-Ray +> Sensor 1-minute average product (GOES 16–19). diff --git a/data/histograms/World/20150610.data.nc.bz2 b/data/histograms/World/20150610.data.nc.bz2 new file mode 100644 index 0000000..d2bad66 Binary files /dev/null and b/data/histograms/World/20150610.data.nc.bz2 differ diff --git a/data/histograms/World/20170701.data.nc.bz2 b/data/histograms/World/20170701.data.nc.bz2 deleted file mode 100644 index 4b3822a..0000000 Binary files a/data/histograms/World/20170701.data.nc.bz2 and /dev/null differ diff --git a/data/spot_csvs/2017-07-01.csv.bz2 b/data/madrigal/rsd2015-06-10.01.hdf5 similarity index 51% rename from data/spot_csvs/2017-07-01.csv.bz2 rename to data/madrigal/rsd2015-06-10.01.hdf5 index 0604e4a..cb707dd 100644 Binary files a/data/spot_csvs/2017-07-01.csv.bz2 and b/data/madrigal/rsd2015-06-10.01.hdf5 differ diff --git a/data/omni/omni_data_reduced.txt.bz2 b/data/omni/omni_data_reduced.txt.bz2 index b9f2e85..e15e390 100644 Binary files a/data/omni/omni_data_reduced.txt.bz2 and b/data/omni/omni_data_reduced.txt.bz2 differ diff --git a/data/omni/update_omni_data.py b/data/omni/update_omni_data.py index 82af5ee..cfdc32d 100755 --- a/data/omni/update_omni_data.py +++ b/data/omni/update_omni_data.py @@ -1,4 +1,25 @@ #!/usr/bin/env python3 +""" +Download solar / geomagnetic indices from NASA OMNIWeb. + +The bundled `omni_data_reduced.txt.bz2` in this directory contains hourly +values of Kp, R (sunspot number), Dst, ap, **F10.7**, AE, AL, and AU from +2000-01-01 onward, used throughout `harc_plot` for solar-cycle overlays and +SW-activity filtering. + +Data source and suggested acknowledgement: + + Solar F10.7 cm radio flux data were obtained from the NASA OMNIWeb + service (https://omniweb.gsfc.nasa.gov) as part of the OMNI2 hourly + compiled dataset (Papitashvili & King, 2020), which ingests F10.7 + observations made at the Dominion Radio Astrophysical Observatory, + Penticton, British Columbia, Canada, by Natural Resources Canada. + +The Kp, Dst, ap, AE/AL/AU indices in the same file originate with the +International Service of Geomagnetic Indices (ISGI) / GFZ Potsdam / Kyoto +WDC; consult the OMNI2 documentation for attribution details before using +them. +""" import os import datetime from collections import OrderedDict diff --git a/example_plots/20150610.0000UT-20150611.0000UT.data.nc.slt_mid.spot_density.png b/example_plots/20150610.0000UT-20150611.0000UT.data.nc.slt_mid.spot_density.png new file mode 100644 index 0000000..9de3cb9 Binary files /dev/null and b/example_plots/20150610.0000UT-20150611.0000UT.data.nc.slt_mid.spot_density.png differ diff --git a/example_plots/20150610.0000UT-20150611.0000UT.data.nc.ut_hrs.spot_density.png b/example_plots/20150610.0000UT-20150611.0000UT.data.nc.ut_hrs.spot_density.png new file mode 100644 index 0000000..63e833c Binary files /dev/null and b/example_plots/20150610.0000UT-20150611.0000UT.data.nc.ut_hrs.spot_density.png differ diff --git a/example_plots/20170701.0000UT-20170702.0000UT.data.nc.slt_mid.spot_density.png b/example_plots/20170701.0000UT-20170702.0000UT.data.nc.slt_mid.spot_density.png deleted file mode 100644 index f1d6242..0000000 Binary files a/example_plots/20170701.0000UT-20170702.0000UT.data.nc.slt_mid.spot_density.png and /dev/null differ diff --git a/example_plots/20170701.0000UT-20170702.0000UT.data.nc.ut_hrs.spot_density.png b/example_plots/20170701.0000UT-20170702.0000UT.data.nc.ut_hrs.spot_density.png deleted file mode 100644 index cf86a7d..0000000 Binary files a/example_plots/20170701.0000UT-20170702.0000UT.data.nc.ut_hrs.spot_density.png and /dev/null differ diff --git a/harc_plot/gen_lib.py b/harc_plot/gen_lib.py index cb68f03..b74b1ad 100644 --- a/harc_plot/gen_lib.py +++ b/harc_plot/gen_lib.py @@ -11,16 +11,19 @@ import numpy as np import pandas as pd +import tables +import warnings + import pickle import bz2 from . import geopack from . import calcSun -de_prop = {'marker':'^','edgecolor':'k','facecolor':'white'} -dxf_prop = {'marker':'*','color':'blue'} -dxf_leg_size = 150 -dxf_plot_size = 50 +de_prop = {'marker':'^','edgecolor':'k','facecolor':'white'} +dxf_prop = {'marker':'*','color':'blue'} +dxf_leg_size = 150 +dxf_plot_size = 50 rcp = matplotlib.rcParams rcp['figure.titlesize'] = 'xx-large' @@ -289,29 +292,30 @@ def __cb_safe(self): self.cb_safe_dct = fc def cdict_to_cmap(cdict,name='CustomCMAP',vmin=0.,vmax=30.): - norm = matplotlib.colors.Normalize(vmin=vmin,vmax=vmax) - - red = [] - green = [] - blue = [] - - keys = list(cdict.keys()) - keys.sort() - - for x in keys: - r,g,b, = cdict[x] - x = norm(x) - r = r/255. - g = g/255. - b = b/255. - red.append( (x, r, r)) - green.append( (x, g, g)) - blue.append( (x, b, b)) - cdict = {'red' : tuple(red), - 'green' : tuple(green), - 'blue' : tuple(blue)} - cmap = matplotlib.colors.LinearSegmentedColormap(name, cdict) - return cmap + norm = matplotlib.colors.Normalize(vmin=vmin,vmax=vmax) + + red = [] + green = [] + blue = [] + + keys = list(cdict.keys()) + keys.sort() + + for x in keys: + r,g,b, = cdict[x] + x = norm(x) + r = r/255. + g = g/255. + b = b/255. + red.append( (x, r, r)) + green.append( (x, g, g)) + blue.append( (x, b, b)) + + cdict = {'red' : tuple(red), + 'green' : tuple(green), + 'blue' : tuple(blue)} + cmap = matplotlib.colors.LinearSegmentedColormap(name, cdict) + return cmap def sun_pos(dt=None): """This function computes a rough estimate of the coordinates for @@ -445,30 +449,184 @@ def adjust_axes(ax_0,ax_1): ax_0_pos[2] = ax_1_pos[2] ax_0.set_position(ax_0_pos) -def regional_filter(region,df,kind='mids'): - rgnd = regions[region] +def regional_filter(region, df, kind='mids'): + rgnd = regions[region] lat_lim = rgnd['lat_lim'] lon_lim = rgnd['lon_lim'] if kind == 'mids': - tf_md_lat = np.logical_and(df.md_lat >= lat_lim[0], df.md_lat < lat_lim[1]) - tf_md_long = np.logical_and(df.md_long >= lon_lim[0], df.md_long < lon_lim[1]) - tf_0 = np.logical_and(tf_md_lat,tf_md_long) - tf = tf_0 - df = df[tf].copy() + tf = ((df['md_lat'] >= lat_lim[0]) & (df['md_lat'] < lat_lim[1]) & + (df['md_long'] >= lon_lim[0]) & (df['md_long'] < lon_lim[1])) elif kind == 'endpoints': - tf_rx_lat = np.logical_and(df.rx_lat >= lat_lim[0], df.rx_lat < lat_lim[1]) - tf_rx_long = np.logical_and(df.rx_long >= lon_lim[0], df.rx_long < lon_lim[1]) - tf_rx = np.logical_and(tf_rx_lat,tf_rx_long) - - tf_tx_lat = np.logical_and(df.tx_lat >= lat_lim[0], df.tx_lat < lat_lim[1]) - tf_tx_long = np.logical_and(df.tx_long >= lon_lim[0], df.tx_long < lon_lim[1]) - tf_tx = np.logical_and(tf_tx_lat,tf_tx_long) - tf = np.logical_or(tf_rx,tf_tx) + tf_rx = ((df['rx_lat'] >= lat_lim[0]) & (df['rx_lat'] < lat_lim[1]) & + (df['rx_long'] >= lon_lim[0]) & (df['rx_long'] < lon_lim[1])) + tf_tx = ((df['tx_lat'] >= lat_lim[0]) & (df['tx_lat'] < lat_lim[1]) & + (df['tx_long'] >= lon_lim[0]) & (df['tx_long'] < lon_lim[1])) + tf = tf_rx | tf_tx + else: + raise ValueError(f"regional_filter kind must be 'mids' or 'endpoints', got {kind!r}") + + return df[tf].copy() + + +# Madrigal ssrc 3-letter code -> harc_plot integer source code (keep in sync with `sources` dict above). +_MADRIGAL_SSRC_MAP = {'WSP': 1, 'RBN': 2, 'PSK': 3} + +# Band plan: (low_MHz, high_MHz, band_meters). The 'band' column stores wavelength +# in meters to match the existing harc_plot convention (BandData's 'meters' field and +# calculate_histograms' `df["band"] == band["meters"]` filter). Includes non-ham / +# experimental frequencies observed in Madrigal data so nothing is silently dropped; +# BandData decides what actually gets plotted. Out-of-band spots are flagged with -1. +_BAND_EDGES_MHZ = [ + ( 0.135, 0.138, 2200), # 2200 m (LF experimental, ~137 kHz) + ( 1.8, 2.0, 160), # 160 m + ( 3.5, 4.0, 80), # 80 m + ( 5.25, 5.50, 60), # 60 m (widened to catch 5.288 MHz spots) + ( 7.0, 7.3, 40), # 40 m + (10.1, 10.15, 30), # 30 m + (14.0, 14.35, 20), # 20 m + (18.068, 18.168, 17), # 17 m + (21.0, 21.45, 15), # 15 m + (24.89, 24.99, 12), # 12 m + (28.0, 29.7, 10), # 10 m + (50.0, 54.0, 6), # 6 m + (144.0, 148.0, 2), # 2 m +] + +def _freq_MHz_to_band(freq_MHz): + """Vectorized: MHz -> band wavelength in meters (-1 = out-of-band).""" + out = np.full(len(freq_MHz), -1, dtype=np.int16) + for lo, hi, meters in _BAND_EDGES_MHZ: + out[(freq_MHz >= lo) & (freq_MHz < hi)] = meters + return out + +# Output column -> (raw Madrigal field, numpy dtype, converter). 'occurred', 'freq', +# 'band', and 'source' are synthesized from raw fields rather than mapped 1:1. +_MADRIGAL_NUMERIC_COLS = { + 'dist_Km': ('pthlen', np.float32), + 'tx_lat': ('txlat', np.float32), + 'tx_long': ('txlon', np.float32), + 'rx_lat': ('rxlat', np.float32), + 'rx_long': ('rxlon', np.float32), + 'md_lat': ('latcen', np.float32), + 'md_long': ('loncen', np.float32), + 'snr_dB': ('sn', np.float32), +} +_MADRIGAL_STRING_COLS = { + 'call_tx': 'call_sign_tx', + 'call_rx': 'call_sign_rx', + 'mode': 'smode', +} +_MADRIGAL_ALL_COLS = ( + {'occurred', 'freq', 'band', 'source'} + | set(_MADRIGAL_NUMERIC_COLS) | set(_MADRIGAL_STRING_COLS) +) + +def _detect_ut1_offset(t, probe=1000): + """ + Detect per-file offset (seconds) between ut1_unix and the y/m/d/h/m/s columns. + Returns (offset_sec, uniform). Older Madrigal ham files carry a local-timezone + offset in ut1_unix/ut2_unix; the y/m/d/h/m/s columns are authoritative UTC. + `t` is a pytables Table; only the first `probe` rows are read. + """ + n = min(probe, t.nrows) + ymd = pd.DataFrame({ + 'year': np.char.decode(t.read(stop=n, field='year')).astype(int), + 'month': np.char.decode(t.read(stop=n, field='month')).astype(int), + 'day': np.char.decode(t.read(stop=n, field='day')).astype(int), + 'hour': np.char.decode(t.read(stop=n, field='hour')).astype(int), + 'minute': np.char.decode(t.read(stop=n, field='min')).astype(int), + 'second': np.char.decode(t.read(stop=n, field='sec')).astype(int), + }) + ymd_sec = pd.to_datetime(ymd).values.astype('datetime64[s]').astype(np.int64) + diffs = t.read(stop=n, field='ut1_unix').astype(np.int64) - ymd_sec + uniq = np.unique(diffs) + return int(uniq[0]) if len(uniq) == 1 else int(np.median(diffs)), len(uniq) == 1 + +def _parse_occurred_from_ymd(t): + """Full y/m/d/h/m/s parse (fallback for files with non-uniform ut1_unix offset).""" + return pd.to_datetime(pd.DataFrame({ + 'year': np.char.decode(t.read(field='year')).astype(int), + 'month': np.char.decode(t.read(field='month')).astype(int), + 'day': np.char.decode(t.read(field='day')).astype(int), + 'hour': np.char.decode(t.read(field='hour')).astype(int), + 'minute': np.char.decode(t.read(field='min')).astype(int), + 'second': np.char.decode(t.read(field='sec')).astype(int), + }), utc=True) + +def load_madrigal_hdf5(file_path, columns=None): + """ + Load a Madrigal ham-radio HDF5 file (rsdYYYY-MM-DD.01.hdf5) into a pandas DataFrame. - df = df[tf].copy() + columns: iterable of output column names to load. None (default) loads all. + Available: occurred, freq, band, dist_Km, source, tx_lat, tx_long, rx_lat, + rx_long, md_lat, md_long, call_tx, call_rx, snr_dB, mode. + Selective loading is significantly faster and more memory-efficient on large + files: the pipeline-minimal set {occurred, band, dist_Km, md_lat, md_long, + source} fits a 114M-row day in ~3 GB vs ~78 GB for the full recarray read. - return df + Timestamps: prefers the fast ut1_unix column, but detects the known ut1_unix + local-timezone offset per-file and corrects or falls back as needed. + """ + if columns is None: + cols = set(_MADRIGAL_ALL_COLS) + else: + cols = set(columns) + unknown = cols - _MADRIGAL_ALL_COLS + if unknown: + raise ValueError(f"Unknown columns requested: {sorted(unknown)}. " + f"Available: {sorted(_MADRIGAL_ALL_COLS)}") + + data = {} + with tables.open_file(file_path, 'r') as h: + t = h.get_node('/Data/Table Layout') + + if 'occurred' in cols: + offset_sec, uniform = _detect_ut1_offset(t) + if uniform and offset_sec == 0: + data['occurred'] = pd.to_datetime(t.read(field='ut1_unix'), unit='s', utc=True) + elif uniform: + warnings.warn( + f"{os.path.basename(file_path)}: ut1_unix offset = {offset_sec:+d}s " + "(local-timezone offset). Correcting from ut1_unix.", + RuntimeWarning, stacklevel=2) + data['occurred'] = pd.to_datetime( + t.read(field='ut1_unix') - offset_sec, unit='s', utc=True) + else: + warnings.warn( + f"{os.path.basename(file_path)}: non-uniform ut1_unix offset; " + "falling back to y/m/d/h/m/s parse (slow).", + RuntimeWarning, stacklevel=2) + data['occurred'] = _parse_occurred_from_ymd(t) + + if 'freq' in cols or 'band' in cols: + freq_MHz = (t.read(field='tfreq') / 1.0e6).astype(np.float32) + if 'freq' in cols: + data['freq'] = freq_MHz + if 'band' in cols: + data['band'] = _freq_MHz_to_band(freq_MHz) + + if 'source' in cols: + ssrc_str = np.char.decode(t.read(field='ssrc')) + source_mapped = pd.Series(ssrc_str).map(_MADRIGAL_SSRC_MAP) + n_unknown = source_mapped.isna().sum() + if n_unknown: + unknown_vals = sorted(set(ssrc_str) - set(_MADRIGAL_SSRC_MAP)) + warnings.warn( + f"{os.path.basename(file_path)}: {n_unknown} rows have unknown " + f"ssrc values {unknown_vals}; tagging source=0.", + RuntimeWarning, stacklevel=2) + data['source'] = source_mapped.fillna(0).astype(np.int8).values + + for out_col, (field, dtype) in _MADRIGAL_NUMERIC_COLS.items(): + if out_col in cols: + data[out_col] = t.read(field=field).astype(dtype, copy=False) + + for out_col, field in _MADRIGAL_STRING_COLS.items(): + if out_col in cols: + data[out_col] = np.char.decode(t.read(field=field)) + + return pd.DataFrame(data) def load_spots_csv(date_str,data_sources=[1,2],loc_sources=None, rgc_lim=None,filter_region=None,filter_region_kind='mids'): @@ -488,10 +646,16 @@ def load_spots_csv(date_str,data_sources=[1,2],loc_sources=None, E: Estimated using prefix """ - hdf_path = "data/spot_csvs/{}.hdf".format(date_str) + hdf_path = "data/madrigal/rsd{}.01.hdf5".format(date_str) csv_path = "data/spot_csvs/{}.csv.bz2".format(date_str) if os.path.exists(hdf_path): - df = pd.read_hdf(hdf_path) + # Load only the columns the downstream pipeline actually consumes. tx/rx coords + # are needed only for endpoints-kind regional filtering; skipping them on the + # mids path (the default) roughly halves peak memory on large files. + hdf_cols = {'occurred', 'band', 'dist_Km', 'md_lat', 'md_long', 'source'} + if filter_region_kind == 'endpoints': + hdf_cols |= {'tx_lat', 'tx_long', 'rx_lat', 'rx_long'} + df = load_madrigal_hdf5(hdf_path, columns=hdf_cols) elif os.path.exists(csv_path): df = pd.read_csv(csv_path,parse_dates=['occurred']) else: @@ -503,19 +667,24 @@ def load_spots_csv(date_str,data_sources=[1,2],loc_sources=None, # Select spotting networks if data_sources is not None: - tf = df.source.map(lambda x: x in data_sources) - df = df[tf].copy() + df = df[df['source'].isin(data_sources)].copy() if len(df) == 0: return - # Filter location source + # Filter location source (CSV path only; Madrigal HDF5 does not carry loc_source fields) if loc_sources is not None: - tf = df.tx_loc_source.map(lambda x: x in loc_sources) - df = df[tf].copy() + if 'tx_loc_source' not in df.columns: + warnings.warn( + "loc_sources filter requested but tx_loc_source/rx_loc_source " + "columns absent (Madrigal HDF5 path). Skipping.", + RuntimeWarning, stacklevel=2) + else: + tf = df.tx_loc_source.map(lambda x: x in loc_sources) + df = df[tf].copy() - tf = df.rx_loc_source.map(lambda x: x in loc_sources) - df = df[tf].copy() + tf = df.rx_loc_source.map(lambda x: x in loc_sources) + df = df[tf].copy() if len(df) == 0: return @@ -525,24 +694,29 @@ def load_spots_csv(date_str,data_sources=[1,2],loc_sources=None, tf = np.logical_and(df['dist_Km'] >= rgc_lim[0], df['dist_Km'] < rgc_lim[1]) df = df[tf].copy() - + if len(df) == 0: return - midpoints = geopack.midpoint(df["tx_lat"], df["tx_long"], df["rx_lat"], df["rx_long"]) - df['md_lat'] = midpoints[0] - df['md_long'] = midpoints[1] + # Madrigal HDF5 already carries great-circle midpoints (latcen/loncen); only compute + # for the CSV path. + if 'md_lat' not in df.columns or 'md_long' not in df.columns: + midpoints = geopack.midpoint(df["tx_lat"], df["tx_long"], df["rx_lat"], df["rx_long"]) + df['md_lat'] = midpoints[0] + df['md_long'] = midpoints[1] # Regional Filtering if filter_region is not None: - df_raw = df.copy() - df = regional_filter(filter_region,df,kind=filter_region_kind) + df = regional_filter(filter_region, df, kind=filter_region_kind) if len(df) == 0: return - df["ut_hrs"] = df['occurred'].map(lambda x: x.hour + x.minute/60. + x.second/3600.) - df['slt_mid'] = (df['ut_hrs'] + df['md_long']/15.) % 24. + # Vectorized ut_hrs / slt_mid derivation. The prior .map(lambda) was O(N) in Python + # and dominated load time on multi-million-row days. pandas' .dt accessor is compiled. + occ = df['occurred'] + df["ut_hrs"] = (occ.dt.hour + occ.dt.minute / 60. + occ.dt.second / 3600.).astype(np.float32) + df['slt_mid'] = ((df['ut_hrs'] + df['md_long'] / 15.) % 24.).astype(np.float32) return df diff --git a/harc_plot/goes.py b/harc_plot/goes.py index fe600f9..d7e1b07 100755 --- a/harc_plot/goes.py +++ b/harc_plot/goes.py @@ -33,11 +33,16 @@ import logging import os +import re +import shutil import datetime import fnmatch import glob import ftplib import calendar +import urllib.request +import urllib.error +import warnings import matplotlib matplotlib.use('Agg') @@ -244,6 +249,265 @@ def read_goes(sTime,eTime=None,sat_nr=15,data_dir='data/goes'): data_dict['orbit'].index = [dtGreg_to_datetime(x) for x in data_dict['orbit'].index] return data_dict +# NCEI HTTPS archives. read_goes() above targets ftp://satdat.ngdc.noaa.gov, which NOAA +# retired along with anonymous FTP, so it can no longer download anything -- it fails +# closed, printing 'GOES Data ERROR.' and returning None. These two archives replace it +# and between them span 2009 to present. +# +# Two eras with different file layouts and variable names: +# pre-R (GOES 8-15): one netCDF per month, vars A_AVG / B_AVG, time_tag +# GOES-R (GOES 16-19): one netCDF per day, vars xrsa_flux / xrsb_flux, time +# read_goes_ncei() normalizes both to the A_AVG / B_AVG names the plotting and +# flare-finding routines here already expect. +_NCEI_PRE_R_DIR = ('https://www.ncei.noaa.gov/data/goes-space-environment-monitor/' + 'access/avg/{year:04d}/{month:02d}/goes{sat:02d}/netcdf/') +_NCEI_GOES_R_DIR = ('https://data.ngdc.noaa.gov/platforms/solar-space-observing-satellites/' + 'goes/goes{sat:02d}/l2/data/xrsf-l2-avg1m/{year:04d}/{month:02d}/') + +# GOES-R quality flags: 0 good_data, 1 eclipse, 2 bad_data, 4 interpolated. +# Eclipse and bad_data are dropped; interpolated is kept. +_GOES_R_BAD_FLAGS = (1, 2) + +def goes_xrs_candidates(year): + """Prioritized [(sat_nr, era)] candidates carrying XRS for a given year. + + Ordered most- to least-preferred; read_goes_ncei() tries them in turn and takes + the first that actually has data, so these windows only need to be roughly right. + 'era' selects the archive layout: 'pre_r' or 'goes_r'. + + Ordering follows the *primary operational* satellite for each epoch, because that + is the one whose fluxes reproduce the standard flare classifications. This matters: + for the 2017-09-06 flare, GOES-15 (operational primary) peaks at 9.33e-4 W/m^2, + i.e. X9.3, matching the accepted classification, while GOES-16 -- then still in + post-launch test and reporting on a different flux scale -- peaks at 5.17e-4, i.e. + X5.2. GOES-13 reads X10.4. So GOES-15 is preferred through its 2020-03 end of life + even after GOES-16 data becomes available. + + Corollary worth remembering when comparing events across the archive: X-class + numbers are NOT directly comparable between the pre-R and GOES-R eras. Flare + magnitudes from 2015 and 2024 sit on different scales. + """ + if year >= 2025: return [(19,'goes_r'),(18,'goes_r'),(16,'goes_r')] + if year >= 2023: return [(18,'goes_r'),(16,'goes_r')] + # GOES-15 XRS ended 2020-03-04, so 2020 leads with the GOES-R series but keeps 15 + # as a fallback for January through early March. + if year >= 2020: return [(16,'goes_r'),(17,'goes_r'),(18,'goes_r'),(15,'pre_r')] + if year >= 2010: return [(15,'pre_r'),(14,'pre_r'),(13,'pre_r')] + return [(14,'pre_r'),(11,'pre_r'),(12,'pre_r'),(10,'pre_r')] + +def _ncei_listing(url,timeout=60): + """Filenames in an NCEI Apache directory index. [] if unreachable.""" + try: + with urllib.request.urlopen(url,timeout=timeout) as resp: + html = resp.read().decode('utf-8','replace') + except (urllib.error.URLError,OSError) as e: + logging.info('NCEI listing failed for {0}: {1}'.format(url,e)) + return [] + # Apache indexes repeat each name in href and link text; dedupe but keep order. + names,seen = [],set() + for name in re.findall(r'href="([^"?/][^"]*\.nc)"',html): + if name not in seen: + seen.add(name) + names.append(name) + return names + +def _ncei_download(url,dest,timeout=300): + """Download to dest, returning True on success. + + Writes to a .part file and renames, so an interrupted download never leaves a + truncated file that later runs would treat as a valid cache hit. + """ + part = dest + '.part' + try: + with urllib.request.urlopen(url,timeout=timeout) as resp, open(part,'wb') as fh: + shutil.copyfileobj(resp,fh) + os.replace(part,dest) + return True + except (urllib.error.URLError,OSError) as e: + logging.info('NCEI download failed for {0}: {1}'.format(url,e)) + if os.path.exists(part): + os.remove(part) + return False + +def _months_in_range(sTime,eTime): + """First-of-month dates covering [sTime,eTime).""" + months = [datetime.date(sTime.year,sTime.month,1)] + end = datetime.date(eTime.year,eTime.month,1) + while months[-1] < end: + months.append(add_months(months[-1])) + return months + +def _read_pre_r_nc(file_path): + """(DataFrame[A_AVG,B_AVG], metadata dict) from a pre-R monthly XRS file.""" + with netCDF4.Dataset(file_path) as nc: + tt = nc.variables['time_tag'] + index = pd.to_datetime(netCDF4.num2date( + tt[:],tt.units,only_use_cftime_datetimes=False)) + data = {} + for out_col,var in (('A_AVG','A_AVG'),('B_AVG','B_AVG')): + arr = np.ma.filled(nc.variables[var][:].astype(float),np.nan) + arr[arr == -99999] = np.nan + data[out_col] = arr + md = {k:getattr(nc,k,None) for k in + ('institution','satellite_id','instrument','creation_date')} + return pd.DataFrame(data,index=index),md + +def _read_goes_r_nc(file_path): + """(DataFrame[A_AVG,B_AVG], metadata dict) from a GOES-R daily XRS file. + + xrsa_flux/xrsb_flux are renamed to A_AVG/B_AVG so downstream code that predates + the GOES-R series keeps working unchanged. The wavelength bands are equivalent: + XRS-A is 0.05-0.4 nm, XRS-B is 0.1-0.8 nm. + """ + with netCDF4.Dataset(file_path) as nc: + tt = nc.variables['time'] + index = pd.to_datetime(netCDF4.num2date( + tt[:],tt.units,only_use_cftime_datetimes=False)) + data = {} + for out_col,var,flag_var in (('A_AVG','xrsa_flux','xrsa_flag'), + ('B_AVG','xrsb_flux','xrsb_flag')): + arr = np.ma.filled(nc.variables[var][:].astype(float),np.nan) + if flag_var in nc.variables: + flags = np.ma.filled(nc.variables[flag_var][:],0) + arr[np.isin(flags,_GOES_R_BAD_FLAGS)] = np.nan + data[out_col] = arr + md = {'institution': getattr(nc,'institution',None), + 'satellite_id': getattr(nc,'platform',None), + 'instrument': getattr(nc,'instrument',None), + 'creation_date':getattr(nc,'date_created',None)} + return pd.DataFrame(data,index=index),md + +def read_goes_ncei(sTime,eTime=None,sat_nr=None,era=None,data_dir='data/goes'): + """Load 1-minute GOES XRS data from the NOAA NCEI HTTPS archives. + + Drop-in replacement for read_goes() with the same return contract, but backed by + endpoints that still exist and spanning 2009 to present rather than only the + GOES-13/15 era. Downloaded files are cached in data_dir and reused. + + X-ray flux resolves individual solar flares (1-minute cadence, and the 0.1-0.8 nm + B channel is what defines the A/B/C/M/X classes). This is the observable to use for + flare timing -- F10.7 cannot serve that purpose at any cadence, because DRAO + Penticton determines it once per day. + + Parameters + ---------- + sTime : datetime.datetime + Start of the interval. + eTime : Optional[datetime.datetime] + End of the interval, exclusive. Defaults to sTime + 1 day. + sat_nr : Optional[int] + Force a specific GOES satellite. Default None auto-selects per year via + goes_xrs_candidates(), trying each until one yields data. + era : Optional[str] + 'pre_r' or 'goes_r'. Only consulted when sat_nr is given explicitly. + data_dir : Optional[str] + Local cache directory. + + Returns + ------- + Dict with the same shape read_goes() returns: 'xray' (DataFrame indexed by + datetime with A_AVG, B_AVG, ut_hrs) and 'metadata'. None if no data could be + obtained for the interval. + + Note that 'orbit' data, and therefore the 'slt_mid' column, are not provided -- + the GOES-R XRS product does not carry the spacecraft ephemeris that read_goes() + pulled from the pre-R files. + """ + if eTime is None: + eTime = sTime + datetime.timedelta(days=1) + + os.makedirs(data_dir,exist_ok=True) + + if sat_nr is not None: + # An explicit satellite with no era stated: infer from the series number. + candidates = [(sat_nr, era if era is not None + else ('goes_r' if sat_nr >= 16 else 'pre_r'))] + else: + # Candidates are year-keyed; a window spanning a year boundary needs both. + candidates = [] + for year in sorted({d.year for d in _months_in_range(sTime,eTime)}): + for cand in goes_xrs_candidates(year): + if cand not in candidates: + candidates.append(cand) + + frames,metadata = [],{} + used_sats = set() + for cand_sat,cand_era in candidates: + got_any = False + for month in _months_in_range(sTime,eTime): + if cand_era == 'pre_r': + base = _NCEI_PRE_R_DIR.format(year=month.year,month=month.month,sat=cand_sat) + pattern = 'g{0:02d}_xrs_1m_*.nc'.format(cand_sat) + reader = _read_pre_r_nc + wanted = None # monthly file; take whatever is there + else: + base = _NCEI_GOES_R_DIR.format(year=month.year,month=month.month,sat=cand_sat) + pattern = 'dn_xrsf-l2-avg1m_g{0:02d}_d*.nc'.format(cand_sat) + reader = _read_goes_r_nc + # Daily files: fetch only the days actually in range. + wanted = set() + day = datetime.date(sTime.year,sTime.month,sTime.day) + while day < eTime.date() + datetime.timedelta(days=1): + wanted.add(day.strftime('%Y%m%d')) + day += datetime.timedelta(days=1) + + cached = [os.path.basename(p) for p in + glob.glob(os.path.join(data_dir,pattern))] + names = [n for n in _ncei_listing(base) if fnmatch.fnmatch(n,pattern)] + names = sorted(set(names) | set(cached)) + if wanted is not None: + names = [n for n in names + if any('_d{0}_'.format(d) in n for d in wanted)] + + for name in names: + path = os.path.join(data_dir,name) + if not os.path.exists(path): + if not _ncei_download(base + name,path): + continue + try: + df,md = reader(path) + except Exception as e: + # A corrupt cached file should not kill the whole request; drop it + # so the next run re-downloads. + warnings.warn('{0}: unreadable GOES file ({1}: {2}); removing from ' + 'cache'.format(name,type(e).__name__,e),RuntimeWarning) + try: + os.remove(path) + except OSError: + pass + continue + frames.append(df) + metadata[name] = md + got_any = True + + if got_any: + used_sats.add(cand_sat) + break # first candidate with data wins; don't stack satellites + + if not frames: + warnings.warn('No GOES XRS data available from NCEI for {0} to {1} ' + '(tried satellites {2}).'.format( + sTime,eTime,[c[0] for c in candidates]),RuntimeWarning) + return + + df_xray = pd.concat(frames).sort_index() + df_xray = df_xray[~df_xray.index.duplicated(keep='first')] + df_xray = df_xray[(df_xray.index >= sTime) & (df_xray.index < eTime)] + if len(df_xray) == 0: + warnings.warn('GOES XRS files found but no samples fall in {0} to {1}.'.format( + sTime,eTime),RuntimeWarning) + return + + df_xray['ut_hrs'] = df_xray.index.map(ut_hours) + + data_dict = {'xray':df_xray,'metadata':metadata} + data_dict['metadata']['variables'] = { + 'A_AVG': {'long_label':'x-ray (0.05-0.4 nm) irradiance','units':'W/m^2'}, + 'B_AVG': {'long_label':'x-ray (0.1-0.8 nm) irradiance', 'units':'W/m^2'}, + } + data_dict['sat_nrs'] = sorted(used_sats) + return data_dict + def goes_plot_hr(goes_data,ax,var_tags = ['B_AVG'],xkey='ut_hr',xlim=(0,24),ymin=1e-9,ymax=1e-2, legendSize=10,legendLoc=None,labels=None,**kwargs): """Plot GOES X-Ray Data. @@ -588,13 +852,21 @@ def find_flares(goes_data,window_minutes=60,min_class='X1',sTime=None,eTime=None sWin = win - time_delta_half eWin = win + time_delta_half + # An all-NaN window (eclipse, data gap, flagged-bad stretch) makes idxmax + # return pd.NaT, and `NaT is np.nan` is False -- so the old identity check let + # NaT through into `keys` and the b_avg[keys] lookup below raised + # "KeyError: '[NaT] not in index'". Test for emptiness up front, which also + # avoids the deprecated all-NA idxmax call. + win_vals = b_avg[sWin:eWin] + if len(win_vals) == 0 or win_vals.isna().all(): + continue try: - idx_max = b_avg[sWin:eWin].idxmax() - if idx_max is np.nan: - continue - keys.append(idx_max) - except: - pass + idx_max = win_vals.idxmax() + except (ValueError,TypeError): + continue + if pd.isna(idx_max): + continue + keys.append(idx_max) df_win = pd.DataFrame({'B_AVG':b_avg[keys]}) diff --git a/harc_plot/omni.py b/harc_plot/omni.py index 52f51f9..d2530c1 100644 --- a/harc_plot/omni.py +++ b/harc_plot/omni.py @@ -196,7 +196,25 @@ def plot_dst_kp(self,sTime,eTime,ax,xkey='index',xlabels=True, lines.append(tmp) ax.set_ylabel(ylabel) ax.axhline(0,color='k',ls='--') - ax.set_ylim(-200,100) + + # Expand-only y-limits. (-200,100) is a floor, not a cap: intense storms go + # below -200 nT (-223 on 2015-03-17, -412 on 2024-05-10) and were previously + # clipped off the axis entirely. Limits are computed over the plotted window + # only -- the SYM-H branch above hands matplotlib the full multi-decade series + # and relies on set_xlim to window it, so a naive min/max would pick up the + # all-time extreme. Plots whose data already fits are unchanged. + ylim = [-200.,100.] + win = pd.Series(yy,index=xx) + win = win[(win.index >= xlim[0]) & (win.index <= xlim[1])] + if len(win) > 0: + pad = 20. + y_min = np.nanmin(win.values) + y_max = np.nanmax(win.values) + if np.isfinite(y_min): + ylim[0] = min(ylim[0], float(y_min) - pad) + if np.isfinite(y_max): + ylim[1] = max(ylim[1], float(y_max) + pad) + ax.set_ylim(*ylim) ax.set_xlim(xlim) # Kp ################################### diff --git a/harc_plot/visualize_histograms.py b/harc_plot/visualize_histograms.py index a055dd1..401b3e8 100755 --- a/harc_plot/visualize_histograms.py +++ b/harc_plot/visualize_histograms.py @@ -4,6 +4,7 @@ from collections import OrderedDict import string import ast +import warnings import matplotlib as mpl mpl.use('Agg') @@ -214,10 +215,18 @@ def plot(self,group,ax,ls='--',lw=4,color='white'): sza_ax.set_ylim(110,0) class Goeser(object): - def __init__(self,sTime,eTime,sats=[13,15]): + # 1-minute X-ray flux is meaningless on a multi-year axis, and fetching it for one + # would pull thousands of daily files. Above this span the panel is skipped. + MAX_SPAN_DAYS = 120 + + def __init__(self,sTime,eTime,sats=None,max_span_days=None): + """sats=None (default) auto-selects the primary operational satellite for the + epoch via goes.goes_xrs_candidates(). Pass an explicit list to override.""" self.sTime = sTime self.eTime = eTime self.sats = sats + self.max_span_days = (self.MAX_SPAN_DAYS if max_span_days is None + else max_span_days) self.load_goes() @@ -226,42 +235,67 @@ def load_goes(self): eTime = self.eTime sats = self.sats - goes_dcts = OrderedDict() - for sat in sats: - goes_dcts[sat] = {} + self.flares_combined = pd.DataFrame() + self.goes_dcts = OrderedDict() + + span_days = (eTime - sTime).total_seconds() / 86400. + if span_days > self.max_span_days: + warnings.warn('GOES X-ray panel skipped: {0:.0f}-day span exceeds ' + '{1}-day limit.'.format(span_days,self.max_span_days), + RuntimeWarning) + return - flares_combined = pd.DataFrame() - for sat_nr,gd in goes_dcts.items(): - gd['data'] = goes.read_goes(sTime,eTime,sat_nr=sat_nr) + # sats=None lets read_goes_ncei pick per epoch; otherwise honor the request. + # Either way this is one entry, not a stack of satellites: mixing pre-R and + # GOES-R traces on one axis would overlay incompatible flux scales. + requests = [None] if sats is None else list(sats) + + for sat_req in requests: + gd = {} + gd['data'] = goes.read_goes_ncei(sTime,eTime,sat_nr=sat_req) # Skip if no GOES data present if gd['data'] is None: continue + sat_nr = gd['data'].get('sat_nrs',[sat_req]) + sat_nr = sat_nr[0] if sat_nr else sat_req + flares = goes.find_flares(gd['data'],min_class='M1',window_minutes=60) flares['sat'] = sat_nr gd['flares'] = flares - flares_combined = pd.concat([flares_combined,flares]).sort_index() + self.flares_combined = pd.concat([self.flares_combined,flares]).sort_index() gd['var_tags'] = ['B_AVG'] gd['labels'] = ['GOES {!s}'.format(sat_nr)] + self.goes_dcts[sat_nr] = gd - self.flares_combined = flares_combined[~flares_combined.index.duplicated()].sort_index() - self.goes_dcts = goes_dcts + if len(self.flares_combined) > 0: + self.flares_combined = self.flares_combined[ + ~self.flares_combined.index.duplicated()].sort_index() def plot(self,ax,lw=2): sTime = self.sTime eTime = self.eTime + plotted = 0 for sat_nr,gd in self.goes_dcts.items(): if gd['data'] is None: continue goes.goes_plot(gd['data'],sTime,eTime,ax=ax, var_tags=gd['var_tags'],labels=gd['labels'], legendLoc='upper right',lw=lw) + plotted += 1 title = 'NOAA GOES X-Ray (0.1 - 0.8 nm) Irradiance' size = 20 ax.text(0.01,0.05,title,transform=ax.transAxes,ha='left',fontdict={'size':size,'weight':'bold'}) + if plotted == 0: + # Say so explicitly. A blank log-scale panel otherwise reads as "no flares + # occurred" when it actually means no data was retrieved for this interval. + ax.text(0.5,0.5,'No GOES X-ray data available for this interval', + transform=ax.transAxes,ha='center',va='center', + fontdict={'size':size,'weight':'bold'},color='0.4') + class ncLoader(object): def __init__(self,sTime,eTime=None,srcs=None,band_keys=None,**kwargs): @@ -444,8 +478,24 @@ def _format_timeticklabels(self,ax): def plot(self,baseout_dir='output',xlim=None,ylim=None,xunits='datetime', plot_sza=True,subdir=None,geospace_env=None,plot_region=None, - plot_kpsymh=True,plot_goes=True,plot_f107=False,axvlines=None,axvlines_kw={},axvspans=None,time_format={}, + plot_kpsymh=True,plot_goes=True,plot_f107=False,dst_param='SYM-H', + context_panels=('kpsymh','goes','f107'), + axvlines=None,axvlines_kw={},axvspans=None,time_format={}, xkeys=None,log_z=None,**kwargs): + # Validate before the no-data early return, so a mistyped key is caught on every + # call rather than only on the ones that happen to have data. A typo would + # otherwise silently drop a panel, and a duplicate would draw one twice into + # the same row. + _valid_panels = ('kpsymh','goes','f107') + context_panels = tuple(context_panels) + bad = [p for p in context_panels if p not in _valid_panels] + if bad: + raise ValueError('Unknown context_panels entry/entries {0}. ' + 'Valid keys: {1}'.format(bad,list(_valid_panels))) + if len(set(context_panels)) != len(context_panels): + raise ValueError('context_panels contains duplicates: ' + '{0}'.format(list(context_panels))) + if self.datasets is None: return @@ -529,30 +579,37 @@ def plot(self,baseout_dir='output',xlim=None,ylim=None,xunits='datetime', axs_to_adjust = [] - pinx = -1 - if plot_kpsymh: - pinx += 1 + ######################################## + # Context panels, stacked above the per-band panels. Each renders into + # the row it is handed, so the caller controls the order via + # context_panels; the panel letters follow whatever order is given. + def _panel_kpsymh(pinx): ax = plt.subplot2grid((ny,nx),(pinx,col_1),colspan=col_1_span) ax.set_xlim(xlim) ax.set_ylim(ylim) + # dst_param selects the ring-current index: 'SYM-H' is 1-min Kyoto + # data but is only loaded for SymH(years=...) (2016-2017 by default, + # and no Kyoto file exists past 2019); 'Dst_nT' is hourly OMNI and + # spans the whole archive. Outside the SYM-H years the panel plots + # empty, so callers working across arbitrary epochs want 'Dst_nT'. omni_axs = geospace_env.omni.plot_dst_kp(self.sTime,self.eTime,ax,xlabels=True, - kp_markersize=10,dst_lw=2,dst_param='SYM-H') + kp_markersize=10,dst_lw=2,dst_param=dst_param) ax.tick_params(**tick_params) plot_letter(pinx,ax) label_time = axvlines_kw.get('label_time',True) plot_axv(axvlines,omni_axs[0],color='k',label_time=label_time) plot_axvspans(axvspans,omni_axs[0]) - for ax in omni_axs: - self._format_timeticklabels(ax) - ax.set_xlabel('') - axs_to_adjust += omni_axs - - ######################################## - if plot_goes: + for omni_ax in omni_axs: + self._format_timeticklabels(omni_ax) + omni_ax.set_xlabel('') + # extend(), not +=: augmented assignment would rebind the closed-over + # name as a local and raise UnboundLocalError. + axs_to_adjust.extend(omni_axs) + + def _panel_goes(pinx): goeser = Goeser(self.sTime,self.eTime) - pinx +=1 ax = plt.subplot2grid((ny,nx),(pinx,col_1),colspan=col_1_span) goeser.plot(ax) ax.set_xlim(xlim) @@ -564,9 +621,7 @@ def plot(self,baseout_dir='output',xlim=None,ylim=None,xunits='datetime', self._format_timeticklabels(ax) ax.set_xlabel('') - ######################################## - if plot_f107: - pinx +=1 + def _panel_f107(pinx): ax = plt.subplot2grid((ny,nx),(pinx,col_1),colspan=col_1_span) geospace_env.omni.plot_f107(self.sTime,self.eTime,ax,xlabels=True) ax.set_xlim(xlim) @@ -577,7 +632,21 @@ def plot(self,baseout_dir='output',xlim=None,ylim=None,xunits='datetime', axs_to_adjust.append(ax) self._format_timeticklabels(ax) ax.set_xlabel('') - + + panel_funcs = {'kpsymh':_panel_kpsymh, + 'goes': _panel_goes, + 'f107': _panel_f107} + panel_on = {'kpsymh':plot_kpsymh, + 'goes': plot_goes, + 'f107': plot_f107} + + pinx = -1 + for panel_key in context_panels: + if not panel_on.get(panel_key,False): + continue + pinx += 1 + panel_funcs[panel_key](pinx) + map_sum = 0 for inx,freq in enumerate(freqs): plt_row = inx+pinx+1 diff --git a/requirements.txt b/requirements.txt index 49f07ff..789d455 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,6 +3,7 @@ tqdm numpy pandas h5py +tables cartopy ipdb library diff --git a/scripts/basic_histogram/histograms.py b/scripts/basic_histogram/histograms.py index 2bf8d9e..f66449b 100755 --- a/scripts/basic_histogram/histograms.py +++ b/scripts/basic_histogram/histograms.py @@ -12,8 +12,8 @@ plot_dir = os.path.join('output/galleries/histograms',run_name) params = ['spot_density'] xkeys = ['ut_hrs','slt_mid'] -sTime = datetime.datetime(2017,7,1) -eTime = datetime.datetime(2017,7,2) +sTime = datetime.datetime(2015,6,10) +eTime = datetime.datetime(2015,6,11) region = run_name rgc_lim = (0, 10000) @@ -30,6 +30,7 @@ rd['filter_region_kind'] = 'mids' rd['xb_size_min'] = 10. rd['yb_size_km'] = 250. +rd['data_sources'] = [1, 2, 3] # WSPRNet, RBN, PSKReporter rd['reprocess'] = True rd['output_dir'] = data_dir rd['band_obj'] = harc_plot.gl.BandData() diff --git a/scripts/kp_stations/.ipynb_checkpoints/Kp Stations AACGM-checkpoint.ipynb b/scripts/kp_stations/.ipynb_checkpoints/Kp Stations AACGM-checkpoint.ipynb deleted file mode 100644 index 2529b04..0000000 --- a/scripts/kp_stations/.ipynb_checkpoints/Kp Stations AACGM-checkpoint.ipynb +++ /dev/null @@ -1,310 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import aacgmv2\n", - "import datetime\n", - "import numpy as np\n", - "import pandas as pd" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "dtime = datetime.datetime(2017,9,1)\n", - "\n", - "fname = 'kp_stations.csv'\n", - "df = pd.read_csv(fname)\n", - "\n", - "df_lst = []\n", - "for rinx,row in df.iterrows():\n", - " aa = aacgmv2.get_aacgm_coord(row['glat'],row['glon'],0,dtime)\n", - " row['mlat'] = aa[0]\n", - " row['mlon'] = aa[1]\n", - " row['mlt'] = aa[2]\n", - " df_lst.append(row)\n", - " \n", - "df = pd.DataFrame(df_lst)" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - "
codecitycountryglatglonmlatmlonmlt
0LERLerwickScotland60.1358.857.42676279.8517530.532059
1MEAMeanookCanada54.6246.761.334882-50.66810515.830735
2SITSitkaAlaska57.1224.759.503863-76.45736514.111451
3ESKEskdalemuirScotland55.3356.851.93787176.4487230.305190
4UPSUppsalaSweden59.917.456.43896695.0845481.547578
5OTTOttawaCanada45.4284.554.1767443.20862619.422517
6BFEBrorfeldeDenmark55.611.751.73546188.8302241.130623
7HADHartlandEngland51.0355.546.73157674.1168490.149732
8WNGWingstGermany53.89.149.69100386.1052350.948957
9NGKNiemegkGermany52.112.747.67627488.7153241.122963
10FRDFredericksburgUSA38.2282.647.4321110.05330319.212162
11CNBCanberraAustralia-34.7149.0-44.348957-132.75533810.358253
12EYREyrewellNew Zealand-42.6172.4-48.971607-102.90351612.348374
\n", - "
" - ], - "text/plain": [ - " code city country glat glon mlat mlon \\\n", - "0 LER Lerwick Scotland 60.1 358.8 57.426762 79.851753 \n", - "1 MEA Meanook Canada 54.6 246.7 61.334882 -50.668105 \n", - "2 SIT Sitka Alaska 57.1 224.7 59.503863 -76.457365 \n", - "3 ESK Eskdalemuir Scotland 55.3 356.8 51.937871 76.448723 \n", - "4 UPS Uppsala Sweden 59.9 17.4 56.438966 95.084548 \n", - "5 OTT Ottawa Canada 45.4 284.5 54.176744 3.208626 \n", - "6 BFE Brorfelde Denmark 55.6 11.7 51.735461 88.830224 \n", - "7 HAD Hartland England 51.0 355.5 46.731576 74.116849 \n", - "8 WNG Wingst Germany 53.8 9.1 49.691003 86.105235 \n", - "9 NGK Niemegk Germany 52.1 12.7 47.676274 88.715324 \n", - "10 FRD Fredericksburg USA 38.2 282.6 47.432111 0.053303 \n", - "11 CNB Canberra Australia -34.7 149.0 -44.348957 -132.755338 \n", - "12 EYR Eyrewell New Zealand -42.6 172.4 -48.971607 -102.903516 \n", - "\n", - " mlt \n", - "0 0.532059 \n", - "1 15.830735 \n", - "2 14.111451 \n", - "3 0.305190 \n", - "4 1.547578 \n", - "5 19.422517 \n", - "6 1.130623 \n", - "7 0.149732 \n", - "8 0.948957 \n", - "9 1.122963 \n", - "10 19.212162 \n", - "11 10.358253 \n", - "12 12.348374 " - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "df" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "44.34895701713687\n", - "61.334882313036495\n" - ] - } - ], - "source": [ - "print(np.min(np.abs(df['mlat'])))\n", - "print(np.max(np.abs(df['mlat'])))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/scripts/solar_zenith/.ipynb_checkpoints/Astropy Solar Zenith-checkpoint.ipynb b/scripts/solar_zenith/.ipynb_checkpoints/Astropy Solar Zenith-checkpoint.ipynb deleted file mode 100644 index 9cd20b4..0000000 --- a/scripts/solar_zenith/.ipynb_checkpoints/Astropy Solar Zenith-checkpoint.ipynb +++ /dev/null @@ -1,113 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline\n", - "import numpy as np\n", - "import matplotlib as mpl\n", - "from matplotlib import pyplot as plt\n", - "import datetime\n", - "import pandas as pd\n", - "\n", - "import astropy.units as u\n", - "from astropy.time import Time\n", - "from astropy.coordinates import get_sun, EarthLocation, AltAz\n", - "\n", - "from library import gen_lib as gl" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "class SolarZenith(object):\n", - " def __init__(self,sTime,eTime,lon=-100,lat=40,tstep='1M'):\n", - " dt = pd.Timedelta(tstep)\n", - " time_vec = [sTime]\n", - " while time_vec[-1] < eTime:\n", - " time_vec.append(time_vec[-1]+dt)\n", - "\n", - " self.time_vec = time_vec\n", - " self.times_aa = Time(time_vec)\n", - " self.df = pd.DataFrame(index=time_vec)\n", - " self.lat = lat\n", - " self.lon = lon\n", - " \n", - " def calc_astropy(self,alt_km=0):\n", - " loc = EarthLocation(lon=self.lon*u.deg,lat=self.lat*u.deg,height=alt_km*1e3*u.m)\n", - " altaz = AltAz(location=loc,obstime=self.times_aa)\n", - "\n", - " sun = get_sun(self.times_aa)\n", - " zen = sun.transform_to(altaz).zen\n", - " \n", - " key = 'AstroPy {!s} km'.format(alt_km)\n", - " self.df[key] = np.array(zen)\n", - " \n", - " def calc_geopack(self):\n", - " azs, els = gl.sunAzEl(self.time_vec,self.lat,self.lon)\n", - " key = 'Geopack 0 km'\n", - " self.df[key] = els\n", - " \n", - " def plot(self):\n", - " fig = plt.figure(figsize=(15,10))\n", - " ax = fig.add_subplot(1,1,1)\n", - " self.df.plot(ax=ax)\n", - " title = '({:0.1f}$\\degree$ N, {:0.1f}$\\degree$ E)'.format(self.lat,self.lon)\n", - " ax.set_title(title)\n", - " ax.grid()\n", - " plt.show()\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "sTime = datetime.datetime(2017,9,6)\n", - "eTime = datetime.datetime(2017,9,7)\n", - "sz = SolarZenith(sTime,eTime)\n", - "sz.calc_geopack()\n", - "sz.calc_astropy(0)\n", - "sz.calc_astropy(300)\n", - "sz.plot()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/scripts/solarcycle/batch_histograms.py b/scripts/solarcycle/batch_histograms.py new file mode 100755 index 0000000..4ae17d4 --- /dev/null +++ b/scripts/solarcycle/batch_histograms.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python3 +""" +Parallel batch-histogram processor for the full Madrigal ham-radio archive. + +Runs `harc_plot.calculate_histograms.main()` across a date range using a process +pool, producing one `YYYYMMDD.data.nc.bz2` per day under a config-named output +directory. Each worker processes one day at a time; days with existing output +are skipped unless --reprocess is given. + +Networks are kept separate: choose one via --network (WSPRNet, RBN, PSKReporter, +or All for the combined view). For the three-networks-separate workflow, invoke +this script three times. + +Typical usage for the full 5.1 TB pass: + + ./batch_histograms.py --start 2009-01-01 --stop 2027-01-01 \\ + --network WSPRNet --workers 8 +""" +from __future__ import annotations + +import argparse +import datetime +import multiprocessing as mp +import os +import sys +import time +from functools import partial + +import dateutil.parser +import tqdm + +import harc_plot + +# Map canonical network name -> harc_plot source integer (see gen_lib.sources). +NETWORK_TO_SOURCES = { + 'WSPRNet': [1], + 'RBN': [2], + 'PSKReporter': [3], + 'All': [1, 2, 3], +} + + +def _process_day(date, *, data_sources, region, rgc_lim, xkeys, params, + xb_size_min, yb_size_km, data_dir): + """Worker: process a single UTC day. Always reprocess=False here; the driver + handles clearing the output dir once up front so workers do not race on it. + + Exceptions are caught and returned as (date, error_message) so one bad Madrigal + file does not tear down the whole pool. Observed failure modes include truncated + HDF5 files (HDF5ExtError) and files missing /Data/Table Layout (NoSuchNodeError) + — all upstream ingestion issues at MIT Haystack, not bugs in the pipeline.""" + rd = { + 'sDate': date, + 'eDate': date + datetime.timedelta(days=1), + 'params': params, + 'xkeys': xkeys, + 'rgc_lim': rgc_lim, + 'filter_region': region, + 'filter_region_kind': 'mids', + 'xb_size_min': xb_size_min, + 'yb_size_km': yb_size_km, + 'loc_sources': None, + 'data_sources': data_sources, + 'reprocess': False, + 'output_dir': data_dir, + 'band_obj': harc_plot.gl.BandData(), + } + try: + harc_plot.calculate_histograms.main(rd) + return (date, None) + except Exception as e: + return (date, f'{type(e).__name__}: {e}') + + +def _daterange(start, stop): + """Inclusive [start, stop) daily range.""" + d = start + while d < stop: + yield d + d = d + datetime.timedelta(days=1) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('--start', required=True, help='Inclusive start date (YYYY-MM-DD)') + parser.add_argument('--stop', required=True, help='Exclusive stop date (YYYY-MM-DD)') + parser.add_argument('--network', choices=NETWORK_TO_SOURCES, default='All', + help='Network to include (default: All)') + parser.add_argument('--workers', type=int, default=max(1, mp.cpu_count() // 2), + help='Parallel worker processes (default: half of CPU count)') + parser.add_argument('--region', default='World', + help='Regional filter from harc_plot.gl.regions (default: World)') + parser.add_argument('--rgc-min', type=float, default=0.0, help='Min path length, km') + parser.add_argument('--rgc-max', type=float, default=10000.0, help='Max path length, km') + parser.add_argument('--xb-min', type=float, default=60.0, + help='Time bin size, minutes (default: 60 for solarcycle overview)') + parser.add_argument('--yb-km', type=float, default=250.0, help='Range bin size, km') + parser.add_argument('--xkey', default='slt_mid', choices=['ut_hrs', 'slt_mid'], + help='Time-axis coordinate (default: slt_mid)') + parser.add_argument('--output-root', default='data/solarcycle/histograms', + help='Parent dir for per-config output subfolders') + parser.add_argument('--reprocess', action='store_true', + help='Overwrite existing day outputs instead of skipping') + args = parser.parse_args(argv) + + sDate = dateutil.parser.isoparse(args.start) + eDate = dateutil.parser.isoparse(args.stop) + if sDate >= eDate: + parser.error(f'--start ({args.start}) must be before --stop ({args.stop})') + + dates = list(_daterange(sDate, eDate)) + + # Config-keyed output dir: this makes runs for different configs coexist cleanly. + run_name = (f'{args.region}-{args.network}' + f'-{int(args.xb_min)}min-{int(args.yb_km)}km' + f'-{int(args.rgc_min)}-{int(args.rgc_max)}km') + data_dir = os.path.join(args.output_root, run_name) + # Clear the output dir exactly once up front (if --reprocess), then workers + # always run with reprocess=False and simply skip existing per-day outputs. + harc_plot.gl.prep_output({0: data_dir}, clear=args.reprocess) + + worker = partial( + _process_day, + data_sources = NETWORK_TO_SOURCES[args.network], + region = args.region, + rgc_lim = (args.rgc_min, args.rgc_max), + xkeys = [args.xkey], + params = ['spot_density'], + xb_size_min = args.xb_min, + yb_size_km = args.yb_km, + data_dir = data_dir, + ) + + print(f'Processing {len(dates)} days × network={args.network} × region={args.region}', + file=sys.stderr) + print(f' workers: {args.workers}', file=sys.stderr) + print(f' output dir: {data_dir}', file=sys.stderr) + print(f' reprocess: {args.reprocess}', file=sys.stderr) + + t0 = time.time() + failed = [] # list of (date, error_message) + # spawn context avoids fork-safety issues with pytables / HDF5 file handles. + ctx = mp.get_context('spawn') + with ctx.Pool(args.workers) as pool: + it = pool.imap_unordered(worker, dates, chunksize=1) + for result_date, err in tqdm.tqdm(it, total=len(dates), desc='days', unit='day'): + if err is not None: + failed.append((result_date, err)) + dt = time.time() - t0 + print(f'Done in {dt:.1f}s ({dt/max(len(dates),1):.1f}s/day avg).', file=sys.stderr) + if failed: + print(f'{len(failed)} days failed:', file=sys.stderr) + for d, err in sorted(failed): + print(f' {d.strftime("%Y-%m-%d")} {err}', file=sys.stderr) + # Also write a machine-readable list next to the output dir for later triage. + fail_path = os.path.join(data_dir, '_failed_dates.txt') + with open(fail_path, 'w') as fh: + for d, err in sorted(failed): + fh.write(f'{d.strftime("%Y-%m-%d")}\t{err}\n') + print(f'Wrote failure list to {fail_path}', file=sys.stderr) + + +if __name__ == '__main__': + main() diff --git a/scripts/solarcycle/plot_solarcycle.py b/scripts/solarcycle/plot_solarcycle.py new file mode 100755 index 0000000..d4c307e --- /dev/null +++ b/scripts/solarcycle/plot_solarcycle.py @@ -0,0 +1,148 @@ +#!/usr/bin/env python3 +""" +Aggregate per-day histogram NetCDFs produced by batch_histograms.py into a +single multi-year solar-cycle panel plot. + +For each --network, reads every YYYYMMDD.data.nc.bz2 under the matching +config-keyed input dir, concatenates into an xarray Dataset, and renders +time-series + map panels via harc_plot.visualize_histograms.ncLoader. + +The loaded-and-concatenated Dataset is pickled to a cache file so iterating +on plot styling does not re-read all ~4000 daily NetCDFs every time. + +Typical usage: + + ./plot_solarcycle.py --network WSPRNet + ./plot_solarcycle.py --network RBN --start 2009-02-25 + ./plot_solarcycle.py --network PSKReporter --reset-cache +""" +from __future__ import annotations + +import argparse +import bz2 +import datetime +import hashlib +import os +import pickle +import sys + +import dateutil.parser +import harc_plot +from harc_plot import visualize_histograms as vh + + +NETWORK_DEFAULT_START = { + 'WSPRNet': '2015-01-01', # First WSPRNet data in Madrigal + 'RBN': '2009-02-25', # Full RBN archive span + 'PSKReporter': '2015-01-01', # First PSKReporter data in Madrigal + 'All': '2015-01-01', # Limited by network with latest start +} + + +def _hash_rd(rd): + """Short deterministic hash of the run dict for cache keying.""" + key = repr(sorted((k, repr(v)) for k, v in rd.items())).encode() + return hashlib.sha1(key).hexdigest()[:12] + + +def _load_nc_cache(rd, cache_dir, reset_cache=False): + """Load-and-concatenate NetCDFs, caching the xarray result to a pickle.""" + h = _hash_rd(rd) + fname = f"{rd['sTime']:%Y%m%d}-{rd['eTime']:%Y%m%d}.{h}.p.bz2" + fpath = os.path.join(cache_dir, fname) + os.makedirs(cache_dir, exist_ok=True) + + if os.path.exists(fpath) and reset_cache: + os.remove(fpath) + + if os.path.exists(fpath): + print(f'Loading cached concatenation: {fpath}', file=sys.stderr) + with bz2.BZ2File(fpath, 'rb') as fh: + return pickle.load(fh) + + print(f'Building concatenation (no cache at {fpath})', file=sys.stderr) + nc = vh.ncLoader(**rd) + nc.rd_original = rd + with bz2.BZ2File(fpath, 'wb') as fh: + pickle.dump(nc, fh) + print(f'Cached to {fpath}', file=sys.stderr) + return nc + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument('--network', choices=list(NETWORK_DEFAULT_START), + required=True) + parser.add_argument('--start', default=None, + help='Inclusive start YYYY-MM-DD (default: network-dependent)') + parser.add_argument('--stop', default='2026-03-18', + help='Exclusive stop YYYY-MM-DD (default: 2026-03-18)') + parser.add_argument('--region', default='World') + parser.add_argument('--xb-min', type=float, default=60.0) + parser.add_argument('--yb-km', type=float, default=250.0) + parser.add_argument('--rgc-min', type=float, default=0.0) + parser.add_argument('--rgc-max', type=float, default=10000.0) + parser.add_argument('--data-root', default='/scratch/w2naf/harc_plot_solarcycle', + help='Parent dir containing per-network config subfolders') + parser.add_argument('--plot-root', default='/scratch/w2naf/harc_plot_solarcycle/plots', + help='Output dir for plots') + parser.add_argument('--bands', nargs='+', type=int, + default=[28, 21, 14, 7, 3, 1], + help='Band keys (int MHz from BandData)') + parser.add_argument('--xkey', default='slt_mid', choices=['slt_mid', 'ut_hrs']) + parser.add_argument('--reset-cache', action='store_true') + args = parser.parse_args(argv) + + if args.start is None: + args.start = NETWORK_DEFAULT_START[args.network] + + sTime = dateutil.parser.isoparse(args.start) + eTime = dateutil.parser.isoparse(args.stop) + + run_name = (f'{args.region}-{args.network}' + f'-{int(args.xb_min)}min-{int(args.yb_km)}km' + f'-{int(args.rgc_min)}-{int(args.rgc_max)}km') + data_dir = os.path.join(args.data_root, run_name) + plot_dir = os.path.join(args.plot_root, run_name) + cache_dir = os.path.join(args.plot_root, 'cache') + + if not os.path.isdir(data_dir): + parser.error(f'input dir not found: {data_dir} — did the batch run complete?') + + # ncLoader matches filenames against the date range, so giving it the full + # glob is fine; out-of-range days are simply ignored. + srcs = os.path.join(data_dir, '*.data.nc.bz2') + + rd = { + 'srcs': srcs, + 'sTime': sTime, + 'eTime': eTime, + 'band_keys': args.bands, + 'xkeys': [args.xkey], + 'baseout_dir': plot_dir, + 'plot_region': args.region, + 'geospace_env': None, + 'plot_sza': False, + 'plot_trend': True, + 'plot_kpsymh': False, # Reference data is 2017-era only; turn back on once #4 is done. + 'plot_goes': False, # Ditto. + 'plot_f107': True, # F10.7 is the main solar-activity axis for this paper. + 'log_z': False, + } + + nc = _load_nc_cache(rd, cache_dir=cache_dir, reset_cache=args.reset_cache) + if nc.datasets is None: + print(f'No data loaded for {run_name} in [{sTime}, {eTime}). ' + f'Check the input dir; aborting.', file=sys.stderr) + return 1 + + paths = nc.plot(**rd) + if paths: + print('Wrote plots:', file=sys.stderr) + for p in paths: + print(f' {p}', file=sys.stderr) + + +if __name__ == '__main__': + sys.exit(main() or 0)