Merge pull request #277 from pnuu/debug-kdtree-segfault

Fix calculating area slices for flipped projections

pyresample/geometry.py

@@ -1961,9 +1961,14 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None):
 
             xstart = 0 if x[0] is np.ma.masked else x[0]
             ystart = 0 if y[1] is np.ma.masked else y[1]
-            xstop = self.width if x[1] is np.ma.masked else x[1] + 1
-            ystop = self.height if y[0] is np.ma.masked else y[0] + 1
-
+            if self.area_extent[0] > self.area_extent[2]:
+                xstop = self.width if x[0] is np.ma.masked else x[0] + 1
+            else:
+                xstop = self.width if x[1] is np.ma.masked else x[1] + 1
+            if self.area_extent[1] > self.area_extent[3]:
+                ystop = self.height if y[1] is np.ma.masked else y[1] + 1
+            else:
+                ystop = self.height if y[0] is np.ma.masked else y[0] + 1
             return slice(xstart, xstop), slice(ystart, ystop)
 
         if self.proj_dict.get('proj') != 'geos':
@@ -1986,7 +1991,6 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None):
             raise NotImplementedError
         x, y = self.get_xy_from_lonlat(np.rad2deg(intersection.lon),
                                        np.rad2deg(intersection.lat))
-
         x_slice = slice(np.ma.min(x), np.ma.max(x) + 1)
         y_slice = slice(np.ma.min(y), np.ma.max(y) + 1)
         if shape_divisible_by is not None:

pyresample/test/test_geometry.py

@@ -1073,6 +1073,26 @@ def test_get_area_slices(self):
         self.assertEqual(slice_x, slice(1610, 2343))
         self.assertEqual(slice_y, slice(158, 515, None))
 
+        # The same as source area, but flipped in X and Y
+        area_id = 'cover'
+        area_name = 'Area to cover'
+        proj_id = 'test'
+        x_size = 3712
+        y_size = 3712
+        area_extent = (5567248.074173927, 5570248.477339745, -5570248.477339745, -5561247.267842293)
+        proj_dict = {'a': 6378169.0, '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(0, x_size, None), slice_x)
+        self.assertEqual(slice(0, y_size, None), slice_y)
+
         # totally different area
         projections = [{"init": 'EPSG:4326'}]
         if utils.is_pyproj2():