Merge pull request #123 from djhoese/bugfix-area-slicing

Fix get_area_slices edge case when target is larger than source

docs/source/swath.rst

@@ -301,7 +301,7 @@ Function for resampling using bilinear interpolation for irregular source grids.
  ...                                     reduce_data=True, segments=None,
  ...                                     epsilon=0)
 
-The **target_area** needs to be an area definition with **proj4_string**
+The **target_area** needs to be an area definition with **proj_str**
 attribute.
 
 ..

pyresample/bilinear/__init__.py

@@ -235,7 +235,7 @@ def get_bil_info(source_geo_def, target_area_def, radius=50e3, neighbours=32,
     idx_ref = np.where(index_mask, 0, idx_ref)
 
     # Get output projection as pyproj object
-    proj = Proj(target_area_def.proj4_string)
+    proj = Proj(target_area_def.proj_str)
 
     # Get output x/y coordinates
     out_x, out_y = _get_output_xy(target_area_def, proj)

pyresample/bilinear/xarr.py

@@ -126,7 +126,7 @@ def get_bil_info(self):
         index_array = da.where(index_mask, 0, index_array)
 
         # Get output projection as pyproj object
-        proj = Proj(self.target_geo_def.proj4_string)
+        proj = Proj(self.target_geo_def.proj_str)
 
         # Get output x/y coordinates
         out_x, out_y = _get_output_xy_dask(self.target_geo_def, proj)

pyresample/ewa/__init__.py

@@ -115,7 +115,7 @@ def ll2cr(swath_def, area_def, fill=np.nan, copy=True):
         lats = lats.astype(np.float64)
 
     # Break the input area up in to the expected parameters for ll2cr
-    p = area_def.proj4_string
+    p = area_def.proj_str
     cw = area_def.pixel_size_x
     # cell height must be negative for this to work as expected
     ch = -abs(area_def.pixel_size_y)

pyresample/geometry.py

@@ -81,6 +81,15 @@ def __init__(self, lons=None, lats=None, nprocs=1):
         self.cartesian_coords = None
         self.hash = None
 
+    def __getitem__(self, key):
+        """Slice a 2D geographic definition."""
+        y_slice, x_slice = key
+        return self.__class__(
+            lons=self.lons[y_slice, x_slice],
+            lats=self.lats[y_slice, x_slice],
+            nprocs=self.nprocs
+        )
+
     def __hash__(self):
         """Compute the hash of this object."""
         if self.hash is None:
@@ -978,7 +987,7 @@ def colrow2lonlat(self, cols, rows):
         To be used with scarse data points instead of slices
         (see get_lonlats).
         """
-        p = Proj(self.proj4_string)
+        p = Proj(self.proj_str)
         x = self.projection_x_coords
         y = self.projection_y_coords
         return p(y[y.size - cols], x[x.size - rows], inverse=True)
@@ -1024,7 +1033,7 @@ def get_xy_from_lonlat(self, lon, lat):
             if lon.shape != lat.shape:
                 raise ValueError("lon and lat is not of the same shape!")
 
-        pobj = Proj(self.proj4_string)
+        pobj = Proj(self.proj_str)
         xm_, ym_ = pobj(lon, lat)
 
         return self.get_xy_from_proj_coords(xm_, ym_)
@@ -1078,15 +1087,15 @@ def get_xy_from_proj_coords(self, xm, ym):
         y__ = (ym - upl_y) / yscale
 
         if isinstance(x__, np.ndarray) and isinstance(y__, np.ndarray):
-            mask = (((x__ < 0) | (x__ > self.x_size)) |
-                    ((y__ < 0) | (y__ > self.y_size)))
+            mask = (((x__ < 0) | (x__ >= self.x_size)) |
+                    ((y__ < 0) | (y__ >= self.y_size)))
             return (np.ma.masked_array(x__.astype('int'), mask=mask,
                                        fill_value=-1, copy=False),
                     np.ma.masked_array(y__.astype('int'), mask=mask,
                                        fill_value=-1, copy=False))
         else:
-            if ((x__ < 0 or x__ > self.x_size) or
-                    (y__ < 0 or y__ > self.y_size)):
+            if ((x__ < 0 or x__ >= self.x_size) or
+                    (y__ < 0 or y__ >= self.y_size)):
                 raise ValueError('Point outside area:( %f %f)' % (x__, y__))
             return int(x__), int(y__)
 
@@ -1324,6 +1333,8 @@ def get_lonlats(self, nprocs=None, data_slice=None, cache=False, dtype=None):
     @property
     def proj4_string(self):
         """Return projection definition as Proj.4 string."""
+        warnings.warn("'proj4_string' is deprecated, please use 'proj_str' "
+                      "instead.", DeprecationWarning)
         return utils.proj4_dict_to_str(self.proj_dict)
 
     def get_area_slices(self, area_to_cover):
@@ -1335,7 +1346,7 @@ def get_area_slices(self, area_to_cover):
             raise NotImplementedError('Only geos supported')
 
         # Intersection only required for two different projections
-        if area_to_cover.proj_dict.get('proj') == self.proj_dict['proj']:
+        if area_to_cover.proj_str == self.proj_str:
             logger.debug('Projections for data and slice areas are'
                          ' identical: %s', area_to_cover.proj_dict['proj'])
             # Get xy coordinates
@@ -1350,8 +1361,12 @@ def get_area_slices(self, area_to_cover):
             return slice(xstart, xstop), slice(ystart, ystop)
 
         data_boundary = Boundary(*get_geostationary_bounding_box(self))
+        if area_to_cover.proj_dict['proj'] == 'geos':
+            area_boundary = Boundary(
+                *get_geostationary_bounding_box(area_to_cover))
+        else:
+            area_boundary = AreaDefBoundary(area_to_cover, 100)
 
-        area_boundary = AreaDefBoundary(area_to_cover, 100)
         intersection = data_boundary.contour_poly.intersection(
             area_boundary.contour_poly)
         if intersection is None:
@@ -1360,7 +1375,8 @@ def get_area_slices(self, area_to_cover):
         x, y = self.get_xy_from_lonlat(np.rad2deg(intersection.lon),
                                        np.rad2deg(intersection.lat))
 
-        return slice(min(x), max(x) + 1), slice(min(y), max(y) + 1)
+        return (slice(np.ma.min(x), np.ma.max(x) + 1),
+                slice(np.ma.min(y), np.ma.max(y) + 1))
 
     def crop_around(self, other_area):
         """Crop this area around `other_area`."""
@@ -1418,8 +1434,8 @@ def get_geostationary_bounding_box(geos_area, nb_points=50):
 
     # generate points around the north hemisphere in satellite projection
     # make it a bit smaller so that we stay inside the valid area
-    x = np.cos(np.linspace(-np.pi, 0, nb_points / 2)) * (xmax - 0.0001)
-    y = -np.sin(np.linspace(-np.pi, 0, nb_points / 2)) * (ymax - 0.0001)
+    x = np.cos(np.linspace(-np.pi, 0, int(nb_points / 2))) * (xmax - 0.0001)
+    y = -np.sin(np.linspace(-np.pi, 0, int(nb_points / 2))) * (ymax - 0.0001)
 
     ll_x, ll_y, ur_x, ur_y = geos_area.area_extent
 
@@ -1577,7 +1593,14 @@ def squeeze(self):
     @property
     def proj4_string(self):
         """Returns projection definition as Proj.4 string"""
-        return self.defs[0].proj4_string
+        warnings.warn("'proj4_string' is deprecated, please use 'proj_str' "
+                      "instead.", DeprecationWarning)
+        return self.defs[0].proj_str
+
+    @property
+    def proj_str(self):
+        """Returns projection definition as Proj.4 string"""
+        return self.defs[0].proj_str
 
     def update_hash(self, the_hash=None):
         for areadef in self.defs:

pyresample/test/test_bilinear.py

@@ -154,19 +154,19 @@ def test_solve_quadratic(self):
         self.assertAlmostEqual(res[0], 0.5, 5)
 
     def test_get_output_xy(self):
-        proj = Proj(self.target_def.proj4_string)
+        proj = Proj(self.target_def.proj_str)
         out_x, out_y = bil._get_output_xy(self.target_def, proj)
         self.assertTrue(out_x.all())
         self.assertTrue(out_y.all())
 
     def test_get_input_xy(self):
-        proj = Proj(self.target_def.proj4_string)
+        proj = Proj(self.target_def.proj_str)
         in_x, in_y = bil._get_output_xy(self.swath_def, proj)
         self.assertTrue(in_x.all())
         self.assertTrue(in_y.all())
 
     def test_get_bounding_corners(self):
-        proj = Proj(self.target_def.proj4_string)
+        proj = Proj(self.target_def.proj_str)
         out_x, out_y = bil._get_output_xy(self.target_def, proj)
         in_x, in_y = bil._get_input_xy(self.swath_def, proj,
                                        self.input_idxs, self.idx_ref)

pyresample/test/test_geometry.py

@@ -27,16 +27,6 @@ class Test(unittest.TestCase):
 
     """Unit testing the geometry and geo_filter modules"""
 
-    def assert_raises(self, exception, call_able, *args):
-        """assertRaises() has changed from py2.6 to 2.7! Here is an attempt to
-        cover both"""
-        import sys
-        if sys.version_info < (2, 7):
-            self.assertRaises(exception, call_able, *args)
-        else:
-            with self.assertRaises(exception):
-                call_able(*args)
-
     def test_lonlat_precomp(self):
         area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
                                            {'a': '6378144.0',
@@ -730,8 +720,8 @@ def test_get_xy_from_lonlat(self):
         self.assertEqual(y__, 0)
 
         lon, lat = p__(999000, -10, inverse=True)
-        self.assert_raises(ValueError, area_def.get_xy_from_lonlat, lon, lat)
-        self.assert_raises(ValueError, area_def.get_xy_from_lonlat, 0., 0.)
+        self.assertRaises(ValueError, area_def.get_xy_from_lonlat, lon, lat)
+        self.assertRaises(ValueError, area_def.get_xy_from_lonlat, 0., 0.)
 
         # Test getting arrays back:
         lons = [lon_ll + eps_lonlat, lon_ur - eps_lonlat]
@@ -746,22 +736,8 @@ def test_get_xy_from_lonlat(self):
     def test_get_area_slices(self):
         """Check area slicing."""
         from pyresample import utils
-        area_id = 'cover'
-        area_name = 'Area to cover'
-        proj_id = 'test'
-        x_size = 3712
-        y_size = 3712
-        area_extent = (-5570248.477339261, -5567248.074173444, 5567248.074173444, 5570248.477339261)
-        proj_dict = {'a': 6378169.5, 'b': 6356583.8, 'h': 35785831.0,
-                     'lon_0': 0.0, 'proj': 'geos', 'units': 'm'}
-
-        area_to_cover = utils.get_area_def(area_id,
-                                           area_name,
-                                           proj_id,
-                                           proj_dict,
-                                           x_size, y_size,
-                                           area_extent)
 
+        # The area of our source data
         area_id = 'orig'
         area_name = 'Test area'
         proj_id = 'test'
@@ -777,9 +753,41 @@ def test_get_area_slices(self):
                                       x_size, y_size,
                                       area_extent)
 
+        # An area that is a subset of the original one
+        area_to_cover = utils.get_area_def(
+            'cover_subset',
+            'Area to cover',
+            'test',
+            proj_dict,
+            1000, 1000,
+            area_extent=(area_extent[0] + 10000,
+                         area_extent[1] + 10000,
+                         area_extent[2] - 10000,
+                         area_extent[3] - 10000))
         slice_x, slice_y = area_def.get_area_slices(area_to_cover)
-        self.assertEqual(slice_x, slice(0, 3712))
-        self.assertEqual(slice_y, slice(0, 3712))
+        self.assertEqual(slice(3, 3709, None), slice_x)
+        self.assertEqual(slice(3, 3709, None), slice_y)
+
+        # An area similar to the source data but not the same
+        area_id = 'cover'
+        area_name = 'Area to cover'
+        proj_id = 'test'
+        x_size = 3712
+        y_size = 3712
+        area_extent = (-5570248.477339261, -5567248.074173444, 5567248.074173444, 5570248.477339261)
+        proj_dict = {'a': 6378169.5, 'b': 6356583.8, 'h': 35785831.0,
+                     'lon_0': 0.0, 'proj': 'geos', 'units': 'm'}
+
+        area_to_cover = utils.get_area_def(area_id,
+                                           area_name,
+                                           proj_id,
+                                           proj_dict,
+                                           x_size, y_size,
+                                           area_extent)
+        slice_x, slice_y = area_def.get_area_slices(area_to_cover)
+        self.assertEqual(slice(46, 3667, None), slice_x)
+        self.assertEqual(slice(52, 3663, None), slice_y)
+
 
         area_to_cover = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
                                                 {'a': 6378144.0,
@@ -798,7 +806,7 @@ def test_get_area_slices(self):
         self.assertEqual(slice_x, slice(1610, 2343))
         self.assertEqual(slice_y, slice(158, 515, None))
 
-    def test_proj4_string(self):
+    def test_proj_str(self):
         from collections import OrderedDict
         proj_dict = OrderedDict()
         proj_dict['proj'] = 'stere'
@@ -811,11 +819,11 @@ def test_proj4_string(self):
                                        proj_dict, 10, 10,
                                        [-1370912.72, -909968.64, 1029087.28,
                                         1490031.36])
-        self.assertEquals(area.proj4_string,
-                          '+proj=stere +a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0')
+        self.assertEquals(area.proj_str,
+                          '+a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0 +proj=stere')
         proj_dict['no_rot'] = ''
-        self.assertEquals(area.proj4_string,
-                          '+proj=stere +a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0 +no_rot')
+        self.assertEquals(area.proj_str,
+                          '+a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0 +no_rot +proj=stere')
 
 
 def assert_np_dict_allclose(dict1, dict2):
@@ -844,6 +852,17 @@ def test_swath(self):
         self.assertFalse(id(lons1) != id(lons2) or id(lats1) != id(lats2),
                          msg='Caching of swath coordinates failed')
 
+    def test_slice(self):
+        """Test that SwathDefinitions can be sliced."""
+        lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
+        lats1 = np.fromfunction(
+            lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
+
+        swath_def = geometry.SwathDefinition(lons1, lats1)
+        new_swath_def = swath_def[1000:4000, 20:40]
+        self.assertTupleEqual(new_swath_def.lons.shape, (3000, 20))
+        self.assertTupleEqual(new_swath_def.lats.shape, (3000, 20))
+
     def test_concat_1d(self):
         lons1 = np.array([1, 2, 3])
         lats1 = np.array([1, 2, 3])

pyresample/test/test_grid.py

@@ -201,7 +201,7 @@ def test_single_lonlat(self):
             lat, 52.566998432390619, msg='Resampling of single lat failed')
 
     def test_proj4_string(self):
-        proj4_string = self.area_def.proj4_string
+        proj4_string = self.area_def.proj_str
         expected_string = '+a=6378144.0 +b=6356759.0 +lat_ts=50.00 +lon_0=8.00 +proj=stere +lat_0=50.00'
         self.assertEqual(
             frozenset(proj4_string.split()), frozenset(expected_string.split()))