diff --git a/pyproject.toml b/pyproject.toml index 4a2265d1..b2812ae3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,6 +10,7 @@ build-backend = "setuptools.build_meta" [project] dependencies = [ + 'PyAstronomy', 'astropy', 'healpy', 'matplotlib', diff --git a/skyreader/constants.py b/skyreader/constants.py index 37c79ef1..e7b94636 100644 --- a/skyreader/constants.py +++ b/skyreader/constants.py @@ -1,5 +1,10 @@ # LAT 14-year Source Catalog (4FGL-DR4) # fermi.gsfc.nasa.gov/ssc/data/access/lat -CATALOG_NAME = "gll_psc_v35.fit" +# CATALOG_NAME = "gll_psc_v35.fit" + +# LAT 16-year Source List (FL16Y) +# https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl16y/ +CATALOG_NAME = "gll_psc_v41.fit" + BASE_DIR = "/cvmfs/icecube.opensciencegrid.org/users/" CATALOG_PATH = f"{BASE_DIR}azegarelli/realtime/catalogs/{CATALOG_NAME}" diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 51dedda3..c0f8c12f 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -26,7 +26,10 @@ plot_catalog ) -# LAT 14-year Source Catalog (4FGL-DR4 in FITS format) ; https://fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/ + +from PyAstronomy import pyasl # type: ignore[import] + +# Fermi catalog from skyreader.constants import CATALOG_PATH from ..utils.areas import calculate_area, get_contour_areas @@ -195,7 +198,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 - # it is actually a valid object before accessing it + # it is actually a valid object before accessing i if isinstance(cb.solids, matplotlib.collections.QuadMesh): cb.solids.set_edgecolor("face") @@ -276,7 +279,7 @@ def create_plot_zoomed( extra_radius=np.nan, systematics=False, plot_bounding_box=False, - plot_4fgl=False, + plot_fermi_sources=False, circular=False, circular_err50=0.2, circular_err90=0.7, @@ -286,10 +289,8 @@ def create_plot_zoomed( catalog_path=CATALOG_PATH ): """Uses healpy to plot a map.""" - def bounding_box(ra, dec, theta, phi): shift = ra-180 - ra_plus = np.max((np.degrees(phi)-shift) % 360) - 180 ra_minus = np.min((np.degrees(phi)-shift) % 360) - 180 dec_plus = (np.pi/2-np.min(theta))*180./np.pi - dec @@ -391,7 +392,7 @@ def bounding_box(ra, dec, theta, phi): contour.T[1] ) - # Find the rough extent of the contours to bound the plot + # Find the rough extent of the contours to bound the plo contours = contours_by_level[-1] ra = min_ra * 180./np.pi dec = min_dec * 180./np.pi @@ -544,9 +545,8 @@ def bounding_box(ra, dec, theta, phi): rot=(lon, lat, 0), bounds=(lower_lon, upper_lon, lower_lat, upper_lat) ) - if plot_4fgl: - # Overlay 4FGL sources - plot_catalog( + if plot_fermi_sources: + ra_fermi_sources, dec_fermi_sources, name_fermi_sources = plot_catalog( equatorial_map, cmap, lower_ra, upper_ra, lower_dec, upper_dec, catalog_path ) @@ -568,17 +568,16 @@ def bounding_box(ra, dec, theta, phi): "dec_minus": dec_minus } # Optional: Print or log the results if needed - contain_txt = "Approximating the " + percentages[l_index] + \ - "% error region as a rectangle, we get:" + " \n" + \ - "\t RA = {0:.2f} + {1:.2f} - {2:.2f}".format( - ra, ra_plus, np.abs(ra_minus)) + " \n" + \ - "\t Dec = {0:.2f} + {1:.2f} - {2:.2f}".format( - dec, dec_plus, np.abs(dec_minus)) + contain_txt = ( + f"Approximating the {percentages[l_index]}% error region as a rectangle, we get:\n" + f"\t RA = {ra:.2f} + {ra_plus:.2f} - {abs(ra_minus):.2f}\n" + f"\t Dec = {dec:.2f} + {dec_plus:.2f} - {dec_minus:.2f}" + ) print(contain_txt) uncertainty = [(ra_minus, ra_plus), (dec_minus, dec_plus)] uncertainties.append(uncertainty) # This is actually an output and not a logging info. - # TODO: we should wrap this in an object, return and log at + # TODO: we should wrap this in an object, return and log a # the higher level. print( f"Contour Area (50%): {contour_areas[0]}", @@ -589,6 +588,9 @@ def bounding_box(ra, dec, theta, phi): "square degrees (scaled)" ) + 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"]) + if plot_bounding_box: bounding_ras_list, bounding_decs_list = [], [] # lower bound @@ -739,6 +741,74 @@ def bounding_box(ra, dec, theta, phi): self._save_contours(contours_by_level, unique_id) LOGGER.info("done.") plt.close() + if plot_fermi_sources: + return sources_in_90 + + 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"] + ra_max = ra_best_fit + rectangular_errors["ra_plus"] + dec_min = dec_best_fit + rectangular_errors["dec_minus"] + dec_max = dec_best_fit + rectangular_errors["dec_plus"] + name_src = np.asarray(name_src) + ra_src = np.asarray(ra_src, dtype=float) + dec_src = np.asarray(dec_src, dtype=float) + # Dec + dec_ok = (dec_src >= dec_min) & (dec_src <= dec_max) + # RA (gestione wrap-around) + ra_min %= 360. + ra_max %= 360. + ra_src = ra_src % 360. + if ra_min <= ra_max: + ra_ok = (ra_src >= ra_min) & (ra_src <= ra_max) + else: + ra_ok = (ra_src >= ra_min) | (ra_src <= ra_max) + mask = ra_ok & dec_ok + return name_src[mask], ra_src[mask], dec_src[mask] + + 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, + dec_best_fit, + ra_src, + dec_src, + name_src, + rectangular_errors + ) + 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% rectangular error box:\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_contours(self, contours_by_level, unique_id) -> None: # Output contours in RA, dec instead of theta, phi diff --git a/skyreader/plot/plotting_tools.py b/skyreader/plot/plotting_tools.py index bbd5c41b..9d9033ae 100644 --- a/skyreader/plot/plotting_tools.py +++ b/skyreader/plot/plotting_tools.py @@ -128,9 +128,9 @@ def hp_ticklabels(zoom=False, lonra=None, latra=None, rot=None, bounds=None): def plot_catalog(master_map, cmap, lower_ra, upper_ra, lower_dec, upper_dec, catalog_path, cmap_min=0., cmap_max=250.): - """"Plots the 4FGL catalog in a color that contrasts with the background + """"Plots the Fermi catalog in a color that contrasts with the background healpix map.""" - hdu = pyfits.open(catalog_path) # LAT 14-year from skyreader.constants or user-specified + hdu = pyfits.open(catalog_path) # Fermi catalog defined in skyreader.constants fgl = hdu[1] pe = [path_effects.Stroke(linewidth=0.5, foreground=cmap(0.0)), path_effects.Normal()] @@ -169,6 +169,8 @@ def color_filter(lon, lat): path_effects=pe) del fgl + return flon_i[fgl_mask]*180./np.pi, flat_i[fgl_mask]*180./np.pi, fname_i[fgl_mask] + ## # Mollweide axes with phi axis flipped and in hours from 24 to 0 instead of # in degrees from -180 to 180. diff --git a/skyreader/utils/handle_map_data.py b/skyreader/utils/handle_map_data.py index f6cee0a4..7815293c 100644 --- a/skyreader/utils/handle_map_data.py +++ b/skyreader/utils/handle_map_data.py @@ -331,9 +331,9 @@ def find_filled_pixels(uniqs: np.ndarray): nside, next_nside, uniqs ) already_filled_uniqs.append(already_filled_uniqs_nside) - already_filled_uniqs = np.concatenate(already_filled_uniqs) + already_filled_uniqs_array = np.concatenate(already_filled_uniqs) already_filled_indeces = np.array( - [np.where(uniqs == uni)[0][0] for uni in already_filled_uniqs] + [np.where(uniqs == uni)[0][0] for uni in already_filled_uniqs_array] ) return already_filled_indeces