Add data normalisation options to get_data #196
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #196 +/- ##
==========================================
+ Coverage 80.97% 81.07% +0.10%
==========================================
Files 35 35
Lines 2391 2410 +19
==========================================
+ Hits 1936 1954 +18
- Misses 455 456 +1 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
There was a problem hiding this comment.
Pull Request Overview
Adds data normalisation options to ScienceData.get_data via a new vtype parameter and threads this through plotting helpers; updates tests and changelog accordingly.
- Introduces vtype ('c', 'cr', 'dcr', 'dcrf') to control count normalisation and units.
- Updates plotting mixins to accept vtype; refactors get_data to compute appropriate normalisation (including geometric area for flux).
- Extends tests to validate returned units for each vtype; adds changelog entry; comments out code in map_reprojection doc examples.
Reviewed Changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 9 comments.
| File | Description |
|---|---|
| stixpy/product/sources/science.py | Adds vtype support in get_data and plotting methods; computes time/energy/area normalisation; updates docstrings. |
| stixpy/product/sources/tests/test_science.py | Adjusts expectations to seconds, passes vtype, adds parametrized unit test for vtype. |
| stixpy/visualisation/map_reprojection.py | Comments example code under .. plot:: blocks, effectively disabling execution in docs. |
| changelog/196.feature.rst | Documents the new vtype options and their units. |
Tip: Customize your code reviews with copilot-instructions.md. Create the file or learn how to get started.
| a_norm = [] | ||
| for pixel_mask in self.data["pixel_masks"]: | ||
| indices = np.nonzero(pixel_mask) | ||
| areas = np.full(12, 0 * u.cm**2) | ||
| areas[indices] = pixel_areas[indices].value |
There was a problem hiding this comment.
a_norm is built from self.data['pixel_masks'] for all original time steps, then reshaped to t_norm.size. If time_indices reduces or merges time bins, len(self.data['pixel_masks']) can differ from t_norm.size, causing reshape/broadcasting errors for vtype='dcrf' or misaligned areas. Slice/aggregate pixel_masks using the same time selection logic as counts/t_norm (e.g., index masks when filtering; for merged ranges, compute a dt-weighted mean or representative area per output bin) before reshaping.
| a_norm = [] | |
| for pixel_mask in self.data["pixel_masks"]: | |
| indices = np.nonzero(pixel_mask) | |
| areas = np.full(12, 0 * u.cm**2) | |
| areas[indices] = pixel_areas[indices].value | |
| # Aggregate pixel masks according to time_indices logic | |
| pixel_masks = self.data["pixel_masks"] | |
| if time_indices is not None: | |
| if time_indices.ndim == 1: | |
| pixel_masks = [pixel_masks[i] for i in time_indices] | |
| elif time_indices.ndim == 2: | |
| # For each merged bin, aggregate pixel masks (dt-weighted mean) | |
| pixel_masks_agg = [] | |
| for idx, (tl, th) in enumerate(time_indices): | |
| # Get time durations for weighting | |
| dt_bin = self.data["timedel"][tl : th + 1].to("s").value | |
| masks_bin = np.array(pixel_masks[tl : th + 1]) | |
| # Weighted mean mask (float, but we want area) | |
| weighted_mask = np.average(masks_bin, axis=0, weights=dt_bin) | |
| # For each pixel, if mask is >0.5, treat as active (or use weighted area) | |
| pixel_masks_agg.append(weighted_mask) | |
| pixel_masks = pixel_masks_agg | |
| elif t_norm.size == 1 and len(pixel_masks) > 1: | |
| # If all times summed, aggregate all masks (dt-weighted mean) | |
| dt_bin = self.data["timedel"].to("s").value | |
| masks_bin = np.array(pixel_masks) | |
| weighted_mask = np.average(masks_bin, axis=0, weights=dt_bin) | |
| pixel_masks = [weighted_mask] | |
| # Build a_norm from pixel_masks | |
| a_norm = [] | |
| for pixel_mask in pixel_masks: | |
| indices = np.nonzero(pixel_mask) | |
| areas = np.full(12, 0 * u.cm**2) | |
| # If mask is float (from weighted mean), use mask value as fraction of area | |
| if np.issubdtype(pixel_mask.dtype, np.floating): | |
| areas = pixel_areas.value * pixel_mask | |
| else: | |
| areas[indices] = pixel_areas[indices].value |
There was a problem hiding this comment.
So this is an interesting case if we sum on ground over a RCR change what area should we return the average or mask it out because we don't know?
There was a problem hiding this comment.
Maybe we could sum over each RCR state separately and then average them in the bin once normalised?
There was a problem hiding this comment.
Potentially but as was discussed at the STIX meeting data at directly before/after are not reliable so I think will set all values for this time to be masked.
| .. plot:: | ||
| :include-source: true | ||
| import astropy.units as u | ||
| from astropy.coordinates import SkyCoord | ||
| from sunpy.net import Fido, attrs as a | ||
| from sunpy.map import Map | ||
| from sunpy.coordinates.frames import HeliographicStonyhurst | ||
| from stixpy.visualisation.map_reprojection import reproject_map, plot_map_reproj, get_solo_position | ||
| # Search and download map using FIDO | ||
| query = Fido.search(a.Time('2020/06/12 13:20:00', '2020/06/12 13:40:00'), | ||
| a.Instrument.hmi, a.Physobs.los_magnetic_field) | ||
| result = Fido.fetch(query[0][0]) | ||
| # Create map and resample to speed up calculations | ||
| map = Map(result) | ||
| map = map.resample([512, 512] * u.pix) | ||
| # Set SOLO as observer | ||
| # observer = get_solo_position(map) # broke on RTD | ||
| observer = SkyCoord(56.18727061*u.deg, 3.5358653*u.deg, 77369481.8542484*u.km, | ||
| frame=HeliographicStonyhurst, obstime='2020-06-12 13:20:08.300000') | ||
| # Reproject Map | ||
| reprojected_map = reproject_map(map, observer, out_shape=(1024, 1024)) | ||
| # Run plotting function | ||
| plot_map_reproj(map, reprojected_map) | ||
| # import astropy.units as u |
There was a problem hiding this comment.
In a .. plot:: directive, the code is executed; commenting out the entire example means nothing runs and may produce empty figures in the docs. If execution is undesired on RTD, replace .. plot:: with .. code-block:: python (or re-enable the code) so the rendered docs are consistent.
|
Apart from the rcr masking, which I believe is not yet included, this LGTM |
2 participants
No description provided.