Merge pull request #1208 from woodmd/master

Bug fix in gammapy.maps interpolation

gammapy/maps/geom.py

@@ -215,15 +215,15 @@ def axes_pix_to_coord(axes, pix):
     return coords
 
 
-def coord_to_idx(edges, x, bounded=False):
+def coord_to_idx(edges, x, clip=False):
     """Convert axis coordinates ``x`` to bin indices.
 
     Returns -1 for values below/above the lower/upper edge.
     """
     x = np.array(x, ndmin=1)
     ibin = np.digitize(x, edges) - 1
 
-    if bounded:
+    if clip:
         ibin[x < edges[0]] = 0
         ibin[x > edges[-1]] = len(edges) - 1
     else:
@@ -517,14 +517,14 @@ def coord_to_pix(self, coord):
         pix = coord_to_pix(self._nodes, coord, interp=self._interp)
         return np.array(pix + self._pix_offset, ndmin=1)
 
-    def coord_to_idx(self, coord, bounded=False):
+    def coord_to_idx(self, coord, clip=False):
         """Transform from axis coordinate to bin index.
 
         Parameters
         ----------
         coord : `~numpy.ndarray`
             Array of axis coordinate values.
-        bounded : bool
+        clip : bool
             Choose whether to clip the index to the valid range of the
             axis.  If false then indices for values outside the axis
             range will be set -1.
@@ -534,7 +534,7 @@ def coord_to_idx(self, coord, bounded=False):
         idx : `~numpy.ndarray`
             Array of bin indices.
         """
-        return coord_to_idx(self.edges, coord, bounded)
+        return coord_to_idx(self.edges, coord, clip)
 
     def coord_to_idx_interp(self, coord):
         """Compute indices of two nearest bins to the given coordinate.
@@ -545,8 +545,8 @@ def coord_to_idx_interp(self, coord):
             Array of axis coordinate values.
         """
 
-        return (coord_to_idx(self.center[:-1], coord, bounded=True),
-                coord_to_idx(self.center[:-1], coord, bounded=True) + 1,)
+        return (coord_to_idx(self.center[:-1], coord, clip=True),
+                coord_to_idx(self.center[:-1], coord, clip=True) + 1,)
 
     def slice(self, idx):
         """Create a new axis object by extracting a slice from this axis.
@@ -828,7 +828,7 @@ def coord_to_pix(self, coords):
         """
         pass
 
-    def coord_to_idx(self, coords):
+    def coord_to_idx(self, coords, clip=False):
         """Convert map coordinates to pixel indices.
 
         Parameters
@@ -839,16 +839,21 @@ def coord_to_idx(self, coords):
             If passed as a tuple then the ordering should be
             (longitude, latitude, c_0, ..., c_N) where c_i is the
             coordinate vector for axis i.
+        clip : bool
+            Choose whether to clip indices to the valid range of the
+            geometry.  If false then indices for coordinates outside
+            the geometry range will be set -1.
 
         Returns
         -------
         pix : tuple
             Tuple of pixel indices in image and band dimensions.
             Elements set to -1 correspond to coordinates outside the
             map.
+
         """
         pix = self.coord_to_pix(coords)
-        return self.pix_to_idx(pix)
+        return self.pix_to_idx(pix, clip=clip)
 
     @abc.abstractmethod
     def pix_to_coord(self, pix):
@@ -867,14 +872,18 @@ def pix_to_coord(self, pix):
         pass
 
     @abc.abstractmethod
