diff --git a/.github/workflows/wipac-cicd.yml b/.github/workflows/wipac-cicd.yml index ceb174fc..ecb53664 100644 --- a/.github/workflows/wipac-cicd.yml +++ b/.github/workflows/wipac-cicd.yml @@ -27,7 +27,9 @@ jobs: steps: - uses: actions/checkout@v3 - uses: actions/setup-python@v3 - - uses: WIPACrepo/wipac-dev-flake8-action@v1.0 + - uses: WIPACrepo/wipac-dev-flake8-action@v1.1 + with: + max-function-length: 400 mypy: needs: [py-versions] diff --git a/dependencies-examples.log b/dependencies-examples.log index e0a52aa8..8a950460 100644 --- a/dependencies-examples.log +++ b/dependencies-examples.log @@ -8,7 +8,7 @@ astropy==6.1.4 # via # healpy # icecube-skyreader (setup.py) -astropy-iers-data==0.2024.9.30.0.32.59 +astropy-iers-data==0.2024.10.7.0.32.46 # via astropy cachetools==5.5.0 # via wipac-rest-tools @@ -16,7 +16,7 @@ certifi==2024.8.30 # via requests cffi==1.17.1 # via cryptography -charset-normalizer==3.3.2 +charset-normalizer==3.4.0 # via requests contourpy==1.3.0 # via matplotlib @@ -36,7 +36,7 @@ matplotlib==3.9.2 # via icecube-skyreader (setup.py) meander==0.0.3 # via icecube-skyreader (setup.py) -numpy==2.1.1 +numpy==2.1.2 # via # astropy # contourpy @@ -93,9 +93,9 @@ urllib3==2.2.3 # via # requests # wipac-rest-tools -wipac-dev-tools==1.12.1 +wipac-dev-tools==1.13.0 # via # icecube-skyreader (setup.py) # wipac-rest-tools -wipac-rest-tools==1.7.9 +wipac-rest-tools==1.8.0 # via icecube-skyreader (setup.py) diff --git a/dependencies-mypy.log b/dependencies-mypy.log index f9de1f54..357a35d4 100644 --- a/dependencies-mypy.log +++ b/dependencies-mypy.log @@ -8,7 +8,7 @@ astropy==6.1.4 # via # healpy # icecube-skyreader (setup.py) -astropy-iers-data==0.2024.9.30.0.32.59 +astropy-iers-data==0.2024.10.7.0.32.46 # via astropy cachetools==5.5.0 # via wipac-rest-tools @@ -16,7 +16,7 @@ certifi==2024.8.30 # via requests cffi==1.17.1 # via cryptography -charset-normalizer==3.3.2 +charset-normalizer==3.4.0 # via requests contourpy==1.3.0 # via matplotlib @@ -38,7 +38,7 @@ matplotlib==3.9.2 # via icecube-skyreader (setup.py) meander==0.0.3 # via icecube-skyreader (setup.py) -numpy==2.1.1 +numpy==2.1.2 # via # astropy # contourpy @@ -104,9 +104,9 @@ urllib3==2.2.3 # via # requests # wipac-rest-tools -wipac-dev-tools==1.12.1 +wipac-dev-tools==1.13.0 # via # icecube-skyreader (setup.py) # wipac-rest-tools -wipac-rest-tools==1.7.9 +wipac-rest-tools==1.8.0 # via icecube-skyreader (setup.py) diff --git a/dependencies-tests.log b/dependencies-tests.log index c5946011..2ed026dc 100644 --- a/dependencies-tests.log +++ b/dependencies-tests.log @@ -8,11 +8,11 @@ astropy==6.1.4 # via # healpy # icecube-skyreader (setup.py) -astropy-iers-data==0.2024.9.30.0.32.59 +astropy-iers-data==0.2024.10.7.0.32.46 # via astropy certifi==2024.8.30 # via requests -charset-normalizer==3.3.2 +charset-normalizer==3.4.0 # via requests contourpy==1.3.0 # via matplotlib @@ -32,7 +32,7 @@ matplotlib==3.9.2 # via icecube-skyreader (setup.py) meander==0.0.3 # via icecube-skyreader (setup.py) -numpy==2.1.1 +numpy==2.1.2 # via # astropy # contourpy @@ -83,5 +83,5 @@ tzdata==2024.2 # via pandas urllib3==2.2.3 # via requests -wipac-dev-tools==1.12.1 +wipac-dev-tools==1.13.0 # via icecube-skyreader (setup.py) diff --git a/dependencies.log b/dependencies.log index d1663d0a..fa83610f 100644 --- a/dependencies.log +++ b/dependencies.log @@ -8,11 +8,11 @@ astropy==6.1.4 # via # healpy # icecube-skyreader (setup.py) -astropy-iers-data==0.2024.9.30.0.32.59 +astropy-iers-data==0.2024.10.7.0.32.46 # via astropy certifi==2024.8.30 # via requests -charset-normalizer==3.3.2 +charset-normalizer==3.4.0 # via requests contourpy==1.3.0 # via matplotlib @@ -30,7 +30,7 @@ matplotlib==3.9.2 # via icecube-skyreader (setup.py) meander==0.0.3 # via icecube-skyreader (setup.py) -numpy==2.1.1 +numpy==2.1.2 # via # astropy # contourpy @@ -72,5 +72,5 @@ tzdata==2024.2 # via pandas urllib3==2.2.3 # via requests -wipac-dev-tools==1.12.1 +wipac-dev-tools==1.13.0 # via icecube-skyreader (setup.py) diff --git a/setup.cfg b/setup.cfg index e13907e9..60b0e44d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,7 @@ [wipac:cicd_setup_builder] pypi_name = icecube-skyreader python_min = 3.9 +python_max = 3.12 author = WIPAC Developers author_email = developers@icecube.wisc.edu keywords_spaced = "skymap scanner" skymap HEALPix neutrino reconstruction diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py index 0346ec33..8f00a0b3 100644 --- a/skyreader/plot/plot.py +++ b/skyreader/plot/plot.py @@ -2,14 +2,14 @@ # fmt: off # pylint: skip-file -# flake8: noqa import logging import pickle from pathlib import Path -from typing import Any, Dict, Final, List, Optional, Tuple, TypedDict, Union +from typing import List, Union import healpy # type: ignore[import] +# import mhealpy import matplotlib # type: ignore[import] import meander # type: ignore[import] import numpy as np @@ -17,21 +17,23 @@ from matplotlib import patheffects from matplotlib import pyplot as plt from matplotlib import text -from pathlib import Path +from matplotlib.projections import projection_registry # type: ignore[import] from .plotting_tools import ( AstroMollweideAxes, DecFormatter, - RaFormatter, format_fits_header, hp_ticklabels, - plot_catalog, + plot_catalog ) +from ..utils.areas import calculate_area, get_contour_areas +from ..utils.handle_map_data import extract_map, get_contour_levels from ..result import SkyScanResult LOGGER = logging.getLogger("skyreader.plot") + class SkyScanPlotter: PLOT_SIZE_Y_IN: float = 3.85 PLOT_SIZE_X_IN: float = 6 @@ -40,25 +42,19 @@ class SkyScanPlotter: PLOT_COLORMAP = matplotlib.colormaps['plasma_r'] def __init__(self, output_dir: Path = Path(".")): - # Set here plotting parameters and things that do not depend on the individual scan. + # Set here plotting parameters and things that + # do not depend on the individual scan. self.output_dir = output_dir - - @staticmethod - # Calculates are using Gauss-Green theorem / shoelace formula - # TODO: vectorize using numpy. - # Note: in some cases the argument is not a np.ndarray so one has to convert the data series beforehand. - def calculate_area(vs) -> float: - a = 0 - x0, y0 = vs[0] - for [x1,y1] in vs[1:]: - dx = np.cos(x1)-np.cos(x0) - dy = y1-y0 - a += 0.5*(y0*dx - np.cos(x0)*dy) - x0 = x1 - y0 = y1 - return a - - def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: + projection_registry.register(AstroMollweideAxes) + + def create_plot( + self, + result: SkyScanResult, + dozoom: bool = False, + systematics: bool = False, + llh_map: bool = True, + angular_error_floor: Union[None, float] = None, + ) -> 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 @@ -66,8 +62,9 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: """ 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 + # number of grid points along RA coordinate + xsize = int(self.PLOT_SIZE_X_IN * dpi) + ysize = int(xsize // 2) # number of grid points along dec coordinate dec = np.linspace(-np.pi/2., np.pi/2., ysize) ra = np.linspace(0., 2.*np.pi, xsize) # project the map to a rectangular matrix xsize x ysize @@ -82,101 +79,106 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: event_metadata = result.get_event_metadata() 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}" + run_str = f"Run: {event_metadata.run_id}" + evt_str = f"Event: {event_metadata.event_id}" + typ_str = f"Type: {event_metadata.event_type}" + mjd_str = f"MJD: {event_metadata.mjd}" + plot_title = f"{run_str} {evt_str} {typ_str} {mjd_str}" - plot_filename = f"{unique_id}.{'plot_zoomed_legacy.' if dozoom else ''}pdf" + if dozoom: + addition_to_filename = 'plot_zoomed_legacy.' + else: + addition_to_filename = '' + plot_filename = f"{unique_id}.{addition_to_filename}pdf" LOGGER.info(f"saving plot to {plot_filename}") - min_llh, max_llh = np.nan, np.nan - min_ra, min_dec = 0., 0. - - grid_map, grid_pix = None, None - - # now plot maps above each other - for nside in nsides: - LOGGER.info(f"constructing map for nside {nside}...") - # grid_pix = healpy.ang2pix(nside, THETA, PHI) - grid_pix = healpy.ang2pix(nside, np.pi/2. - DEC, RA) - this_map = np.ones(healpy.nside2npix(nside))*np.inf - - for pixel_data in result.result[f'nside-{nside}']: - pixel = pixel_data['index'] - # show 2*delta_LLH - value = 2*pixel_data['llh'] - if np.isfinite(value): - if np.isnan(min_llh) or value < min_llh: - minCoDec, min_ra = healpy.pix2ang(nside, pixel) - min_dec = np.pi/2 - minCoDec - min_llh = value - if np.isnan(max_llh) or value > max_llh: - max_llh = value - this_map[pixel] = value - - if grid_map is None: - grid_map = this_map[grid_pix] - else: - grid_map = np.where( np.isfinite(this_map[grid_pix]), this_map[grid_pix], grid_map) - - del this_map - - LOGGER.info(f"Completed map for nside {nside}.") + ( + grid_value, grid_ra, grid_dec, equatorial_map + ) = extract_map( + result, llh_map, angular_error_floor, remove_min_val=not llh_map + ) - # clean up - if grid_pix is not None: - del grid_pix + grid_pix = healpy.ang2pix(max(nsides), np.pi/2. - DEC, RA) + plotting_map = equatorial_map[grid_pix] - if grid_map is None: - # create an "empty" map if there are no pixels at all - grid_pix = healpy.ang2pix(8, np.pi/2 - DEC, RA) - this_map = np.ones(healpy.nside2npix(8))*np.inf - grid_map = this_map[grid_pix] - del this_map - del grid_pix + min_value = grid_value[0] # for probability map, this is actually + # the max_value + min_dec = grid_dec[0] + min_ra = grid_ra[0] - LOGGER.info(f"min RA: {min_ra *180./np.pi} deg, {min_ra*12./np.pi} hours.") + LOGGER.info( + f"min RA: {min_ra * 180./np.pi} deg, {min_ra*12./np.pi} hours." + ) 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 - - grid_map = np.ma.masked_invalid(grid_map) + plotting_map = plotting_map - min_value + equatorial_map = equatorial_map - min_value + vmin = 0. + vmax = 50 + map_to_plot = plotting_map + if llh_map: + cmap = self.PLOT_COLORMAP + text_colorbar = r"$-2 \ln(L)$" + vmin = np.nanmin(equatorial_map) + vmax = np.nanmax(equatorial_map) + map_to_plot = plotting_map + else: + cmap = matplotlib.colormaps[ + self.PLOT_COLORMAP.name.rstrip('_r') + ] + text_colorbar = r"log10$(p)$" + vmin = np.min(np.log10(equatorial_map)) + vmax = np.max(np.log10(equatorial_map)) + map_to_plot = np.log10(plotting_map) + equatorial_map = np.ma.masked_invalid(equatorial_map) + map_to_plot = np.ma.masked_invalid(map_to_plot) LOGGER.info(f"Preparing plot: {plot_filename}...") - # the color map to use - cmap = self.PLOT_COLORMAP - cmap.set_under(alpha=0.) # make underflows transparent - cmap.set_bad(alpha=1., color=(1.,0.,0.)) # make NaNs bright red + # features of the color map to use + 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 - fig = matplotlib.pyplot.figure(figsize=(self.PLOT_SIZE_X_IN,self.PLOT_SIZE_Y_IN)) + fig = matplotlib.pyplot.figure( + figsize=(self.PLOT_SIZE_X_IN, self.PLOT_SIZE_Y_IN) + ) ax = None if dozoom: - ax = fig.add_subplot(111) #,projection='cartesian') + ax = fig.add_subplot(111) # ,projection='cartesian') else: cmap.set_over(alpha=0.) # make underflows transparent - ax = fig.add_subplot(111,projection='astro mollweide') + ax = fig.add_subplot(111, projection='astro mollweide') # rasterized makes the map bitmap while the labels remain vectorial # flip longitude to the astro convention - image = ax.pcolormesh(ra, dec, grid_map, vmin=min_llh, vmax=max_llh, rasterized=True, cmap=cmap) + image = ax.pcolormesh( + ra, + dec, + map_to_plot, + vmin=vmin, + vmax=vmax, + rasterized=True, + cmap=cmap + ) # ax.set_xlim(np.pi, -np.pi) + ( + contour_levels, contour_labels, contour_colors + ) = get_contour_levels(equatorial_map, llh_map, systematics) - - contour_levels = (np.array([1.39, 4.61, 11.83, 28.74])+min_llh)[:2] - contour_labels = [r'50%', r'90%', r'3$\sigma$', r'5$\sigma$'][:2] - contour_colors=['k', 'r', 'g', 'b'][:2] - leg_element=[] + leg_element = [] cs_collections = [] for level, color in zip(contour_levels, contour_colors): - contour_set = ax.contour(ra, dec, grid_map, levels=[level], colors=[color]) + if not llh_map: + level = np.log10(level) + contour_set = ax.contour( + ra, dec, map_to_plot, levels=[level], colors=[color] + ) cs_collections.append(contour_set.collections[0]) e, _ = contour_set.legend_elements() leg_element.append(e[0]) @@ -187,53 +189,85 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: # 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)$") + cb = fig.colorbar( + image, + orientation='horizontal', + shrink=.6, + pad=0.05, + ticks=[vmin, vmax], + ) + cb.ax.xaxis.set_label_text(text_colorbar) else: ax.set_xlabel('right ascension') ax.set_ylabel('declination') - cb = fig.colorbar(image, orientation='horizontal', shrink=.6, pad=0.13) + cb = fig.colorbar( + image, orientation='horizontal', shrink=.6, pad=0.13 + ) cb.ax.xaxis.set_label_text(r"$-2 \Delta \ln (L)$") leg_labels = [] for i in range(len(contour_labels)): vs = cs_collections[i].get_paths()[0].vertices # Compute area enclosed by vertices. - # Take absolute values to be independent of orientation of the boundary integral. - contour_area = abs(self.calculate_area(vs)) # will be in square-radians - contour_area_sqdeg = contour_area*(180.*180.)/(np.pi*np.pi) # convert to square-degrees - - leg_labels.append(f'{contour_labels[i]} - area: {contour_area_sqdeg:.2f}sqdeg') + # Take absolute values to be independent of orientation of + # the boundary integral. + contour_area = abs(calculate_area(vs)) # in square-radians + # convert to square-degrees + contour_area_sqdeg = contour_area*(180.*180.)/(np.pi*np.pi) + + area_string = f"area: {contour_area_sqdeg:.2f}sqdeg" + leg_labels.append( + f'{contour_labels[i]} - {area_string}' + ) - ax.scatter(min_ra, min_dec, s=20, marker='*', color='black', label=r'scan best-fit', zorder=2) - ax.legend(leg_element, leg_labels, loc='lower right', fontsize=8, scatterpoints=1, ncol=2) + ax.scatter( + min_ra, + min_dec, + s=20, + marker='*', + color='black', + label=r'scan best-fit', + zorder=2 + ) + ax.legend( + leg_element, + leg_labels, + loc='lower right', + fontsize=8, + scatterpoints=1, + ncol=2 + ) - LOGGER.info(f"Contour Area (90%): {contour_area_sqdeg} degrees (cartesian) {contour_area_sqdeg * np.cos(min_dec)**2} degrees (scaled)") - + LOGGER.info(f"Contour Area (90%): {contour_area_sqdeg} " + f"degrees (cartesian) " + f"{contour_area_sqdeg * np.cos(min_dec)**2} " + "degrees (scaled)") x_width = 1.6 * np.sqrt(contour_area_sqdeg) LOGGER.info(f"x width is {x_width}") - if np.isnan(x_width): - # this get called only when contour_area / x_width is NaN so possibly never invoked in typical situations - - raise RuntimeError("Estimated area / width is NaN and the fallback logic for this scenario is no longer supported. If you encounter this error raise an issue to SkyReader.") - - # mypy error: "QuadContourSet" has no attribute "allsegs" [attr-defined] - # this attribute is likely deprecated but this scenario is rarely (if ever) hit - # original code is kept commented for the time being - - # note: contour_set is re-assigned at every iteration of the loop on - # contour_levels, contour_colors, so this effectively corresponds to the - # the last contour_set - - # x_width = 1.6*(max(contour_set.allsegs[i][0][:,0]) - min(contour_set.allsegs[i][0][:,0])) + # this get called only when contour_area / x_width is + # NaN so possibly never invoked in typical situations + raise RuntimeError( + "Estimated area / width is NaN and the fallback logic " + "for this scenario is no longer supported. If you " + "encounter this error raise an issue to SkyReader." + ) + # mypy error: "QuadContourSet" has no attribute "allsegs" + # [attr-defined]. This attribute is likely deprecated but + # this scenario is rarely (if ever) hit original code is + # kept commented for the time being + # note: contour_set is re-assigned at every iteration of + # the loop on contour_levels, contour_colors, so this + # effectively corresponds to the last contour_set + # x_width = 1.6*(max(contour_set.allsegs[i][0][:,0]) - + # min(contour_set.allsegs[i][0][:,0])) y_width = 0.5 * x_width - lower_x = max(min_ra - x_width*np.pi/180., 0.) - upper_x = min(min_ra + x_width*np.pi/180., 2 * np.pi) - lower_y = max(min_dec -y_width*np.pi/180., -np.pi/2.) + lower_x = max(min_ra - x_width*np.pi/180., 0.) + upper_x = min(min_ra + x_width*np.pi/180., 2 * np.pi) + lower_y = max(min_dec - y_width*np.pi/180., -np.pi/2.) upper_y = min(min_dec + y_width*np.pi/180., np.pi/2.) ax.set_xlim(upper_x, lower_x) @@ -246,22 +280,26 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: factor = 0.25*(np.pi/180.) while (upper_x - lower_x)/factor > 6: - factor *= 2. + factor *= 2. tick_label_grid = factor - ax.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(base=tick_label_grid)) - ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(base=tick_label_grid)) + 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 - # mypy compliance: since cb.solids could be None, we check that it is actually - # a valid object before accessing it + # mypy compliance: since cb.solids could be None, we check that + # it is actually a valid object before accessing it if isinstance(cb.solids, matplotlib.collections.QuadMesh): cb.solids.set_edgecolor("face") if dozoom: ax.set_aspect('equal') - + ax.tick_params(axis='x', labelsize=10) ax.tick_params(axis='y', labelsize=10) @@ -269,18 +307,32 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: ax.grid(True, color='k', alpha=0.5) # Otherwise, add the path effects. - # mypy requires set_path_effects() to take a list of AbstractPathEffect - effects: List[patheffects.AbstractPathEffect] = [patheffects.withStroke(linewidth=1.1, foreground='w')] + # mypy requires set_path_effects() to take a list of AbstractPathEffect + effects: List[patheffects.AbstractPathEffect] = [ + patheffects.withStroke(linewidth=1.1, foreground='w') + ] for artist in ax.findobj(text.Text): - # mypy error: Argument 1 to "set_path_effects" of "Artist" has incompatible type "list[withStroke]"; expected "list[AbstractPathEffect]" [arg-type] + # mypy error: Argument 1 to "set_path_effects" of "Artist" + # has incompatible type "list[withStroke]"; expected + # "list[AbstractPathEffect]" [arg-type] artist.set_path_effects(effects) # 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) + 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) + fig.subplots_adjust( + bottom=spacing, + top=0.92-spacing, + left=spacing+0.1, + right=1.-spacing + ) # set the title fig.suptitle(plot_title) @@ -293,231 +345,273 @@ def create_plot(self, result: SkyScanResult, dozoom: bool = False) -> None: @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, - 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): + """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 + comp1 = np.cos(phi)*np.sin(ra)*np.sin(delta) + comp2 = np.cos(ra)*np.cos(delta)*np.sin(dec) + comp3 = np.cos(dec)*np.sin(delta)*np.sin(phi) + comp4 = np.cos(delta)*np.sin(dec)*np.sin(ra) + comp5 = np.sin(delta)*np.cos(ra)*np.cos(phi) + comp6 = np.sin(delta)*np.cos(dec)*np.sin(ra)*np.sin(phi) + comp7 = np.cos(dec)*np.cos(delta) + comp8 = np.sin(dec)*np.sin(delta)*np.sin(phi) + vx = comp1 + comp2 + comp3 + vy = comp4 - comp5 + comp6 + vz = comp7 - comp8 + 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, + 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, + angular_error_floor=None, # if not None, sigma of the + # gaussian to convolute the map with in deg. + llh_map=True + ): """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 + 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 dec_minus = (np.pi/2-np.max(theta))*180./np.pi - dec return ra_plus, ra_minus, dec_plus, dec_minus dpi = self.PLOT_DPI_ZOOMED - lonra=[-10.,10.] - latra=[-10.,10.] + lonra = [-10., 10.] + latra = [-10., 10.] event_metadata = result.get_event_metadata() 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}" + run_str = f"Run: {event_metadata.run_id}" + evt_str = f"Event: {event_metadata.event_id}" + typ_str = f"Type: {event_metadata.event_type}" + mjd_str = f"MJD: {event_metadata.mjd}" + plot_title = f"{run_str} {evt_str} {typ_str} {mjd_str}" nsides = result.nsides + max_nside = max(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" - LOGGER.info(f"saving plot to {plot_filename}") - nsides = result.nsides - LOGGER.info(f"available nsides: {nsides}") - - grid_map = dict() - max_nside = max(nsides) - equatorial_map = np.full(healpy.nside2npix(max_nside), np.nan) - - for nside in nsides: - LOGGER.info(f"constructing map for nside {nside}...") - npix = healpy.nside2npix(nside) - - map_data = result.result[f'nside-{nside}'] - pixels = map_data['index'] - values = map_data['llh'] - this_map = np.full(npix, np.nan) - this_map[pixels] = values - if nside < max_nside: - this_map = healpy.ud_grade(this_map, max_nside) - mask = np.logical_and(~np.isnan(this_map), np.isfinite(this_map)) - equatorial_map[mask] = this_map[mask] - - for pixel_data in result.result[f"nside-{nside}"]: - pixel = pixel_data['index'] - value = pixel_data['llh'] - if np.isfinite(value) and not np.isnan(value): - tmp_theta, tmp_phi = healpy.pix2ang(nside, pixel) - tmp_dec = np.pi/2 - tmp_theta - tmp_ra = tmp_phi - grid_map[(tmp_dec, tmp_ra)] = value - LOGGER.info(f"done with map for nside {nside}...") - - grid_dec_list, grid_ra_list, grid_value_list = [], [], [] - - for (dec, ra), value in grid_map.items(): - grid_dec_list.append(dec); grid_ra_list.append(ra) - grid_value_list.append(value) - grid_dec: np.ndarray = np.asarray(grid_dec_list) - grid_ra: np.ndarray = np.asarray(grid_ra_list) - grid_value: np.ndarray = np.asarray(grid_value_list) - - sorting_indices = np.argsort(grid_value) - grid_value = grid_value[sorting_indices] - grid_dec = grid_dec[sorting_indices] - grid_ra = grid_ra[sorting_indices] - - min_value = grid_value[0] + ( + grid_value, grid_ra, grid_dec, equatorial_map + ) = extract_map(result, llh_map, angular_error_floor) min_dec = grid_dec[0] min_ra = grid_ra[0] - LOGGER.info(f"min RA: {min_ra *180./np.pi} deg, {min_ra*12./np.pi} hours.") + LOGGER.info( + f"min RA: {min_ra * 180./np.pi} deg, {min_ra*12./np.pi} hours." + ) LOGGER.info(f"min Dec: {min_dec * 180./np.pi} deg") - # renormalize - grid_value = grid_value - min_value - min_value = 0. + ( + contour_levels, contour_labels, contour_colors + ) = get_contour_levels(equatorial_map, llh_map, systematics) - # show 2 * delta_LLH - grid_value = grid_value * 2. + sample_points = np.array([np.pi/2 - grid_dec, grid_ra]).T - # Do same for the healpy map - equatorial_map[np.isinf(equatorial_map)] = np.nan - equatorial_map -= np.nanmin(equatorial_map) - equatorial_map *= 2. + if not circular: + if llh_map: + grid_values_for_contours = grid_value + contour_levels_for_contours = contour_levels + else: + grid_values_for_contours = np.log(grid_value) + contour_levels_for_contours = np.log(contour_levels) + # Get contours from healpix map + contours_by_level = meander.spherical_contours( + sample_points, + grid_values_for_contours, + contour_levels_for_contours, + ) + LOGGER.info(f"saving plot to {plot_filename}") LOGGER.info(f"preparing plot: {plot_filename}...") - cmap = self.PLOT_COLORMAP - cmap.set_under('w') - cmap.set_bad(alpha=1., color=(1.,0.,0.)) # make NaNs bright red - - # Calculate the contours - if systematics: - # from Pan-Starrs event 127852 - contour_levels = (np.array([22.2, 64.2])+min_value) # these are values determined from MC by Will on the TS (2*LLH) - contour_labels = [r'50% (IC160427A syst.)', r'90% (IC160427A syst.)'] - contour_colors=['k', 'r'] - else: - # Wilks - contour_levels = (np.array([1.39, 4.61, 11.83, 28.74])+min_value)[:3] - contour_labels = [r'50%', r'90%', r'3$\sigma$', r'5$\sigma$'][:3] - contour_colors=['k', 'r', 'g', 'b'][:3] - - sample_points = np.array([np.pi/2 - grid_dec, grid_ra]).T - - # Call meander module to find contours - if not circular: - contours_by_level = meander.spherical_contours(sample_points, - grid_value, contour_levels - ) + # In case it is desired, just draw circular contours over the ts map 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) - contour50 = np.dstack((Theta50,Phi50)) - contour90 = np.dstack((Theta90,Phi90)) + Theta50, Phi50 = self.circular_contour( + min_ra, min_dec, sigma50, max_nside + ) + Theta90, Phi90 = self.circular_contour( + min_ra, min_dec, sigma90, max_nside + ) + contour50 = np.dstack((Theta50, Phi50)) + contour90 = np.dstack((Theta90, Phi90)) contours_by_level = [contour50, contour90] - + + # Calculate areas using Gauss-Green's theorem for a spherical space + contour_areas = get_contour_areas(contours_by_level, min_ra) + # Check for RA values that are out of bounds for level in contours_by_level: for contour in level: - contour.T[1] = np.where(contour.T[1] < 0., - contour.T[1] + 2.*np.pi, contour.T[1] - ) - + contour.T[1] = np.where( + contour.T[1] < 0., + contour.T[1] + 2.*np.pi, + contour.T[1] + ) # Find the rough extent of the contours to bound the plot contours = contours_by_level[-1] ra = min_ra * 180./np.pi dec = min_dec * 180./np.pi theta, phi = np.concatenate(contours_by_level[-1]).T - ra_plus, ra_minus, dec_plus, dec_minus = bounding_box(ra, dec, theta, phi) + ra_plus, ra_minus, dec_plus, dec_minus = bounding_box( + ra, dec, theta, phi + ) ra_bound = min(15, max(3, max(-ra_minus, ra_plus))) dec_bound = min(15, max(2, max(-dec_minus, dec_plus))) lonra = [-ra_bound, ra_bound] latra = [-dec_bound, dec_bound] - #Begin the figure + # Begin the figure plt.clf() # Rotate into healpy coordinates lon, lat = np.degrees(min_ra), np.degrees(min_dec) - healpy.cartview(map=equatorial_map, title=plot_title, - min=0., #min 2DeltaLLH value for colorscale - max=40., #max 2DeltaLLH value for colorscale - rot=(lon,lat,0.), cmap=cmap, hold=True, - cbar=None, lonra=lonra, latra=latra, - unit=r"$-2 \Delta \ln (L)$", + if llh_map: + cmap = self.PLOT_COLORMAP + cmap.set_under('w') + # make NaNs bright red + cmap.set_bad(alpha=1., color=(1., 0., 0.)) + healpy.cartview( + map=equatorial_map, + title=plot_title, + min=0., # min 2DeltaLLH value for colorscale + max=40., # max 2DeltaLLH value for colorscale + rot=(lon, lat, 0.), + cmap=cmap, + hold=True, + cbar=None, + lonra=lonra, + latra=latra, + unit=r"$-2 \Delta \ln (L)$", ) + ticks = None + format = None + cb_label = r"$-2 \Delta \ln (L)$" + else: + cmap = matplotlib.colormaps[ + self.PLOT_COLORMAP.name.rstrip('_r') + ] + cmap.set_under('w') + # make NaNs bright red + cmap.set_bad(alpha=1., color=(1., 0., 0.)) + max_prob = np.nanmax(equatorial_map) + if len(contour_levels) >= 3: + min_prob = contour_levels[2] + else: + min_prob = max_prob/1e8 + healpy.cartview( + map=equatorial_map.clip(1e-12, None), + title=plot_title, + min=min_prob, # min prob value for colorscale + max=max_prob, # max prob value for colorscale + rot=(lon, lat, 0.), + cmap=cmap, + hold=True, + cbar=None, + lonra=lonra, + latra=latra, + norm='log', + unit="Probability", + ) + ticks = [ + min_prob, + min_prob*(max_prob/min_prob)**(1/5), + min_prob*(max_prob/min_prob)**(2/5), + min_prob*(max_prob/min_prob)**(3/5), + min_prob*(max_prob/min_prob)**(4/5), + max_prob + ] + format = "{x:.1e}" + cb_label = "Probability" fig = plt.gcf() ax = plt.gca() image = ax.get_images()[0] # Place colorbar by hand - cb = fig.colorbar(image, ax=ax, orientation='horizontal', aspect=50) - cb.ax.xaxis.set_label_text(r"$-2 \Delta \ln (L)$") + cb = fig.colorbar( + image, + ax=ax, + orientation='horizontal', + aspect=50, + ticks=ticks, + format=format + ) + cb.ax.xaxis.set_label_text(cb_label) # Plot the best-fit location # This requires some more coordinate transformations - healpy.projplot(np.pi/2 - min_dec, min_ra, - '*', ms=5, label=r'scan best fit', color='black', zorder=2) + healpy.projplot( + np.pi/2 - min_dec, + min_ra, + '*', + ms=5, + label=r'scan best fit', + color='black', + zorder=2 + ) # Plot the contours - contour_areas=[] - for contour_label, contour_color, contours in zip( - contour_labels, contour_colors, contours_by_level): - contour_area = 0. - for contour in contours: - _ = contour.copy() - _[:,1] += np.pi-np.radians(ra) - _[:,1] %= 2*np.pi - contour_area += self.calculate_area(_) - contour_area_sqdeg = abs(contour_area) * (180.*180.)/(np.pi*np.pi) # convert to square-degrees - contour_areas.append(contour_area_sqdeg) + for ( + contour_area_sqdeg, + contour_label, + contour_color, + contours + ) in zip( + contour_areas, + contour_labels, + contour_colors, + contours_by_level + ): contour_label = contour_label + ' - area: {0:.2f} sqdeg'.format( contour_area_sqdeg) first = True for contour in contours: theta, phi = contour.T if first: - healpy.projplot(theta, phi, linewidth=2, c=contour_color, - label=contour_label) + healpy.projplot( + theta, + phi, + linewidth=2, + c=contour_color, + label=contour_label + ) else: healpy.projplot(theta, phi, linewidth=2, c=contour_color) first = False @@ -539,48 +633,82 @@ def bounding_box(ra, dec, theta, phi): upper_lat = max(tmp_lower_lat, tmp_upper_lat) # Label the axes - hp_ticklabels(zoom=True, lonra=lonra, latra=latra, - rot=(lon,lat,0), - bounds=(lower_lon, upper_lon, lower_lat, upper_lat)) + hp_ticklabels( + zoom=True, + lonra=lonra, + latra=latra, + rot=(lon, lat, 0), + bounds=(lower_lon, upper_lon, lower_lat, upper_lat) + ) if plot_4fgl: # Overlay 4FGL sources - plot_catalog(equatorial_map, cmap, lower_ra, upper_ra, lower_dec, upper_dec) + plot_catalog( + equatorial_map, cmap, lower_ra, upper_ra, lower_dec, upper_dec + ) # Approximate contours as rectangles ra = min_ra * 180./np.pi dec = min_dec * 180./np.pi - for l, contours in enumerate(contours_by_level[:2]): + percentages = ["50", "90"] + for l_index, contours in enumerate(contours_by_level[:2]): ra_plus = None theta, phi = np.concatenate(contours).T - ra_plus, ra_minus, dec_plus, dec_minus = bounding_box(ra, dec, theta, phi) - contain_txt = "Approximating the {0}% error region as a rectangle, we get:".format(["50", "90"][l]) + " \n" + \ + ra_plus, ra_minus, dec_plus, dec_minus = bounding_box( + ra, dec, theta, phi + ) + 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)) # This is actually an output and not a logging info. - # TODO: we should wrap this in an object, return and log at the higher level. + # TODO: we should wrap this in an object, return and log at + # the higher level. print(contain_txt) + print( + "Contour Area (50%):", + contour_areas[0], + "square degrees (scaled)" + ) + print( + "Contour Area (90%):", + contour_areas[1], + "square degrees (scaled)" + ) + if plot_bounding_box: bounding_ras_list, bounding_decs_list = [], [] # lower bound - bounding_ras_list.extend(list(np.linspace(ra+ra_minus, - ra+ra_plus, 10))) + bounding_ras_list.extend(list(np.linspace( + ra+ra_minus, + ra+ra_plus, + 10 + ))) bounding_decs_list.extend([dec+dec_minus]*10) # right bound bounding_ras_list.extend([ra+ra_plus]*10) - bounding_decs_list.extend(list(np.linspace(dec+dec_minus, - dec+dec_plus, 10))) + bounding_decs_list.extend(list(np.linspace( + dec+dec_minus, + dec+dec_plus, + 10 + ))) # upper bound - bounding_ras_list.extend(list(np.linspace(ra+ra_plus, - ra+ra_minus, 10))) + bounding_ras_list.extend(list(np.linspace( + ra+ra_plus, + ra+ra_minus, + 10 + ))) bounding_decs_list.extend([dec+dec_plus]*10) # left bound bounding_ras_list.extend([ra+ra_minus]*10) - bounding_decs_list.extend(list(np.linspace(dec+dec_plus, - dec+dec_minus,10))) + bounding_decs_list.extend(list(np.linspace( + dec+dec_plus, + dec+dec_minus, + 10 + ))) # join end to beginning bounding_ras_list.append(bounding_ras_list[0]) bounding_decs_list.append(bounding_decs_list[0]) @@ -591,12 +719,19 @@ def bounding_box(ra, dec, theta, phi): bounding_theta = np.pi/2 - np.radians(bounding_decs) bounding_contour = np.array([bounding_theta, bounding_phi]) bounding_contour_area = 0. - bounding_contour_area = abs(self.calculate_area(bounding_contour.T)) - bounding_contour_area *= (180.*180.)/(np.pi*np.pi) # convert to square-degrees - contour_label = r'90% Bounding rectangle' + ' - area: {0:.2f} sqdeg'.format( - bounding_contour_area) - healpy.projplot(bounding_theta, bounding_phi, linewidth=0.75, - c='r', linestyle='dashed', label=contour_label) + 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' + \ + ' - area: {0:.2f} sqdeg'.format(bounding_contour_area) + healpy.projplot( + bounding_theta, + bounding_phi, + linewidth=0.75, + c='r', + linestyle='dashed', + label=contour_label + ) # Output contours in RA, dec instead of theta, phi saving_contours: list = [] @@ -617,27 +752,35 @@ 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 {self.output_dir / savename}") - ascii.write(tab, self.output_dir / savename, overwrite=True) - except OSError as err: - LOGGER.error("OS Error prevented contours from being written, maybe a memory issue. Error is:\n{err}") + LOGGER.info("Dumping to {savename}") + ascii.write(tab, savename, overwrite=True) + except OSError: + LOGGER.error( + "OS Error prevented contours from being written, " + "maybe a memory issue. Error is:\n{err}") uncertainty = [(ra_minus, ra_plus), (dec_minus, dec_plus)] fits_header = format_fits_header( - (event_metadata.run_id, event_metadata.event_id, event_metadata.event_type), - 0, + ( + event_metadata.run_id, + event_metadata.event_id, + event_metadata.event_type + ), + event_metadata.mjd, np.degrees(min_ra), np.degrees(min_dec), uncertainty, + llh_map, ) mmap_nside = healpy.get_nside(equatorial_map) # Plot the original online reconstruction location 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") - + # 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") extra_ra_rad = np.radians(extra_ra) extra_dec_rad = np.radians(extra_dec) @@ -645,39 +788,90 @@ def bounding_box(ra, dec, theta, phi): 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)) + 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_scale, cont_col, cont_sty in zip( + [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] 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) + healpy.projplot( + np.degrees(spline_lon), + np.degrees(spline_lat), + lonlat=True, + linewidth=2., + color=cont_col, + linestyle=cont_sty + ) plt.legend(fontsize=6, loc="lower left") - - print("Contour Area (50%):", contour_areas[0], "degrees (scaled)") - print("Contour Area (90%):", contour_areas[1], "degrees (scaled)") - - # Dump the whole contour - pickle_path = self.output_dir / (unique_id + ".contour.pkl") - LOGGER.info(f"Saving contour to {pickle_path}") - with open(pickle_path, "wb") as f: + path = unique_id + ".contour.pkl" + print("Saving contour to", path) + with open(path, "wb") as f: pickle.dump(saving_contours, f) - healpy.write_map(self.output_dir / f"{unique_id}.skymap_nside_{mmap_nside}.fits.gz", - equatorial_map, coord = 'C', column_names = ['2LLH'], - extra_header = fits_header, overwrite=True) + # logic for multiordermaps -> for future developements + + # save multiorder version of the map + # multiorder_map = mhealpy.HealpixMap( + # grid_value, uniq_array + # ) + # multiorder_map.write_map( + # f"{unique_id}.skymap_nside_{mmap_nside}.multiorder.fits.gz", + # column_names=["PROBABILITY"], + # extra_header=fits_header, + # overwrite=True, + # ) + + if llh_map: + column_names = ['2DLLH'] + else: + # avoid excessively heavy data format for the flattened map + equatorial_map[ + equatorial_map < 1e-16 + ] = np.mean( + equatorial_map[equatorial_map < 1e-16] + ) + column_names = ["PROBABILITY"] + + # save flattened map + if llh_map: + type_map = "llh" + else: + type_map = "probability" + healpy.write_map( + f"{unique_id}.skymap_nside_{mmap_nside}_{type_map}.fits.gz", + equatorial_map, + coord='C', + column_names=column_names, + extra_header=fits_header, + overwrite=True + ) # Save the figure LOGGER.info(f"Saving: {plot_filename}...") - #ax.invert_xaxis() - fig.savefig(self.output_dir / plot_filename, dpi=dpi, transparent=True, bbox_inches='tight') + # ax.invert_xaxis() + fig.savefig( + self.output_dir / plot_filename, + dpi=dpi, + transparent=True, + bbox_inches='tight' + ) LOGGER.info("done.") diff --git a/skyreader/plot/plotting_tools.py b/skyreader/plot/plotting_tools.py index 4c49c30c..29a6ef33 100644 --- a/skyreader/plot/plotting_tools.py +++ b/skyreader/plot/plotting_tools.py @@ -10,18 +10,22 @@ import numpy as np from matplotlib import pyplot as plt from matplotlib.axes import Axes # type: ignore[import] -from matplotlib.projections import projection_registry # type: ignore[import] from matplotlib.projections.geo import MollweideAxes # type: ignore[import] from matplotlib.ticker import FixedLocator, Formatter # type: ignore[import] from matplotlib.transforms import Affine2D # type: ignore[import] matplotlib.use('agg') -def format_fits_header(event_id_tuple, mjd, ra, dec, uncertainty): +def format_fits_header(event_id_tuple, mjd, ra, dec, uncertainty, llh_map): """Prepare some of the relevant event information for a fits file header.""" run_id, event_id, event_type = event_id_tuple + if llh_map: + uncertainty_comment = 'Change in 2LLH based on Wilks theorem' + else: + uncertainty_comment = 'Probability-per-pixel (all pixels with same area)' + header = [ ('RUNID', run_id), ('EVENTID', event_id), @@ -39,7 +43,7 @@ def format_fits_header(event_id_tuple, mjd, ra, dec, uncertainty): ('DEC_ERR_MINUS', np.round(np.abs(uncertainty[1][0]),2), '90% containment error low'), ('COMMENTS', '50%(90%) uncertainty location' \ - + ' => Change in 2LLH based on Wilks theorem'), + + ' => ' + uncertainty_comment), ('NOTE', 'Please ignore pixels with infinite or NaN values.' \ + ' They are rare cases of the minimizer failing to converge') ] @@ -111,7 +115,7 @@ def plot_catalog(master_map, cmap, lower_ra, upper_ra, lower_dec, upper_dec, cmap_min=0., cmap_max=250.): """"Plots the 4FGL catalog in a color that contrasts with the background healpix map.""" - hdu = pyfits.open('/cvmfs/icecube.opensciencegrid.org/users/azegarelli/realtime/catalogs/gll_psc_v34.fit') ## LAT 14-year Source Catalog (4FGL-DR4 in FITS format) ; https://fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/ + hdu = pyfits.open('/cvmfs/icecube.opensciencegrid.org/users/azegarelli/realtime/catalogs/gll_psc_v34.fit') ## LAT 14-year Source Catalog (4FGL-DR4 in FITS format) ; https://fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/') fgl = hdu[1] pe = [path_effects.Stroke(linewidth=0.5, foreground=cmap(0.0)), path_effects.Normal()] @@ -150,7 +154,6 @@ def color_filter(lon, lat): path_effects=pe) del fgl - ## # Mollweide axes with phi axis flipped and in hours from 24 to 0 instead of # in degrees from -180 to 180. @@ -215,7 +218,7 @@ def set_longitude_grid(self, degrees): FixedLocator( np.linspace(0, 2*np.pi, number, True)[1:-1])) self._longitude_degrees = degrees - self.xaxis.set_major_formatter(self.RaFormatter(degrees)) + self.xaxis.set_major_formatter(self.RaFormatter(degrees)) def _set_lim_and_transforms(self): # Copied from matplotlib.geo.GeoAxes._set_lim_and_transforms and modified @@ -250,6 +253,4 @@ def _get_affine_transform(self): _, yscale = transform.transform_point((0, np.pi / 2.0)) return Affine2D() \ .scale(0.5 / xscale, 0.5 / yscale) \ - .translate(0.5, 0.5) - -projection_registry.register(AstroMollweideAxes) + .translate(0.5, 0.5) \ No newline at end of file diff --git a/skyreader/utils/__init__.py b/skyreader/utils/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/skyreader/utils/areas.py b/skyreader/utils/areas.py new file mode 100644 index 00000000..9b5b3135 --- /dev/null +++ b/skyreader/utils/areas.py @@ -0,0 +1,35 @@ +import numpy as np +from typing import List + + +# Calculates are using Gauss-Green theorem / shoelace formula +# TODO: vectorize using numpy. +# Note: in some cases the argument is not a np.ndarray so one has +# to convert the data series beforehand. +def calculate_area(vs) -> float: + a = 0 + x0, y0 = vs[0] + for [x1, y1] in vs[1:]: + dx = np.cos(x1)-np.cos(x0) + dy = y1-y0 + a += 0.5*(y0*dx - np.cos(x0)*dy) + x0 = x1 + y0 = y1 + return a + + +# Given a list of contours by level, returns the areas of the contours +def get_contour_areas(contours_by_level_list, min_ra) -> List[float]: + contour_areas = [] + ra = min_ra * 180./np.pi + for contours in contours_by_level_list: + contour_area = 0. + for contour in contours: + _ = contour.copy() + _[:,1] += np.pi-np.radians(ra) + _[:,1] %= 2*np.pi + contour_area += abs(calculate_area(_)) + # convert to square-degrees + contour_area_sqdeg = abs(contour_area) * (180.*180.)/(np.pi*np.pi) + contour_areas.append(contour_area_sqdeg) + return contour_areas diff --git a/skyreader/utils/handle_map_data.py b/skyreader/utils/handle_map_data.py new file mode 100644 index 00000000..929012c9 --- /dev/null +++ b/skyreader/utils/handle_map_data.py @@ -0,0 +1,202 @@ +from typing import Union +import logging + +import healpy # type: ignore[import] +import numpy as np + +from ..result import SkyScanResult + +LOGGER = logging.getLogger("skyreader.extract_map") + + +def extract_map( + result: SkyScanResult, + llh_map: bool = True, + angular_error_floor: Union[None, float] = None, + remove_min_val: bool = True, +): + """ + Extract from the output of skymap_scanner the healpy map + args: + - result: SkyScanResult. The output of the Skymap Scanner + - llh_map: bool = True. If True the likelihood will be plotted, + otherwise the probability + - angular_error_floor: Union[None, float] = None. if not None, + sigma of the gaussian to convolute the map with in deg. + - remove_min_val: bool = True. Remove minimum value from -llh + no effect if probability map. + + returns: + - grid_value: value-per-scanned-pixel (pixels with + different nsides) + - grid_ra: right ascension for each pixel in grid_value + - grid_dec: declination for each pixel in grid_value + - equatorial_map: healpix map with maximum nside (all pixels + with same nside) + """ + + grid_map = dict() + + nsides = result.nsides + max_nside = max(nsides) + equatorial_map = np.full(healpy.nside2npix(max_nside), np.nan) + + for nside in nsides: + LOGGER.info(f"constructing map for nside {nside}...") + npix = healpy.nside2npix(nside) + map_data = result.result[f'nside-{nside}'] + pixels = map_data['index'] + values = map_data['llh'] + this_map = np.full(npix, np.nan) + this_map[pixels] = values + if nside < max_nside: + this_map = healpy.ud_grade(this_map, max_nside) + mask = np.logical_and(~np.isnan(this_map), np.isfinite(this_map)) + equatorial_map[mask] = this_map[mask] + + for pixel_data in result.result[f"nside-{nside}"]: + pixel = pixel_data['index'] + value = pixel_data['llh'] + if np.isfinite(value) and not np.isnan(value): + tmp_theta, tmp_phi = healpy.pix2ang(nside, pixel) + tmp_dec = np.pi/2 - tmp_theta + tmp_ra = tmp_phi + grid_map[(tmp_dec, tmp_ra)] = value + # nested_pixel = healpy.ang2pix( + # nside, tmp_theta, tmp_phi, nest=True + # ) + # uniq = 4*nside*nside + nested_pixel + # uniq_list.append(uniq) + LOGGER.info(f"done with map for nside {nside}...") + + grid_dec_list, grid_ra_list, grid_value_list = [], [], [] + + for (dec, ra), value in grid_map.items(): + grid_dec_list.append(dec) + grid_ra_list.append(ra) + grid_value_list.append(value) + grid_dec: np.ndarray = np.asarray(grid_dec_list) + grid_ra: np.ndarray = np.asarray(grid_ra_list) + grid_value: np.ndarray = np.asarray(grid_value_list) + # uniq_array: np.ndarray = np.asarray(uniq_list) + + sorting_indices = np.argsort(grid_value) + grid_value = grid_value[sorting_indices] + grid_dec = grid_dec[sorting_indices] + grid_ra = grid_ra[sorting_indices] + # uniq_array = uniq_array[sorting_indices] + + min_value = grid_value[0] + + if remove_min_val or (not llh_map): + # renormalize + grid_value = grid_value - min_value + min_value = 0. + + # renormalize + equatorial_map[np.isinf(equatorial_map)] = np.nan + equatorial_map -= np.nanmin(equatorial_map) + + if llh_map: + # show 2 * delta_LLH + grid_value = grid_value * 2. + equatorial_map *= 2. + else: + # Convert to probability + equatorial_map = np.exp(-1. * equatorial_map) + equatorial_map = equatorial_map / np.nansum(equatorial_map) + + # nan values are a problem for the convolution and the contours + min_map = np.nanmin(equatorial_map) + equatorial_map[np.isnan(equatorial_map)] = min_map + + if angular_error_floor is not None: + # convolute with a gaussian. angular_error_floor is the + # sigma in deg. + equatorial_map = healpy.smoothing( + equatorial_map, + sigma=np.deg2rad(angular_error_floor), + ) + + # normalize map + min_map = np.nanmin(equatorial_map[equatorial_map > 0.0]) + equatorial_map[np.isnan(equatorial_map)] = min_map + equatorial_map = equatorial_map.clip(min_map, None) + normalization = np.nansum(equatorial_map) + equatorial_map = equatorial_map / normalization + + # obtain values for grid map + grid_value = healpy.get_interp_val( + equatorial_map, np.pi/2 - grid_dec, grid_ra + ) + grid_value[np.isnan(grid_value)] = min_map + grid_value = grid_value.clip(min_map, None) + + return grid_value, grid_ra, grid_dec, equatorial_map + + +def get_contour_levels( + equatorial_map: np.ndarray, + llh_map: bool = True, + systematics: bool = False, +): + """ + get contour levels for the desired map + + args: + - equatorial_map: np.ndarray. Necessary in case of + probability map. + - llh_map: bool. If True llh levels, otherwise probability + levels. + - systematics: bool. Only for llh maps. If True include + recalibrated llh values from Pan-Starrs event 127852 + (IC160427A syst.) + + returns: + - contour_levels, levels of the contours + - contour_labels, respective labels for the contours + - contour_colors, respective colors for the contours + """ + + # Calculate the contour levels + if llh_map: # likelihood map + min_value = np.min(equatorial_map) + if systematics: + # from Pan-Starrs event 127852 + # these are values determined from MC by Will on the TS (2*LLH) + # Not clear yet how to translate this for the probability map + contour_levels = (np.array([22.2, 64.2])+min_value) + contour_labels = [ + r'50% (IC160427A syst.)', r'90% (IC160427A syst.)' + ] + contour_colors = ['k', 'r'] + # Wilks + else: + contour_levels = ( + np.array([1.39, 4.61, 11.83, 28.74])+min_value + )[:3] + contour_labels = [ + r'50%', r'90%', r'3$\sigma$', r'5$\sigma$' + ][:3] + contour_colors = ['k', 'r', 'g', 'b'][:3] + else: # probability map + if systematics: + raise AssertionError( + "No corrected values for contours in probability maps" + ) + else: + sorted_values = np.sort(equatorial_map)[::-1] + probability_levels = ( + np.array([0.5, 0.9, 1-1.35e-3, 1-2.87e-7]) + )[:3] + contour_levels = list() + for prob in probability_levels: + level_index = ( + np.nancumsum(sorted_values) >= prob + ).tolist().index(True) + level = sorted_values[level_index] + contour_levels.append(level) + contour_labels = [r'50%', r'90%', r'3$\sigma$', r'5$\sigma$'][:3] + contour_colors = ['k', 'r', 'g', 'b'][:3] + + return contour_levels, contour_labels, contour_colors