diff --git a/stixpy/imaging/tools.py b/stixpy/imaging/tools.py new file mode 100644 index 00000000..191d43c9 --- /dev/null +++ b/stixpy/imaging/tools.py @@ -0,0 +1,178 @@ +import astropy.units as u +import sunpy.map +import sunpy.time +import xrayvision.imaging +from astropy.time import Time +from sunpy.coordinates import HeliographicStonyhurst +from xrayvision.visibility import Visibilities + +import stixpy.calibration.visibility +import stixpy.coordinates.transforms +from stixpy.coordinates.transforms import STIXImaging + + +def estimate_stix_flare_location(cpd_sci, time_range=None, energy_range=None, subcollimators=None, + cpd_bkg=None, time_range_bkg=None): + """ + Estimates flare location from STIX compressed pixel data. + + FLare location is assumed to be the location of the brightest pixel in the backprojection + map calculated from the input STIX pixel data. + + Parameters + ---------- + cpd_sci: `stixpy.product.Product` + The STIX pixel data. Assumed to be already background subtracted. + time_range: `sunpy.time.TimeRange` (optional) + The time range over which to estimate the flare location. + Default is all times in cpd_sci. + energy_range: `astropy.units.Quantity` in spectral units (optional) + Length-2 quantity giving the lower and upper bounds of the energy range to use for imaging. + Default is all finite energies in cpd_sci. + subcollimators: `iterable` of `str` (optional) + The labels of the subcollimators to include in the output Visibilities object, e.g. + ``['10a', '10b', '10c',...]`` + Default is all subcollimators in cpd_sci + cpd_bkg: stixpy.product.Product` (optional) + The background to subtract from the pixel data before determining the flare location. + Passed to construct_stix_calibrated_visibilities(). + time_range_bkg: `sunpy.time.TimeRange` + The time range within cpd_bkg to use for the background subtraction. + Passed to construct_stix_calibrated_visibilities(). + + Returns + ------- + flare_loc: `astropy.coordinates.SkyCoord` + The estimated flare location. + map_bp: `sunpy.map.Map` + The backprojection map from whose brightest pixel the flare location was estimated. + """ + if time_range is None: + time_range = cpd_sci.time_range + if subcollimators is None: + subcollimators = ["10a", "10b", "10c", "9a", "9b", "9c", + "8a", "8b", "8c", "7a", "7b", "7c"] + # Construct STIX location and centre of FOV. + roll, solo_xyz, pointing = stixpy.coordinates.transforms.get_hpc_info(time_range.start, + time_range.end) + solo = HeliographicStonyhurst(*solo_xyz, obstime=time_range.center, + representation_type="cartesian") + fov_centre = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=time_range.start, + obstime_end=time_range.end, observer=solo) + # Generate calibrated visibilities using coarser subcollimators. + vis = construct_stix_calibrated_visibilities( + cpd_sci, fov_centre, time_range=time_range, energy_range=energy_range, + subcollimators=subcollimators) + # Produce backprojected image and find brightest pixel. Use this for flare location. + imsize = [512, 512] * u.pixel # number of pixels of the map to reconstruct + plate_scale = [10, 10] * u.arcsec / u.pixel # pixel size in arcsec + bp_image = xrayvision.imaging.vis_to_image(vis, imsize, pixel_size=plate_scale) + max_idx = np.argwhere(bp_image == bp_image.max()).ravel() + # Calculate WCS for backprojected image in STIXImaging and HPC frames. + # Recalculate STIX HPC info as slightly different times will have been used + # than input times due to onboard STIX time binning. + vis_tr = sunpy.time.TimeRange(vis.meta["time_range"]) + roll, solo_xyz, pointing = stixpy.coordinates.transforms.get_hpc_info(vis_tr.start, vis_tr.end) + solo = HeliographicStonyhurst(*solo_xyz, obstime=vis_tr.center, representation_type="cartesian") + coord = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=vis_tr.start, obstime_end=vis_tr.end, observer=solo) + header_bp = sunpy.map.make_fitswcs_header(bp_image, coord, telescope="STIX", + observatory="Solar Orbiter", scale=plate_scale) + map_bp = sunpy.map.Map(bp_image, header_bp) + wcs_bp = map_bp.wcs + # Estimate flare location from brightest pixel in backprojection image + flare_loc = wcs_bp.array_index_to_world(*max_idx) + flare_loc = STIXImaging(flare_loc[0].to(u.arcsec), flare_loc[1].to(u.arcsec), + obstime=vis_tr.start, obstime_end=vis_tr.end, observer=solo) + return flare_loc, map_bp + + +def construct_stix_calibrated_visibilities(cpd_sci, flare_location, time_range=None, + energy_range=None, subcollimators=None, + cpd_bkg=None, time_range_bkg=None, **kwargs): + """ + Constructs calibrated STIX visibilities from STIX compressed pixel data. + + Extra kwargs are passed to `stixpy.calibration.visibility.create_meta_pixels` + + Parameters + ---------- + cpd_sci: `stixpy.product.Product` + The STIX pixel data. Assumed to be already background subtracted. + flare_location: `astropy.coordinates.SkyCoord` + The flare location. Frame must be convertible to `stixpy.coordinates.transdforms.STIXImaging`. + time_range: `sunpy.time.TimeRange` (optional) + The time range over which to estimate the flare location. + Default is all times in cpd_sci. + energy_range: `astropy.units.Quantity` in spectral units (optional) + Length-2 quantity giving the lower and upper bounds of the energy range to use for imaging. + Default is all finite energies in cpd_sci. + subcollimators: `iterable` of `str` + The labels of the subcollimators to use in estimating the flare locations, e.g. + ``['10a', '10b', '10c',...]`` + Default is all subcollimators in cpd_sci. + cpd_bkg: stixpy.product.Product` (optional) + The background to subtract from the pixel data before determining the flare location. + If None and time_range_bkg is also None, no background is subtracted. + time_range_bkg: `sunpy.time.TimeRange` + The time range within cpd_bkg to use for the background subtraction. + If None, entire time range of cpd_bkg is used to determine background. + If not None, and cpd_bkg is None, the background is determined from this time range + applied to cpd_sci. + + Returns + ------- + vis: `xrayvision.visibility.Visibilities` + The calibrated STIX visibilities. + """ + # Sanitze inputs. + if time_range is None: + time_range = cpd_sci.time_range + times = Time([time_range.start, time_range.end]) + if energy_range is None: + energy_range = u.Quantity([cpd_sci.energies["e_low"][np.isfinite(cpd_sci.energies["e_low"])][0], + cpd_sci.energies["e_high"][np.isfinite(cpd_sci.energies["e_high"])][-1]]) + no_shadowing = kwargs.pop("no_shadowing", True) + # Generate meta_pixels from pixel_data. + meta_pixels = stixpy.calibration.visibility.create_meta_pixels( + cpd_sci, time_range=times, energy_range=energy_range, + flare_location=flare_location, no_shadowing=no_shadowing, **kwargs) + # Subtract background if background pixel data provided. + if cpd_bkg is None and time_range_bkg is not None: + cpd_bkg = cpd_sci + if cpd_bkg is not None: + if time_range_bkg is None: + time_range_bkg = cpd_bkg.time_range + times_bkg = Time([time_range_bkg.start, time_range_bkg.end]) + meta_pixels_bkg = stixpy.calibration.visibility.create_meta_pixels( + cpd_bkg, time_range=times_bkg, energy_range=energy_range, + flare_location=[0, 0] * u.arcsec, no_shadowing=no_shadowing, **kwargs) + meta_pixels = _subtract_background_from_stix_meta_pixels(meta_pixels, meta_pixels_bkg) + # Generate and calibrate visibilities. + vis = stixpy.calibration.visibility.create_visibility(meta_pixels) + vis = stixpy.calibration.visibility.calibrate_visibility(vis, flare_location=flare_location) + if subcollimators is not None: + idx_subcol = np.argwhere(np.isin(vis.meta["vis_labels"], subcollimator_labels)).ravel() + vis = vis[idx_subcol] + return vis + + +def _subtract_background_from_stix_meta_pixels(meta_pixels_sci, meta_pixels_bkg): + """ + Estimates flare location from STIX pixel data. + + Parameters + ---------- + meta_pixels_sci: `dict` + The STIX meta pixel representing the observations. Format must be + same as output from `stixpy.calibration.visibility.create_meta_pixels`. + meta_pixels_bkg: `dict` + The STIX meta pixels representing the background. Format must be + same as output from `stixpy.calibration.visibility.create_meta_pixels`. + """ + meta_pixels_bkg_sub = { + **meta_pixels_sci, + "abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"], + "abcd_rate_error_kev_cm": np.sqrt(meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2), + } + return meta_pixels_bkg_sub diff --git a/stixpy/map/stix.py b/stixpy/map/stix.py index 9f9ecfaa..0029a991 100644 --- a/stixpy/map/stix.py +++ b/stixpy/map/stix.py @@ -1,7 +1,10 @@ import astropy.units as u from astropy import wcs +from sunpy.coordinates import HeliographicStonyhurst from sunpy.map import GenericMap +import stixpy.coordinates.transforms + __all__ = ["STIXMap"] @@ -72,3 +75,45 @@ def wcs(self): # Validate the WCS here. w2.wcs.set() return w2 + + def to_hpc(self, + reference_coordinate, + reference_pixel: u.Unit('pix') = None): + """ + Return a version of the map in the Helioprojective Cartesian (HPC) coordinate frame. + + This is quicker than a full reprojection, `~stixpy.coordinates.transforms.STIXImaging` + is also a helioprojective frame with a different z-axis and rotation. + Therefore, the new WCS information can be unambiguously reconstructed without + altering the data array. + + Parameters + ---------- + reference_coordinate: `astropy.coordinates.SkyCoord` + The coordinate of the reference pixel. + Must be transformable to `sunpy.coordinates.Helioprojective`. + reference_pixel: `astropy.coordinate.Quantity` length-2 (optional) + The (x, y) pixel index of the reference pixel. + Default is center of the map's field of view. + + Returns + ------- + hpc_map: `sunpy.map.Map` + Map of the STIX image in the HPC coordinate frame. + """ + if reference_pixel is None: + reference_pixel = (np.array(self.data.shape)[::-1] / 2) << u.pix + roll, solo_xyz, pointing = stixpy.coordinates.transforms.get_hpc_info( + reference_coordinate.obstime) + solo = HeliographicStonyhurst(*solo_xyz, obstime=reference_coordinate.obstime, + representation_type="cartesian") + header = make_fitswcs_header( + self.data, + reference_coordinate.transform_to(Helioprojective(obstime=reference_coordinate.obstime, + observer=solo)), + telescope="STIX", + observatory="Solar Orbiter", + scale=u.Quantity(self.scale), + rotation_angle=90 * u.deg + roll,) + return sunpy.map.Map(self.data, header, + mask=self.mask, uncertainty=self.uncertainty, meta=self.meta)