From 843adffeba0ac39a81397c4aeb4472e7cefb22ca Mon Sep 17 00:00:00 2001 From: Angela Zegarelli <159175351+azegarelli@users.noreply.github.com> Date: Thu, 18 Jun 2026 16:55:22 +0200 Subject: [PATCH 1/8] add functions to check sources within the real 90% contour --- skyreader/plot/plot.py | 73 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 72 insertions(+), 1 deletion(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index c0f8c12..59fe7ad 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -589,7 +589,8 @@ def bounding_box(ra, dec, theta, phi): ) if plot_fermi_sources: - sources_in_90 = self._save_sources_inside_90_box(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, rectangular_errors["90"]) + #sources_in_90 = self._save_sources_inside_90_box(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, rectangular_errors["90"]) + sources_in_90 = self._save_sources_inside_90_contour(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, contours_by_levels[1]) if plot_bounding_box: bounding_ras_list, bounding_decs_list = [], [] @@ -743,6 +744,32 @@ def bounding_box(ra, dec, theta, phi): plt.close() if plot_fermi_sources: return sources_in_90 + + def _sources_inside_90_contour(self, ra_src, dec_src, name_src, contours): + # Convert inputs to numpy arrays + name_src = np.asarray(name_src) + ra_src = np.asarray(ra_src, dtype=float) % 360.0 + dec_src = np.asarray(dec_src, dtype=float) + # Convert (RA, Dec) → spherical coordinates used by healpy + # theta = colatitude (0 = North pole, pi = South pole) + # phi = longitude + src_theta = np.pi/2 - np.deg2rad(dec_src) + src_phi = np.deg2rad(ra_src) + # Combine coordinates into point array: (theta, phi) + points = np.column_stack([src_theta, src_phi]) + # Initialize mask: False for all sources initially + inside_mask = np.zeros(len(points), dtype=bool) + for contour in contours: + # contour is (theta, phi) + contour_theta = contour[:, 0] + contour_phi = contour[:, 1] + # Build polygon in (theta, phi) space + poly = np.column_stack([contour_theta, contour_phi]) + # Create a 2D path object for point-in-polygon test + path = Path(poly) + # Update mask: mark sources inside this contour as True + inside_mask |= path.contains_points(points) + return name_src[inside_mask], ra_src[inside_mask], dec_src[inside_mask] def _sources_inside_90_box(self, ra_best_fit, dec_best_fit, ra_src, dec_src, name_src, rectangular_errors): ra_min = ra_best_fit + rectangular_errors["ra_minus"] @@ -765,6 +792,50 @@ def _sources_inside_90_box(self, ra_best_fit, dec_best_fit, ra_src, dec_src, nam mask = ra_ok & dec_ok return name_src[mask], ra_src[mask], dec_src[mask] + def _save_sources_inside_90_contour(self, ra_best_fit, dec_best_fit, ra_src, dec_src, name_src, contour): + src_inside_90_name, src_inside_90_ra, src_inside_90_dec = self._sources_inside_90_contour( + ra_src, + dec_src, + name_src, + contour + ) + ang_dist = np.zeros(len(src_inside_90_name), dtype=float) + for i in range(len(src_inside_90_name)): + ang_dist[i] = pyasl.getAngDist( + ra_best_fit, + dec_best_fit, + src_inside_90_ra[i], + src_inside_90_dec[i] + ) + sources_in_90 = [ + { + "name": name, + "ra": ra_src, + "dec": dec_src, + "ang_dist": dist + } + for name, ra_src, dec_src, dist in zip( + src_inside_90_name, + src_inside_90_ra, + src_inside_90_dec, + ang_dist + ) + ] + sources_in_90_sorted = sorted( + sources_in_90, + key=lambda s: s["ang_dist"] + ) + print("\nFermi sources inside 90% contour:\n") + for s in sources_in_90_sorted: + print( + f"- {s['name'].strip()}\n" + f" RA: {s['ra']:.3f} deg\n" + f" Dec: {s['dec']:.3f} deg\n" + f" ang_dist: {s['ang_dist']:.3f} deg\n" + ) + return sources_in_90 + + def _save_sources_inside_90_box(self, ra_best_fit, dec_best_fit, ra_src, dec_src, name_src, rectangular_errors): src_inside_90_name, src_inside_90_ra, src_inside_90_dec = self._sources_inside_90_box( ra_best_fit, From 19cf40e1324f56e8a1750e00f107bda02485f1aa Mon Sep 17 00:00:00 2001 From: Angela Zegarelli <159175351+azegarelli@users.noreply.github.com> Date: Thu, 18 Jun 2026 16:58:15 +0200 Subject: [PATCH 2/8] fix flake8 --- skyreader/plot/plot.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 59fe7ad..076612e 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -589,8 +589,8 @@ def bounding_box(ra, dec, theta, phi): ) if plot_fermi_sources: - #sources_in_90 = self._save_sources_inside_90_box(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, rectangular_errors["90"]) - sources_in_90 = self._save_sources_inside_90_contour(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, contours_by_levels[1]) + # sources_in_90 = self._save_sources_inside_90_box(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, rectangular_errors["90"]) + sources_in_90 = self._save_sources_inside_90_contour(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, contours_by_level[1]) if plot_bounding_box: bounding_ras_list, bounding_decs_list = [], [] @@ -835,7 +835,6 @@ def _save_sources_inside_90_contour(self, ra_best_fit, dec_best_fit, ra_src, dec ) return sources_in_90 - def _save_sources_inside_90_box(self, ra_best_fit, dec_best_fit, ra_src, dec_src, name_src, rectangular_errors): src_inside_90_name, src_inside_90_ra, src_inside_90_dec = self._sources_inside_90_box( ra_best_fit, From 712811f6244e6d7fd81bcd07f88b534940fb3231 Mon Sep 17 00:00:00 2001 From: Angela Zegarelli <159175351+azegarelli@users.noreply.github.com> Date: Thu, 18 Jun 2026 17:00:38 +0200 Subject: [PATCH 3/8] fix flake8 From 7c778295d2afb02755a63be6e6b16120eacab654 Mon Sep 17 00:00:00 2001 From: azegarelli Date: Thu, 18 Jun 2026 17:03:31 +0200 Subject: [PATCH 4/8] fix flake8 --- skyreader/plot/plot.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 076612e..a2b949d 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -15,8 +15,8 @@ import numpy as np from astropy.io import ascii # type: ignore[import] from matplotlib import patheffects -from matplotlib import pyplot as plt -from matplotlib import text +from matplotlib import pyplot as pl +from matplotlib import tex from matplotlib.projections import projection_registry # type: ignore[import] from .plotting_tools import ( @@ -39,7 +39,7 @@ prepare_flattened_map, prepare_multiorder_map, ) -from ..result import SkyScanResult +from ..result import SkyScanResul LOGGER = logging.getLogger("skyreader.plot") @@ -52,7 +52,7 @@ class SkyScanPlotter: PLOT_COLORMAP = matplotlib.colormaps['plasma_r'] def __init__(self, output_dir: Path = Path(".")): - # Set here plotting parameters and things that + # Set here plotting parameters and things tha # do not depend on the individual scan. self.output_dir = output_dir projection_registry.register(AstroMollweideAxes) @@ -66,7 +66,7 @@ def create_plot( ) -> None: """Creates a full-sky plot using a meshgrid at fixed resolution. Optionally creates a zoomed-in plot. Resolutions are defined in - PLOT_DPI_STANDARD and PLOT_DPI_ZOOMED. Zoomed mode is very inefficient + PLOT_DPI_STANDARD and PLOT_DPI_ZOOMED. Zoomed mode is very inefficien as the meshgrid is created for the full sky. """ dpi = self.PLOT_DPI_STANDARD @@ -144,7 +144,7 @@ def create_plot( LOGGER.info(f"Preparing plot: {plot_filename}...") # features of the color map to use - cmap.set_under(alpha=0.) # make underflows transparent + cmap.set_under(alpha=0.) # make underflows transparen cmap.set_bad(alpha=1., color=(1., 0., 0.)) # make NaNs bright red # prepare the figure canvas @@ -154,7 +154,7 @@ def create_plot( ax = None - cmap.set_over(alpha=0.) # make underflows transparent + cmap.set_over(alpha=0.) # make underflows transparen ax = fig.add_subplot(111, projection='astro mollweide') # rasterized makes the map bitmap while the labels remain vectorial @@ -197,7 +197,7 @@ def create_plot( # cb.ax.xaxis.labelpad = -8 # workaround for issue with viewers, see colorbar docstring - # mypy compliance: since cb.solids could be None, we check that + # mypy compliance: since cb.solids could be None, we check tha # it is actually a valid object before accessing i if isinstance(cb.solids, matplotlib.collections.QuadMesh): cb.solids.set_edgecolor("face") @@ -209,7 +209,7 @@ def create_plot( ax.grid(True, color='k', alpha=0.5) # Otherwise, add the path effects. - # mypy requires set_path_effects() to take a list of AbstractPathEffect + # mypy requires set_path_effects() to take a list of AbstractPathEffec effects: List[patheffects.AbstractPathEffect] = [ patheffects.withStroke(linewidth=1.1, foreground='w') ] @@ -477,7 +477,7 @@ def bounding_box(ra, dec, theta, phi): orientation='horizontal', aspect=50, ticks=ticks, - format=format + format=forma ) cb.ax.xaxis.set_label_text(cb_label) @@ -634,7 +634,7 @@ def bounding_box(ra, dec, theta, phi): bounding_contour_area = abs(calculate_area(bounding_contour.T)) # convert to square-degrees bounding_contour_area *= (180.*180.)/(np.pi*np.pi) - contour_label = r'90% Bounding rectangle' + \ + contour_label = r'90% Bounding rectangle' + f' - area: {bounding_contour_area:.2f} sqdeg' healpy.projplot( bounding_theta, @@ -744,7 +744,7 @@ def bounding_box(ra, dec, theta, phi): plt.close() if plot_fermi_sources: return sources_in_90 - + def _sources_inside_90_contour(self, ra_src, dec_src, name_src, contours): # Convert inputs to numpy arrays name_src = np.asarray(name_src) @@ -765,7 +765,7 @@ def _sources_inside_90_contour(self, ra_src, dec_src, name_src, contours): contour_phi = contour[:, 1] # Build polygon in (theta, phi) space poly = np.column_stack([contour_theta, contour_phi]) - # Create a 2D path object for point-in-polygon test + # Create a 2D path object for point-in-polygon tes path = Path(poly) # Update mask: mark sources inside this contour as True inside_mask |= path.contains_points(points) @@ -812,13 +812,13 @@ def _save_sources_inside_90_contour(self, ra_best_fit, dec_best_fit, ra_src, dec "name": name, "ra": ra_src, "dec": dec_src, - "ang_dist": dist + "ang_dist": dis } for name, ra_src, dec_src, dist in zip( src_inside_90_name, src_inside_90_ra, src_inside_90_dec, - ang_dist + ang_dis ) ] sources_in_90_sorted = sorted( @@ -834,7 +834,7 @@ def _save_sources_inside_90_contour(self, ra_best_fit, dec_best_fit, ra_src, dec f" ang_dist: {s['ang_dist']:.3f} deg\n" ) return sources_in_90 - + def _save_sources_inside_90_box(self, ra_best_fit, dec_best_fit, ra_src, dec_src, name_src, rectangular_errors): src_inside_90_name, src_inside_90_ra, src_inside_90_dec = self._sources_inside_90_box( ra_best_fit, @@ -857,13 +857,13 @@ def _save_sources_inside_90_box(self, ra_best_fit, dec_best_fit, ra_src, dec_src "name": name, "ra": ra_src, "dec": dec_src, - "ang_dist": dist + "ang_dist": dis } for name, ra_src, dec_src, dist in zip( src_inside_90_name, src_inside_90_ra, src_inside_90_dec, - ang_dist + ang_dis ) ] sources_in_90_sorted = sorted( From 2cf747b7ed60de25540eecef8aa8bba2088758dd Mon Sep 17 00:00:00 2001 From: azegarelli Date: Thu, 18 Jun 2026 17:13:34 +0200 Subject: [PATCH 5/8] restore missing t --- skyreader/plot/plot.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index a2b949d..cd0e092 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -15,8 +15,8 @@ import numpy as np from astropy.io import ascii # type: ignore[import] from matplotlib import patheffects -from matplotlib import pyplot as pl -from matplotlib import tex +from matplotlib import pyplot as plt +from matplotlib import text from matplotlib.projections import projection_registry # type: ignore[import] from .plotting_tools import ( @@ -39,7 +39,7 @@ prepare_flattened_map, prepare_multiorder_map, ) -from ..result import SkyScanResul +from ..result import SkyScanResult LOGGER = logging.getLogger("skyreader.plot") @@ -52,7 +52,7 @@ class SkyScanPlotter: PLOT_COLORMAP = matplotlib.colormaps['plasma_r'] def __init__(self, output_dir: Path = Path(".")): - # Set here plotting parameters and things tha + # Set here plotting parameters and things that # do not depend on the individual scan. self.output_dir = output_dir projection_registry.register(AstroMollweideAxes) @@ -66,7 +66,7 @@ def create_plot( ) -> None: """Creates a full-sky plot using a meshgrid at fixed resolution. Optionally creates a zoomed-in plot. Resolutions are defined in - PLOT_DPI_STANDARD and PLOT_DPI_ZOOMED. Zoomed mode is very inefficien + PLOT_DPI_STANDARD and PLOT_DPI_ZOOMED. Zoomed mode is very inefficient as the meshgrid is created for the full sky. """ dpi = self.PLOT_DPI_STANDARD @@ -144,7 +144,7 @@ def create_plot( LOGGER.info(f"Preparing plot: {plot_filename}...") # features of the color map to use - cmap.set_under(alpha=0.) # make underflows transparen + cmap.set_under(alpha=0.) # make underflows transparent cmap.set_bad(alpha=1., color=(1., 0., 0.)) # make NaNs bright red # prepare the figure canvas @@ -154,7 +154,7 @@ def create_plot( ax = None - cmap.set_over(alpha=0.) # make underflows transparen + cmap.set_over(alpha=0.) # make underflows transparent ax = fig.add_subplot(111, projection='astro mollweide') # rasterized makes the map bitmap while the labels remain vectorial @@ -197,7 +197,7 @@ def create_plot( # cb.ax.xaxis.labelpad = -8 # workaround for issue with viewers, see colorbar docstring - # mypy compliance: since cb.solids could be None, we check tha + # mypy compliance: since cb.solids could be None, we check that # it is actually a valid object before accessing i if isinstance(cb.solids, matplotlib.collections.QuadMesh): cb.solids.set_edgecolor("face") @@ -209,7 +209,7 @@ def create_plot( ax.grid(True, color='k', alpha=0.5) # Otherwise, add the path effects. - # mypy requires set_path_effects() to take a list of AbstractPathEffec + # mypy requires set_path_effects() to take a list of AbstractPathEffect effects: List[patheffects.AbstractPathEffect] = [ patheffects.withStroke(linewidth=1.1, foreground='w') ] @@ -477,7 +477,7 @@ def bounding_box(ra, dec, theta, phi): orientation='horizontal', aspect=50, ticks=ticks, - format=forma + format=format ) cb.ax.xaxis.set_label_text(cb_label) @@ -634,7 +634,7 @@ def bounding_box(ra, dec, theta, phi): bounding_contour_area = abs(calculate_area(bounding_contour.T)) # convert to square-degrees bounding_contour_area *= (180.*180.)/(np.pi*np.pi) - contour_label = r'90% Bounding rectangle' + + contour_label = r'90% Bounding rectangle' + \ f' - area: {bounding_contour_area:.2f} sqdeg' healpy.projplot( bounding_theta, @@ -765,7 +765,7 @@ def _sources_inside_90_contour(self, ra_src, dec_src, name_src, contours): contour_phi = contour[:, 1] # Build polygon in (theta, phi) space poly = np.column_stack([contour_theta, contour_phi]) - # Create a 2D path object for point-in-polygon tes + # Create a 2D path object for point-in-polygon test path = Path(poly) # Update mask: mark sources inside this contour as True inside_mask |= path.contains_points(points) @@ -812,7 +812,7 @@ def _save_sources_inside_90_contour(self, ra_best_fit, dec_best_fit, ra_src, dec "name": name, "ra": ra_src, "dec": dec_src, - "ang_dist": dis + "ang_dist": dist } for name, ra_src, dec_src, dist in zip( src_inside_90_name, @@ -857,13 +857,13 @@ def _save_sources_inside_90_box(self, ra_best_fit, dec_best_fit, ra_src, dec_src "name": name, "ra": ra_src, "dec": dec_src, - "ang_dist": dis + "ang_dist": dist } for name, ra_src, dec_src, dist in zip( src_inside_90_name, src_inside_90_ra, src_inside_90_dec, - ang_dis + ang_dist ) ] sources_in_90_sorted = sorted( From c13d3758e824fdc1209f1761053dbde899aa98c7 Mon Sep 17 00:00:00 2001 From: azegarelli Date: Thu, 18 Jun 2026 17:14:54 +0200 Subject: [PATCH 6/8] restore missing t --- skyreader/plot/plot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index cd0e092..678192f 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -818,7 +818,7 @@ def _save_sources_inside_90_contour(self, ra_best_fit, dec_best_fit, ra_src, dec src_inside_90_name, src_inside_90_ra, src_inside_90_dec, - ang_dis + ang_dist ) ] sources_in_90_sorted = sorted( From 1f9aa0db22da2bf870d423ba304d5b05eccf53e9 Mon Sep 17 00:00:00 2001 From: Angela Zegarelli <159175351+azegarelli@users.noreply.github.com> Date: Thu, 18 Jun 2026 17:23:48 +0200 Subject: [PATCH 7/8] Replace Path with PolyPath for contour polygon --- skyreader/plot/plot.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 678192f..1983895 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -17,6 +17,7 @@ from matplotlib import patheffects from matplotlib import pyplot as plt from matplotlib import text +from matplotlib.path import Path as PolyPath from matplotlib.projections import projection_registry # type: ignore[import] from .plotting_tools import ( @@ -766,7 +767,7 @@ def _sources_inside_90_contour(self, ra_src, dec_src, name_src, contours): # Build polygon in (theta, phi) space poly = np.column_stack([contour_theta, contour_phi]) # Create a 2D path object for point-in-polygon test - path = Path(poly) + path = PolyPath(poly) # Update mask: mark sources inside this contour as True inside_mask |= path.contains_points(points) return name_src[inside_mask], ra_src[inside_mask], dec_src[inside_mask] From 11b5524f9c40c3e11f3d13025187b9475a060264 Mon Sep 17 00:00:00 2001 From: Angela Zegarelli <159175351+azegarelli@users.noreply.github.com> Date: Thu, 18 Jun 2026 17:59:31 +0200 Subject: [PATCH 8/8] Add check for Fermi sources box in plot function --- skyreader/plot/plot.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 1983895..5f0b032 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -281,6 +281,7 @@ def create_plot_zoomed( systematics=False, plot_bounding_box=False, plot_fermi_sources=False, + check_fermi_sources_box=False, circular=False, circular_err50=0.2, circular_err90=0.7, @@ -590,8 +591,10 @@ def bounding_box(ra, dec, theta, phi): ) if plot_fermi_sources: - # sources_in_90 = self._save_sources_inside_90_box(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, rectangular_errors["90"]) - sources_in_90 = self._save_sources_inside_90_contour(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, contours_by_level[1]) + if check_fermi_sources_box: + sources_in_90 = self._save_sources_inside_90_box(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, rectangular_errors["90"]) + else: + sources_in_90 = self._save_sources_inside_90_contour(ra, dec, ra_fermi_sources, dec_fermi_sources, name_fermi_sources, contours_by_level[1]) if plot_bounding_box: bounding_ras_list, bounding_decs_list = [], []