From 2e3b99067ead41da3b780ae9052941f21fcfa47a Mon Sep 17 00:00:00 2001 From: Massimiliano Lincetto Date: Mon, 30 Oct 2023 10:17:17 +0100 Subject: [PATCH 1/8] discontinue dozoom mode for create_plot --- skyreader/plot/plot.py | 60 +++++++++++++++++++++++------------------- 1 file changed, 33 insertions(+), 27 deletions(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index c9411842..71566a27 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -58,13 +58,13 @@ 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 +84,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 +138,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 +160,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 +184,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 +256,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 +265,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 +283,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) From 7303a45c17b382ddebd725aef4107f7c055af230 Mon Sep 17 00:00:00 2001 From: wipacdevbot Date: Mon, 30 Oct 2023 09:18:43 +0000 Subject: [PATCH 2/8] update dependencies*.log files(s) --- dependencies-examples.log | 10 +++++----- dependencies-mypy.log | 12 ++++++------ dependencies-tests.log | 8 ++++---- dependencies.log | 6 +++--- 4 files changed, 18 insertions(+), 18 deletions(-) diff --git a/dependencies-examples.log b/dependencies-examples.log index 86bb52d6..7b8079a9 100644 --- a/dependencies-examples.log +++ b/dependencies-examples.log @@ -8,17 +8,17 @@ 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.1 # via requests contourpy==1.1.1 # via matplotlib -cryptography==41.0.4 +cryptography==41.0.5 # via pyjwt cycler==0.12.1 # via matplotlib @@ -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..c644d60f 100644 --- a/dependencies-mypy.log +++ b/dependencies-mypy.log @@ -8,17 +8,17 @@ 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.1 # via requests contourpy==1.1.1 # via matplotlib -cryptography==41.0.4 +cryptography==41.0.5 # via pyjwt cycler==0.12.1 # via matplotlib @@ -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..095d1607 100644 --- a/dependencies-tests.log +++ b/dependencies-tests.log @@ -10,7 +10,7 @@ astropy==5.3.4 # icecube-skyreader (setup.py) certifi==2023.7.22 # via requests -charset-normalizer==3.3.0 +charset-normalizer==3.3.1 # via requests contourpy==1.1.1 # via matplotlib @@ -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..c8af1d50 100644 --- a/dependencies.log +++ b/dependencies.log @@ -10,7 +10,7 @@ astropy==5.3.4 # icecube-skyreader (setup.py) certifi==2023.7.22 # via requests -charset-normalizer==3.3.0 +charset-normalizer==3.3.1 # via requests contourpy==1.1.1 # via matplotlib @@ -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) From b26afbe8b2ed94a97d9ddccf9fee4b07619bd0bf Mon Sep 17 00:00:00 2001 From: Massimiliano Lincetto Date: Fri, 3 Nov 2023 19:45:49 +0100 Subject: [PATCH 3/8] add contour class --- skyreader/plot/contours.py | 68 +++++++++++++++++++++++++++++++ skyreader/plot/plot.py | 82 +++++++++++++++----------------------- 2 files changed, 101 insertions(+), 49 deletions(-) create mode 100644 skyreader/plot/contours.py diff --git a/skyreader/plot/contours.py b/skyreader/plot/contours.py new file mode 100644 index 00000000..c48f65d0 --- /dev/null +++ b/skyreader/plot/contours.py @@ -0,0 +1,68 @@ +from dataclasses import dataclass +import healpy +import numpy as np +from typing import ClassVar + +@dataclass +class GaussianContour: + ra_deg: float + dec_deg: float + radius_50: float + levels: list[float] = [50., 90.] + title: str = "" + color: str = 'm' + marker: str = 'x' + marker_size: int = 20 + + + CONTOUR_STYLES: ClassVar[dict[int, str]] = { 50 : '-', 90: '--' } + SCALE_FACTORS: ClassVar[dict[int, float]] = { 50: 1.177, 90: 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] + + 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.""" + dec = np.pi/2. - self.dec_rad + ra = self.ra_rad + radius_deg = self.sigma_deg * self.sigma2radius(containment=containment) + radius_rad = np.rad2deg(radius_deg) + delta, step, bins = 0, 0, 0 + delta = radius_rad/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 \ No newline at end of file diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 71566a27..48c6549a 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -28,6 +28,8 @@ plot_catalog, ) +from .contours import GaussianContour + from ..result import SkyScanResult LOGGER = logging.getLogger("skyreader.plot") @@ -297,40 +299,15 @@ def create_plot(self, result: SkyScanResult) -> 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): @@ -357,10 +334,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 @@ -638,30 +616,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") From 24e23f75e4c94162fd989d8a5c7d62bb0322ba66 Mon Sep 17 00:00:00 2001 From: wipacdevbot Date: Fri, 3 Nov 2023 18:47:51 +0000 Subject: [PATCH 4/8] update dependencies*.log files(s) --- dependencies-examples.log | 8 ++++---- dependencies-mypy.log | 8 ++++---- dependencies-tests.log | 8 ++++---- dependencies.log | 8 ++++---- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/dependencies-examples.log b/dependencies-examples.log index 7b8079a9..092afdf7 100644 --- a/dependencies-examples.log +++ b/dependencies-examples.log @@ -14,15 +14,15 @@ certifi==2023.7.22 # via requests cffi==1.16.0 # via cryptography -charset-normalizer==3.3.1 +charset-normalizer==3.3.2 # via requests -contourpy==1.1.1 +contourpy==1.2.0 # via matplotlib 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) diff --git a/dependencies-mypy.log b/dependencies-mypy.log index c644d60f..e48aba0c 100644 --- a/dependencies-mypy.log +++ b/dependencies-mypy.log @@ -14,15 +14,15 @@ certifi==2023.7.22 # via requests cffi==1.16.0 # via cryptography -charset-normalizer==3.3.1 +charset-normalizer==3.3.2 # via requests -contourpy==1.1.1 +contourpy==1.2.0 # via matplotlib 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) diff --git a/dependencies-tests.log b/dependencies-tests.log index 095d1607..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.1 +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) diff --git a/dependencies.log b/dependencies.log index c8af1d50..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.1 +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) From b4f51d3e447752ba996e69f59204f32a137afff1 Mon Sep 17 00:00:00 2001 From: Massimiliano Lincetto Date: Wed, 8 Nov 2023 09:33:39 +0100 Subject: [PATCH 5/8] split class --- skyreader/plot/contours.py | 93 +++++++++++++++++++++++--------------- 1 file changed, 57 insertions(+), 36 deletions(-) diff --git a/skyreader/plot/contours.py b/skyreader/plot/contours.py index c48f65d0..0f93a264 100644 --- a/skyreader/plot/contours.py +++ b/skyreader/plot/contours.py @@ -3,36 +3,75 @@ 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: +class GaussianContour(CircularContour): ra_deg: float dec_deg: float radius_50: float - levels: list[float] = [50., 90.] + levels: list[float] = [50.0, 90.0] title: str = "" - color: str = 'm' - marker: str = 'x' + color: str = "m" + marker: str = "x" marker_size: int = 20 - - CONTOUR_STYLES: ClassVar[dict[int, str]] = { 50 : '-', 90: '--' } - SCALE_FACTORS: ClassVar[dict[int, float]] = { 50: 1.177, 90: 2.1459 } + CONTOUR_STYLES: ClassVar[dict[int, str]] = {50: "-", 90: "--"} + SCALE_FACTORS: ClassVar[dict[int, float]] = {50: 1.177, 90: 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] - + 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. @@ -41,28 +80,10 @@ def sigma2radius(self, containment: float): 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.""" - dec = np.pi/2. - self.dec_rad - ra = self.ra_rad - radius_deg = self.sigma_deg * self.sigma2radius(containment=containment) - radius_rad = np.rad2deg(radius_deg) - delta, step, bins = 0, 0, 0 - delta = radius_rad/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 \ No newline at end of file + """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 + ) + return theta, phi From b4532634cc673cb069b1280f18dfd03ca02d042b Mon Sep 17 00:00:00 2001 From: Massimiliano Lincetto Date: Wed, 8 Nov 2023 09:45:10 +0100 Subject: [PATCH 6/8] intermediate step --- skyreader/plot/plot.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 48c6549a..c3b2d4c6 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -28,7 +28,7 @@ plot_catalog, ) -from .contours import GaussianContour +from .contours import GaussianContour, CircularContour from ..result import SkyScanResult @@ -433,8 +433,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] From 160de3007bd229806589daf7a0cb5299c0a5a659 Mon Sep 17 00:00:00 2001 From: Massimiliano Lincetto Date: Wed, 8 Nov 2023 09:57:55 +0100 Subject: [PATCH 7/8] mypy consistency --- skyreader/plot/contours.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/skyreader/plot/contours.py b/skyreader/plot/contours.py index 0f93a264..47b69739 100644 --- a/skyreader/plot/contours.py +++ b/skyreader/plot/contours.py @@ -1,5 +1,5 @@ from dataclasses import dataclass -import healpy +import healpy # type: ignore[import] import numpy as np from typing import ClassVar @@ -54,8 +54,8 @@ class GaussianContour(CircularContour): marker: str = "x" marker_size: int = 20 - CONTOUR_STYLES: ClassVar[dict[int, str]] = {50: "-", 90: "--"} - SCALE_FACTORS: ClassVar[dict[int, float]] = {50: 1.177, 90: 2.1459} + 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): @@ -69,6 +69,10 @@ def ra_rad(self): 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] @@ -84,6 +88,6 @@ def generate_contour(self, nside: int, containment: float): 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 + ra_rad=self.ra_rad, dec_rad=self.dec_rad, radius_rad=radius_rad, nside=nside ) return theta, phi From a3cbd2c66d33279f2a4dbd00f3a327961384b99f Mon Sep 17 00:00:00 2001 From: Massimiliano Lincetto Date: Wed, 8 Nov 2023 14:11:42 +0100 Subject: [PATCH 8/8] fix f-string --- skyreader/plot/plot.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index c3b2d4c6..6ffa4be3 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -62,9 +62,6 @@ def calculate_area(vs) -> float: 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 @@ -601,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}")