set nan to min values only if convolution is applied, or to contour c…

…alculation or saving of flatted probability map
icecube · Dec 10, 2024 · 8a47317 · 8a47317
1 parent a482ffd
commit 8a47317
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 5 deletions.
diff --git a/skyreader/plot/plot.py b/skyreader/plot/plot.py
@@ -448,16 +448,21 @@ def bounding_box(ra, dec, theta, phi):
         sample_points = np.array([np.pi/2 - grid_dec, grid_ra]).T
 
         if not circular:
+            grid_values_for_contours = copy.copy(grid_value)
             if llh_map:
                 # get rid of nan values only for the contours
                 # this avoids crashes during plotting
-                grid_values_for_contours = copy.copy(grid_value)
                 grid_values_for_contours[
                     np.isnan(grid_values_for_contours)
                 ] = np.nanmax(equatorial_map)
                 contour_levels_for_contours = contour_levels
             else:
-                grid_values_for_contours = np.log(grid_value)
+                grid_values_for_contours[
+                    np.isnan(grid_values_for_contours)
+                ] = np.nanmin(equatorial_map)
+                grid_values_for_contours = np.log(
+                    grid_values_for_contours
+                )
                 contour_levels_for_contours = np.log(contour_levels)
             # Get contours from healpix map
             contours_by_level = meander.spherical_contours(
@@ -827,6 +832,9 @@ def bounding_box(ra, dec, theta, phi):
             column_names = ['2DLLH']
         else:
             # avoid excessively heavy data format for the flattened map
+            equatorial_map[np.isnan(equatorial_map)] = np.nanmin(
+                equatorial_map
+            )
             equatorial_map[equatorial_map < 1e-16] = np.mean(
                 equatorial_map[equatorial_map < 1e-16]
             )

diff --git a/skyreader/utils/handle_map_data.py b/skyreader/utils/handle_map_data.py
@@ -116,14 +116,16 @@ def extract_map(
         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
+        # if the probability is less than ~1e-300 this is set to 0
+        # because of the float64 data format
         min_map = np.nanmin(equatorial_map[equatorial_map > 0.])
-        equatorial_map[np.isnan(equatorial_map)] = min_map
         equatorial_map[equatorial_map == 0.] = min_map
 
         if angular_error_floor is not None:
             # convolute with a gaussian. angular_error_floor is the
             # sigma in deg.
+            # nan values are a problem for the convolution and the contours
+            equatorial_map[np.isnan(equatorial_map)] = min_map
             equatorial_map = healpy.smoothing(
                 equatorial_map,
                 sigma=np.deg2rad(angular_error_floor),
@@ -140,7 +142,6 @@ def extract_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, uniq_array