diff --git a/dependencies-examples.log b/dependencies-examples.log index 86bb52d6..092afdf7 100644 --- a/dependencies-examples.log +++ b/dependencies-examples.log @@ -8,21 +8,21 @@ astropy==5.3.4 # via # healpy # icecube-skyreader (setup.py) -cachetools==5.3.1 +cachetools==5.3.2 # via wipac-rest-tools certifi==2023.7.22 # via requests cffi==1.16.0 # via cryptography -charset-normalizer==3.3.0 +charset-normalizer==3.3.2 # via requests -contourpy==1.1.1 +contourpy==1.2.0 # via matplotlib -cryptography==41.0.4 +cryptography==41.0.5 # via pyjwt cycler==0.12.1 # via matplotlib -fonttools==4.43.1 +fonttools==4.44.0 # via matplotlib healpy==1.16.6 # via icecube-skyreader (setup.py) @@ -30,7 +30,7 @@ idna==3.4 # via requests kiwisolver==1.4.5 # via matplotlib -matplotlib==3.8.0 +matplotlib==3.8.1 # via # healpy # icecube-skyreader (setup.py) @@ -50,7 +50,7 @@ packaging==23.2 # via # astropy # matplotlib -pandas==2.1.1 +pandas==2.1.2 # via icecube-skyreader (setup.py) pillow==10.1.0 # via matplotlib @@ -99,7 +99,7 @@ urllib3==2.0.7 # via # requests # wipac-rest-tools -wipac-dev-tools==1.7.0 +wipac-dev-tools==1.7.1 # via # icecube-skyreader (setup.py) # wipac-rest-tools diff --git a/dependencies-mypy.log b/dependencies-mypy.log index 8c46f0ec..e48aba0c 100644 --- a/dependencies-mypy.log +++ b/dependencies-mypy.log @@ -8,21 +8,21 @@ astropy==5.3.4 # via # healpy # icecube-skyreader (setup.py) -cachetools==5.3.1 +cachetools==5.3.2 # via wipac-rest-tools certifi==2023.7.22 # via requests cffi==1.16.0 # via cryptography -charset-normalizer==3.3.0 +charset-normalizer==3.3.2 # via requests -contourpy==1.1.1 +contourpy==1.2.0 # via matplotlib -cryptography==41.0.4 +cryptography==41.0.5 # via pyjwt cycler==0.12.1 # via matplotlib -fonttools==4.43.1 +fonttools==4.44.0 # via matplotlib healpy==1.16.6 # via icecube-skyreader (setup.py) @@ -32,7 +32,7 @@ iniconfig==2.0.0 # via pytest kiwisolver==1.4.5 # via matplotlib -matplotlib==3.8.0 +matplotlib==3.8.1 # via # healpy # icecube-skyreader (setup.py) @@ -53,7 +53,7 @@ packaging==23.2 # astropy # matplotlib # pytest -pandas==2.1.1 +pandas==2.1.2 # via icecube-skyreader (setup.py) pillow==10.1.0 # via matplotlib @@ -71,7 +71,7 @@ pyparsing==3.1.1 # via matplotlib pypng==0.20220715.0 # via qrcode -pytest==7.4.2 +pytest==7.4.3 # via # icecube-skyreader (setup.py) # pytest-mock @@ -110,7 +110,7 @@ urllib3==2.0.7 # via # requests # wipac-rest-tools -wipac-dev-tools==1.7.0 +wipac-dev-tools==1.7.1 # via # icecube-skyreader (setup.py) # wipac-rest-tools diff --git a/dependencies-tests.log b/dependencies-tests.log index aaef8426..d83da2f9 100644 --- a/dependencies-tests.log +++ b/dependencies-tests.log @@ -10,13 +10,13 @@ astropy==5.3.4 # icecube-skyreader (setup.py) certifi==2023.7.22 # via requests -charset-normalizer==3.3.0 +charset-normalizer==3.3.2 # via requests -contourpy==1.1.1 +contourpy==1.2.0 # via matplotlib cycler==0.12.1 # via matplotlib -fonttools==4.43.1 +fonttools==4.44.0 # via matplotlib healpy==1.16.6 # via icecube-skyreader (setup.py) @@ -26,7 +26,7 @@ iniconfig==2.0.0 # via pytest kiwisolver==1.4.5 # via matplotlib -matplotlib==3.8.0 +matplotlib==3.8.1 # via # healpy # icecube-skyreader (setup.py) @@ -47,7 +47,7 @@ packaging==23.2 # astropy # matplotlib # pytest -pandas==2.1.1 +pandas==2.1.2 # via icecube-skyreader (setup.py) pillow==10.1.0 # via matplotlib @@ -57,7 +57,7 @@ pyerfa==2.0.1.1 # via astropy pyparsing==3.1.1 # via matplotlib -pytest==7.4.2 +pytest==7.4.3 # via # icecube-skyreader (setup.py) # pytest-mock @@ -83,5 +83,5 @@ tzdata==2023.3 # via pandas urllib3==2.0.7 # via requests -wipac-dev-tools==1.7.0 +wipac-dev-tools==1.7.1 # via icecube-skyreader (setup.py) diff --git a/dependencies.log b/dependencies.log index 467595e5..d8112a8d 100644 --- a/dependencies.log +++ b/dependencies.log @@ -10,13 +10,13 @@ astropy==5.3.4 # icecube-skyreader (setup.py) certifi==2023.7.22 # via requests -charset-normalizer==3.3.0 +charset-normalizer==3.3.2 # via requests -contourpy==1.1.1 +contourpy==1.2.0 # via matplotlib cycler==0.12.1 # via matplotlib -fonttools==4.43.1 +fonttools==4.44.0 # via matplotlib healpy==1.16.6 # via icecube-skyreader (setup.py) @@ -24,7 +24,7 @@ idna==3.4 # via requests kiwisolver==1.4.5 # via matplotlib -matplotlib==3.8.0 +matplotlib==3.8.1 # via # healpy # icecube-skyreader (setup.py) @@ -44,7 +44,7 @@ packaging==23.2 # via # astropy # matplotlib -pandas==2.1.1 +pandas==2.1.2 # via icecube-skyreader (setup.py) pillow==10.1.0 # via matplotlib @@ -72,5 +72,5 @@ tzdata==2023.3 # via pandas urllib3==2.0.7 # via requests -wipac-dev-tools==1.7.0 +wipac-dev-tools==1.7.1 # via icecube-skyreader (setup.py) diff --git a/skyreader/plot/contours.py b/skyreader/plot/contours.py new file mode 100644 index 00000000..47b69739 --- /dev/null +++ b/skyreader/plot/contours.py @@ -0,0 +1,93 @@ +from dataclasses import dataclass +import healpy # type: ignore[import] +import numpy as np +from typing import ClassVar + + +class CircularContour: + @staticmethod + def circular_contour(ra_rad: float, dec_rad: float, radius_rad: float, nside: int): + """For plotting circular contours on skymaps ra, dec, sigma all + expected in radians.""" + # This transformation is likely reversed later. It can probably be removed. + # Likely legacy of old Skymap Scanner design. + dec_rad = np.pi / 2.0 - dec_rad + step = 1.0 / np.sin(radius_rad) / 10.0 + # bins * step = 360 deg + bins = int(360.0 / step) + Theta = np.zeros(bins + 1, dtype=np.double) + Phi = np.zeros(bins + 1, dtype=np.double) + # Define the contour + for j in range(0, bins): + # phi runs over 0, 2 pi + phi = j * step / 180.0 * np.pi + vx = np.cos(phi) * np.sin(ra_rad) * np.sin(radius_rad) + np.cos(ra_rad) * ( + np.cos(radius_rad) * np.sin(dec_rad) + + np.cos(dec_rad) * np.sin(radius_rad) * np.sin(phi) + ) + vy = np.cos(radius_rad) * np.sin(dec_rad) * np.sin(ra_rad) + np.sin( + radius_rad + ) * ( + -np.cos(ra_rad) * np.cos(phi) + + np.cos(dec_rad) * np.sin(ra_rad) * np.sin(phi) + ) + vz = np.cos(dec_rad) * np.cos(radius_rad) - np.sin(dec_rad) * np.sin( + radius_rad + ) * np.sin(phi) + idx = healpy.vec2pix(nside, vx, vy, vz) + Theta[j], Phi[j] = healpy.pix2ang(nside, idx) + + # Close the contour + Theta[bins], Phi[bins] = Theta[0], Phi[0] + + return Theta, Phi + + +@dataclass +class GaussianContour(CircularContour): + ra_deg: float + dec_deg: float + radius_50: float + levels: list[float] = [50.0, 90.0] + title: str = "" + color: str = "m" + marker: str = "x" + marker_size: int = 20 + + CONTOUR_STYLES: ClassVar[dict[float, str]] = {50.0: "-", 90.0: "--"} + SCALE_FACTORS: ClassVar[dict[float, float]] = {50.0: 1.177, 90.0: 2.1459} + + @property + def dec_rad(self): + return np.deg2rad(self.dec_deg) + + @property + def ra_rad(self): + return np.deg2rad(self.ra_deg) + + @property + def sigma_deg(self): + return self.radius_50 / self.SCALE_FACTORS[50] + + @property + def sigma_rad(self): + return np.rad2deg(self.sigma_deg) + + def get_style(self, containment: float): + return self.CONTOUR_STYLES[containment] + + def sigma2radius(self, containment: float): + # Converts the sigma (standard deviation) to a given containment radius + # given a required % of containment. + # We use hardcoded values but this can be dynamically calculated with + # a chi2 PDF. + return self.SCALE_FACTORS[containment] + + def generate_contour(self, nside: int, containment: float): + """For plotting circular contours on skymaps ra, dec, sigma all + expected in radians.""" + radius_rad = self.sigma_rad * self.sigma2radius(containment=containment) + theta, phi = self.circular_contour( + ra_rad=self.ra_rad, dec_rad=self.dec_rad, radius_rad=radius_rad, nside=nside + ) + return theta, phi diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index c9411842..6ffa4be3 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -28,6 +28,8 @@ plot_catalog, ) +from .contours import GaussianContour, CircularContour + from ..result import SkyScanResult LOGGER = logging.getLogger("skyreader.plot") @@ -58,13 +60,10 @@ def calculate_area(vs) -> float: y0 = y1 return a - def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: + def create_plot(self, result: SkyScanResult) -> 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 - as the meshgrid is created for the full sky. """ - dpi = self.PLOT_DPI_STANDARD if not dozoom else self.PLOT_DPI_ZOOMED + dpi = self.PLOT_DPI_STANDARD # if not dozoom else self.PLOT_DPI_ZOOMED xsize = int(self.PLOT_SIZE_X_IN * dpi) # number of grid points along RA coordinate ysize = int(xsize // 2) # number of grid points along dec coordinate @@ -84,7 +83,8 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: unique_id = f'{str(event_metadata)}_{result.get_nside_string()}' plot_title = f"Run: {event_metadata.run_id} Event {event_metadata.event_id}: Type: {event_metadata.event_type} MJD: {event_metadata.mjd}" - plot_filename = f"{unique_id}.{'plot_zoomed_legacy.' if dozoom else ''}pdf" + # plot_filename = f"{unique_id}.{'plot_zoomed_legacy.' if dozoom else ''}pdf" + plot_filename = f"{unique_id}.pdf" LOGGER.info(f"saving plot to {plot_filename}") min_llh, max_llh = np.nan, np.nan @@ -137,11 +137,13 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: LOGGER.info(f"min Dec: {min_dec * 180./np.pi} deg") # renormalize - if dozoom: - grid_map = grid_map - min_llh - # max_value = max_value - min_value - min_llh = 0. - max_llh = 50 + # plot only delta LLH for zoomed plot + # crosscheck with create_plot_zoomed + # if dozoom: + # grid_map = grid_map - min_llh + # # max_value = max_value - min_value + # min_llh = 0. + # max_llh = 50 grid_map = np.ma.masked_invalid(grid_map) @@ -157,11 +159,11 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: ax = None - if dozoom: - ax = fig.add_subplot(111) #,projection='cartesian') - else: - cmap.set_over(alpha=0.) # make underflows transparent - ax = fig.add_subplot(111,projection='astro mollweide') + # if dozoom: + # ax = fig.add_subplot(111) #,projection='cartesian') + # else: + 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 # flip longitude to the astro convention @@ -181,14 +183,16 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: e, _ = contour_set.legend_elements() leg_element.append(e[0]) - if not dozoom: - # graticule - if isinstance(ax, AstroMollweideAxes): - # mypy guard - ax.set_longitude_grid(30) - ax.set_latitude_grid(30) - cb = fig.colorbar(image, orientation='horizontal', shrink=.6, pad=0.05, ticks=[min_llh, max_llh]) - cb.ax.xaxis.set_label_text(r"$-2 \ln(L)$") + # if not dozoom: + # graticule + if isinstance(ax, AstroMollweideAxes): + # mypy guard + ax.set_longitude_grid(30) + ax.set_latitude_grid(30) + cb = fig.colorbar(image, orientation='horizontal', shrink=.6, pad=0.05, ticks=[min_llh, max_llh]) + cb.ax.xaxis.set_label_text(r"$-2 \ln(L)$") + + """ else: ax.set_xlabel('right ascension') ax.set_ylabel('declination') @@ -251,6 +255,7 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: ax.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(base=tick_label_grid)) ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(base=tick_label_grid)) + """ # cb.ax.xaxis.labelpad = -8 # workaround for issue with viewers, see colorbar docstring @@ -259,8 +264,8 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: if isinstance(cb.solids, matplotlib.collections.QuadMesh): cb.solids.set_edgecolor("face") - if dozoom: - ax.set_aspect('equal') + # if dozoom: + # ax.set_aspect('equal') ax.tick_params(axis='x', labelsize=10) ax.tick_params(axis='y', labelsize=10) @@ -277,10 +282,10 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: # remove white space around figure spacing = 0.01 - if not dozoom: - fig.subplots_adjust(bottom=spacing, top=1.-spacing, left=spacing+0.04, right=1.-spacing) - else: - fig.subplots_adjust(bottom=spacing, top=0.92-spacing, left=spacing+0.1, right=1.-spacing) + # if not dozoom: + fig.subplots_adjust(bottom=spacing, top=1.-spacing, left=spacing+0.04, right=1.-spacing) + # else: + # fig.subplots_adjust(bottom=spacing, top=0.92-spacing, left=spacing+0.1, right=1.-spacing) # set the title fig.suptitle(plot_title) @@ -291,40 +296,15 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: LOGGER.info("done.") - @staticmethod - def circular_contour(ra, dec, sigma, nside): - """For plotting circular contours on skymaps ra, dec, sigma all - expected in radians.""" - dec = np.pi/2. - dec - sigma = np.rad2deg(sigma) - delta, step, bins = 0, 0, 0 - delta= sigma/180.0*np.pi - step = 1./np.sin(delta)/10. - bins = int(360./step) - Theta = np.zeros(bins+1, dtype=np.double) - Phi = np.zeros(bins+1, dtype=np.double) - # define the contour - for j in range(0,bins) : - phi = j*step/180.*np.pi - vx = np.cos(phi)*np.sin(ra)*np.sin(delta) + np.cos(ra)*(np.cos(delta)*np.sin(dec) + np.cos(dec)*np.sin(delta)*np.sin(phi)) - vy = np.cos(delta)*np.sin(dec)*np.sin(ra) + np.sin(delta)*(-np.cos(ra)*np.cos(phi) + np.cos(dec)*np.sin(ra)*np.sin(phi)) - vz = np.cos(dec)*np.cos(delta) - np.sin(dec)*np.sin(delta)*np.sin(phi) - idx = healpy.vec2pix(nside, vx, vy, vz) - DEC, RA = healpy.pix2ang(nside, idx) - Theta[j] = DEC - Phi[j] = RA - Theta[bins] = Theta[0] - Phi[bins] = Phi[0] - return Theta, Phi - def create_plot_zoomed(self, result: SkyScanResult, + # contours: list[ContourSpec], + plot_bounding_box=False, + plot_4fgl=False, extra_ra=np.nan, extra_dec=np.nan, extra_radius=np.nan, systematics=False, - plot_bounding_box=False, - plot_4fgl=False, circular=False, circular_err50=0.2, circular_err90=0.7): @@ -351,10 +331,11 @@ def bounding_box(ra, dec, theta, phi): nsides = result.nsides LOGGER.info(f"available nsides: {nsides}") - if systematics is not True: - plot_filename = unique_id + ".plot_zoomed_wilks.pdf" - else: - plot_filename = unique_id + ".plot_zoomed.pdf" + # if systematics is not True: + # plot_filename = unique_id + ".plot_zoomed_wilks.pdf" + + plot_filename = unique_id + ".plot_zoomed.pdf" + LOGGER.info(f"saving plot to {plot_filename}") nsides = result.nsides @@ -449,8 +430,8 @@ def bounding_box(ra, dec, theta, phi): if circular: sigma50 = np.deg2rad(circular_err50) sigma90 = np.deg2rad(circular_err90) - Theta50, Phi50 = self.circular_contour(min_ra, min_dec, sigma50, nside) - Theta90, Phi90 = self.circular_contour(min_ra, min_dec, sigma90, nside) + Theta50, Phi50 = CircularContour.circular_contour(min_ra, min_dec, sigma50, nside) + Theta90, Phi90 = CircularContour.circular_contour(min_ra, min_dec, sigma90, nside) contour50 = np.dstack((Theta50,Phi50)) contour90 = np.dstack((Theta90,Phi90)) contours_by_level = [contour50, contour90] @@ -617,7 +598,7 @@ def bounding_box(ra, dec, theta, phi): tab = {"ra (rad)": ras, "dec (rad)": decs} savename = unique_id + ".contour_" + val + ".txt" try: - LOGGER.info("Dumping to {savename}") + LOGGER.info(f"Writing contour to {savename}") ascii.write(tab, savename, overwrite=True) except OSError as err: LOGGER.error("OS Error prevented contours from being written, maybe a memory issue. Error is:\n{err}") @@ -632,30 +613,36 @@ def bounding_box(ra, dec, theta, phi): ) mmap_nside = healpy.get_nside(equatorial_map) - # Plot the original online reconstruction location + # Plot additional contours + + if np.sum(np.isnan([extra_ra, extra_dec, extra_radius])) == 0: # dist = angular_distance(minRA, minDec, extra_ra * np.pi/180., extra_dec * np.pi/180.) # print("Millipede best fit is", dist /(np.pi * extra_radius/(1.177 * 180.)), "sigma from reported best fit") + map_nside = healpy.get_nside(equatorial_map) - extra_ra_rad = np.radians(extra_ra) - extra_dec_rad = np.radians(extra_dec) - extra_radius_rad = np.radians(extra_radius) - extra_lon = extra_ra_rad - extra_lat = extra_dec_rad - - healpy.projscatter(np.degrees(extra_lon), np.degrees(extra_lat), - lonlat=True, c='m', marker='x', s=20, label=r'Reported online (50%, 90%)') - for cont_lev, cont_scale, cont_col, cont_sty in zip(['50', '90.'], - [1., 2.1459/1.177], ['m', 'm'], ['-', '--']): - spline_contour = self.circular_contour(extra_ra_rad, extra_dec_rad, - extra_radius_rad*cont_scale, healpy.get_nside(equatorial_map)) - spline_lon = spline_contour[1] + c = GaussianContour(ra_deg=extra_ra, dec_deg=extra_dec, radius_50=extra_radius, title="Online SplineMPE") + + # Plot center + healpy.projscatter(c.ra_deg, + c.dec_deg, + lonlat=True, + c=c.color, + marker=c.marker, + s=c.marker_size, + label=c.title) + + for level in c.levels: + spline_contour = c.generate_contour(nside=map_nside, containment=level) spline_lat = np.pi/2. - spline_contour[0] - healpy.projplot(np.degrees(spline_lon), np.degrees(spline_lat), - lonlat=True, linewidth=2., color=cont_col, - linestyle=cont_sty) + spline_lon = spline_contour[1] + healpy.projplot(np.degrees(spline_lon), + np.degrees(spline_lat), + lonlat=True, linewidth=2., + color=c.color, + linestyle=c.get_style(containment=level)) plt.legend(fontsize=6, loc="lower left")