Merge pull request #170 from ortk95/dev

Propagate NaN values when mapping with spline interpolation

planetmapper/body_xy.py

@@ -906,7 +906,8 @@ def map_img(
         img: np.ndarray,
         degree_interval: float = 1,
         interpolation: Literal['nearest', 'linear', 'quadratic', 'cubic'] = 'linear',
-        warn_nan: bool = True,
+        propagate_nan: bool = True,
+        warn_nan: bool = False,
     ) -> np.ndarray:
         """
         Project an observed image onto a lon/lat grid.
@@ -929,6 +930,13 @@ def map_img(
             interpolation: Interpolation used when mapping. This can either any of
                 `'nearest'`, `'linear'`, `'quadratic'` or `'cubic'`. The default is
                 `'linear'`.
+            propagate_nan: If using spline interpolation, propagate NaN values from the
+                image to the mapped data. If `propagate_nan` is `True` (the default),
+                the interpolation is performed as normal (i.e. with NaN values in the
+                image set to 0), then any mapped locations where the nearest
+                corresponding image pixel is NaN are set to NaN. Note that there may
+                still be very small errors on the boundaries of NaN regions caused by
+                the interpolation.
             warn_nan: Print warning if any values in `img` are NaN when any of the
                 spline interpolations are used.
 
@@ -946,9 +954,9 @@ def map_img(
             'quadratic': 2,
             'cubic': 3,
         }
-        nan_sentinel = -999
 
         if interpolation == 'nearest':
+            nan_sentinel = -999
             x_map = np.asarray(
                 np.nan_to_num(np.round(x_map), nan=nan_sentinel), dtype=int
             )
@@ -963,6 +971,7 @@ def map_img(
                 projected[a, b] = img[y, x]
         elif interpolation in spline_k:
             k = spline_k[interpolation]
+            nans = np.isnan(img)
             if np.any(np.isnan(img)):
                 if warn_nan:
                     print('Warning, image contains NaN values which will be set to 0')
@@ -975,11 +984,16 @@ def map_img(
                 if math.isnan(x):
                     continue
                 y = y_map[a, b]  # y should never be nan when x is not nan
+                if propagate_nan and nans[int(np.round(y)), int(np.round(x))]:
+                    continue
                 projected[a, b] = interpolator(y, x)
         else:
             raise ValueError(f'Unknown interpolation method {interpolation!r}')
         return projected
 
+    def _xy_in_image_frame(self, x: float, y: float) -> bool:
+        return (-0.5 < x < self._nx - 0.5) and (-0.5 < y < self._ny - 0.5)
+
     # Backplane management
     @staticmethod
     def standardise_backplane_name(name: str) -> str:
@@ -1552,7 +1566,7 @@ def _get_xy_map(self, degree_interval: float) -> np.ndarray:
             ra, dec = radec_map[a, b]
             if not math.isnan(ra):
                 x, y = self.radec2xy(ra, dec)
-                if (-0.5 < x < self._nx - 0.5) and (-0.5 < y < self._ny - 0.5):
+                if self._xy_in_image_frame(x, y):
                     out[a, b] = x, y
         return out
 

planetmapper/common.py

@@ -1,3 +1,3 @@
-__version__ = '1.4.4'
+__version__ = '1.4.5'
 __author__ = 'Oliver King'
 __url__ = 'https://github.com/ortk95/planetmapper'

planetmapper/observation.py

@@ -492,6 +492,7 @@ def get_mapped_data(
         self,
         degree_interval: float = 1,
         interpolation: Literal['nearest', 'linear', 'quadratic', 'cubic'] = 'linear',
+        **kw,
     ) -> np.ndarray:
         """
         Projects the observed :attr:`data` onto a lon/lat grid using
@@ -510,6 +511,7 @@ def get_mapped_data(
             interpolation: Interpolation used when mapping. This can either any of
                 `'nearest'`, `'linear'`, `'quadratic'` or `'cubic'`. Passed to
                 :func:`BodyXY.map_img`.
+            **kw: Additional arguments passed to :func:`BodyXY.map_img`.
 
         Returns:
             Array containing a cube of cylindrical map of the values in :attr:`data` at
@@ -518,7 +520,7 @@ def get_mapped_data(
         """
         # Return a copy so that the cached value isn't tainted by any modifications
         return self._get_mapped_data(
-            degree_interval=degree_interval, interpolation=interpolation
+            degree_interval=degree_interval, interpolation=interpolation, **kw
         ).copy()
 
     @_cache_clearable_result_with_args
@@ -527,6 +529,7 @@ def _get_mapped_data(
         self,
         degree_interval: float = 1,
         interpolation: Literal['nearest', 'linear', 'quadratic', 'cubic'] = 'linear',
+        **kw,
     ):
         projected = []
         if interpolation == 'linear' and np.any(np.isnan(self.data)):
@@ -540,7 +543,7 @@ def _get_mapped_data(
             self._update_progress_hook(idx / len(data))
             projected.append(
                 self.map_img(
-                    img, degree_interval=degree_interval, interpolation=interpolation
+                    img, degree_interval=degree_interval, interpolation=interpolation, **kw
                 )
             )
         return np.array(projected)

scratchpad.py

@@ -3,16 +3,18 @@
 """Script for testing stuff during development"""
 import planetmapper
 import matplotlib.pyplot as plt
-import astropy.wcs
-import numpy as np
-import spiceypy as spice
 
 p = '/Users/ortk1/Dropbox/PhD/data/jwst/Uranus_2023jan08/lon2/stage3/d1/Level3_ch4-short_s3d.fits'
 p = '/Users/ortk1/Dropbox/PhD/data/jwst/saturn/SATURN-75N/stage6_flat/d1_fringe_nav/Level3_ch1-short_s3d_nav.fits'
 body = planetmapper.Observation(p)
-body.plot_backplane_map('pixel-x')
-plt.show()
+# body.plot_backplane_map('pixel-x')
+# plt.show()
 
 img = body.data[0]
 mapped = body.map_img(img, interpolation='cubic')
 body.imshow_map(mapped)
+plt.show()
+
+# mapped = body.map_img(img, interpolation='linear', mask_nan=False)
+# body.imshow_map(mapped)
+# plt.show()