-    def pix_to_idx(self, pix):
+    def pix_to_idx(self, pix, clip=False):
         """Convert pixel coordinates to pixel indices.  Returns -1 for pixel
         coordinates that lie outside of the map.
 
         Parameters
         ----------
         pix : tuple
             Tuple of pixel coordinates.
+        clip : bool
+            Choose whether to clip indices to the valid range of the
+            geometry.  If false then indices for coordinates outside
+            the geometry range will be set -1.
 
         Returns
         -------

gammapy/maps/hpx.py

@@ -660,16 +660,24 @@ def pix_to_coord(self, pix):
 
         return coords
 
-    def pix_to_idx(self, pix):
+    def pix_to_idx(self, pix, clip=False):
 
-        # FIXME: Correctly apply bounds on non-spatial pixel
-        # coordinates
+        # FIXME: Look for better method to clip HPX indices
         idx = list(pix_tuple_to_idx(pix, copy=True))
         idx_local = self.global_to_local(idx)
         for i, _ in enumerate(idx):
-            idx[i][(idx_local[i] < 0) | (idx[i] < 0)] = -1
-            if i > 0:
-                idx[i][idx[i] > self.axes[i - 1].nbin - 1] = -1
+
+            if clip:
+                if i > 0:
+                    np.clip(idx[i], 0, self.axes[i - 1].nbin - 1, out=idx[i])
+                else:
+                    np.clip(idx[i], 0, None, out=idx[i])
+            else:
+                if i > 0:
+                    np.putmask(idx[i],
+                               (idx[i] < 0) | (idx[i] >= self.axes[i - 1].nbin), -1)
+                else:
+                    np.putmask(idx[i], (idx_local[i] < 0) | (idx[i] < 0), -1)
 
         return tuple(idx)
 

gammapy/maps/hpxnd.py

@@ -312,8 +312,8 @@ def _interp_by_coords(self, coords, interp):
         import healpy as hp
 
         c = MapCoords.create(coords)
-        pix, wts = self._get_interp_weights(coords,
-                                            self.geom.coord_to_idx(c)[1:])
+        idx_ax = self.geom.coord_to_idx(c, clip=True)[1:]
+        pix, wts = self._get_interp_weights(coords, idx_ax)
 
         if self.geom.ndim == 2:
             return np.sum(self.data.T[pix] * wts, axis=0)
@@ -327,7 +327,7 @@ def _interp_by_coords(self, coords, interp):
             for j, ax in enumerate(self.geom.axes):
 
                 idx = coord_to_idx(ax.center[:-1],
-                                   c[2 + j], bounded=True)  # [None, ...]
+                                   c[2 + j], clip=True)  # [None, ...]
 
                 w = ax.center[idx + 1] - ax.center[idx]
                 if (i & (1 << j)):

gammapy/maps/tests/test_hpx.py

@@ -204,15 +204,19 @@ def test_hpxgeom_coord_to_idx(nside, nested, coordsys, region, axes):
     for i, z in enumerate(zidx):
         assert_allclose(z, idx[i + 1])
 
+    # Test w/ coords outside the geometry
     lon = np.array([0.0, 5.0, 10.0])
     lat = np.array([75.3, 75.3, 74.6])
     coords = make_test_coords(geom, lon, lat)
     zidx = [ax.coord_to_idx(t) for t, ax in zip(coords[2:], geom.axes)]
 
     idx = geom.coord_to_idx(coords)
     if geom.region is not None:
-        assert_allclose(-1 * np.ones(len(coords[0]), dtype=int),
-                        idx[0])
+        assert_allclose(np.full_like(coords[0], -1, dtype=int), idx[0])
+
+    idx = geom.coord_to_idx(coords, clip=True)
+    assert(np.all(np.not_equal(np.full_like(
+        coords[0], -1, dtype=int), idx[0])))
 
 
 def test_hpxgeom_coord_to_pix():
@@ -346,7 +350,7 @@ def test_hpxgeom_contains(nside, nested, coordsys, region, axes):
     geom = HpxGeom(nside, nested, coordsys, region=region, axes=axes)
     coords = geom.get_coords(flat=True)
     assert_allclose(geom.contains(coords),
-                    np.ones(coords[0].shape, dtype=bool))
+                    np.ones_like(coords[0], dtype=bool))
 
     if axes is not None:
         coords = [c[0] for c in coords[:2]] + \

gammapy/maps/tests/test_wcs.py

@@ -72,6 +72,16 @@ def test_wcsgeom_test_coord_to_idx(npix, binsz, coordsys, proj, skydir, axes):
     assert_allclose(geom.get_idx()[0],
                     geom.coord_to_idx(geom.get_coords())[0])
 
+    if not geom.is_allsky:
+        coords = geom.center_coord[:2] + \
+            tuple([ax.center[0] for ax in geom.axes])
+        coords[0][...] += 2.0 * np.max(geom.width[0])
+        idx = geom.coord_to_idx(coords)
+        assert_allclose(np.full_like(coords[0], -1, dtype=int), idx[0])
+        idx = geom.coord_to_idx(coords, clip=True)
+        assert(np.all(np.not_equal(np.full_like(
+            coords[0], -1, dtype=int), idx[0])))
+
 
 @pytest.mark.parametrize(('npix', 'binsz', 'coordsys', 'proj', 'skydir', 'axes'),
                          wcs_test_geoms)

gammapy/maps/wcs.py

@@ -516,24 +516,30 @@ def pix_to_coord(self, pix):
 
         return tuple(coords)
 
-    def pix_to_idx(self, pix):
+    def pix_to_idx(self, pix, clip=False):
         idxs = pix_tuple_to_idx(pix, copy=True)
         if not self.is_regular:
             ibin = [pix[2 + i] for i, ax in enumerate(self.axes)]
             ibin = pix_tuple_to_idx(ibin, copy=True)
-            m = np.all(np.stack([(t > 0) & (t < ax.nbin)
-                                 for t, ax in zip(idxs, self.axes)]), axis=0)
-            for i in range(len(ibin)):
-                ibin[i][~m] = 0
+            for i, ax in enumerate(self.axes):
+                np.clip(ibin[i], 0, ax.nbin - 1, out=ibin[i])
             npix = (self.npix[0][ibin], self.npix[1][ibin])
         else:
             npix = self.npix
 
         for i, idx in enumerate(idxs):
-            if i < 2:
-                idxs[i][(idx < 0) | (idx >= npix[i])] = -1
+            if clip:
+                if i < 2:
+                    np.clip(idxs[i], 0, npix[i], out=idxs[i])
+                else:
+                    np.clip(idxs[i], 0, self.axes[i - 2].nbin - 1, out=idxs[i])
             else:
-                idxs[i][(idx < 0) | (idx >= self.axes[i - 2].nbin)] = -1
+                if i < 2:
+                    np.putmask(idxs[i], (idx < 0) | (idx >= npix[i]), -1)
+                else:
+                    np.putmask(idxs[i], (idx < 0) | (
+                        idx >= self.axes[i - 2].nbin), -1)
+
         return idxs
 
     def contains(self, coords):