Merge pull request #285 from pnuu/bugfix-kdtree-segfault

Check the source area orientation before defining slices

pyresample/geometry.py

@@ -711,9 +711,9 @@ def _compute_omerc_parameters(self, ellipsoid):
 
         # We need to compute alpha-based omerc for geotiff support
         lonc, lat0 = Proj(**proj_dict2points)(0, 0, inverse=True)
-        az1, az2, dist = Geod(**proj_dict2points).inv(lonc, lat0, lon2, lat2)
+        az1, az2, _ = Geod(**proj_dict2points).inv(lonc, lat0, lon2, lat2)
         azimuth = az1
-        az1, az2, dist = Geod(**proj_dict2points).inv(lonc, lat0, lon1, lat1)
+        az1, az2, _ = Geod(**proj_dict2points).inv(lonc, lat0, lon1, lat1)
         if abs(az1 - azimuth) > 1:
             if abs(az2 - azimuth) > 1:
                 logger.warning("Can't find appropriate azimuth.")
@@ -1152,7 +1152,7 @@ def __init__(self, area_id, description, proj_id, projection, width, height,
 
     @property
     def _crs(self):
-        """Helper property for the `crs` property.
+        """Wrap the `crs` property in a helper property.
 
         The :class:`pyproj.crs.CRS` object is not thread-safe. To avoid
         accidentally passing it between threads, we only create it when it
@@ -1982,6 +1982,26 @@ def proj4_string(self):
                       "instead.", DeprecationWarning)
         return proj4_dict_to_str(self.proj_dict)
 
+    def _get_slice_starts_stops(self, area_to_cover):
+        """Get x and y start and stop points for slicing."""
+        llx, lly, urx, ury = area_to_cover.area_extent
+        x, y = self.get_xy_from_proj_coords([llx, urx], [lly, ury])
+
+        if self.area_extent[0] > self.area_extent[2]:
+            xstart = 0 if x[1] is np.ma.masked else x[1]
+            xstop = self.width if x[0] is np.ma.masked else x[0] + 1
+        else:
+            xstart = 0 if x[0] is np.ma.masked else x[0]
+            xstop = self.width if x[1] is np.ma.masked else x[1] + 1
+        if self.area_extent[1] > self.area_extent[3]:
+            ystart = 0 if y[0] is np.ma.masked else y[0]
+            ystop = self.height if y[1] is np.ma.masked else y[1] + 1
+        else:
+            ystart = 0 if y[1] is np.ma.masked else y[1]
+            ystop = self.height if y[0] is np.ma.masked else y[0] + 1
+
+        return xstart, xstop, ystart, ystop
+
     def get_area_slices(self, area_to_cover, shape_divisible_by=None):
         """Compute the slice to read based on an `area_to_cover`."""
         if not isinstance(area_to_cover, AreaDefinition):
@@ -1994,14 +2014,9 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None):
             logger.debug('Projections for data and slice areas are'
                          ' identical: %s',
                          proj_def_to_cover)
-            # Get xy coordinates
-            llx, lly, urx, ury = area_to_cover.area_extent
-            x, y = self.get_xy_from_proj_coords([llx, urx], [lly, ury])
-
-            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
+            # Get slice parameters
+            xstart, xstop, ystart, ystop = self._get_slice_starts_stops(
+                area_to_cover)
 
             return (check_slice_orientation(slice(xstart, xstop)),
                     check_slice_orientation(slice(ystart, ystop)))

pyresample/test/test_geometry.py

@@ -203,7 +203,6 @@ def test_create_areas_def(self):
             self.assertEqual(res['projection'][proj_key],
                              expected['projection'][proj_key])
 
-
         # EPSG
         projections = {'+init=epsg:3006': 'init: epsg:3006'}
         if utils.is_pyproj2():
@@ -1001,6 +1000,71 @@ def test_get_xy_from_lonlat(self):
         self.assertTrue((x__.data == x_expects).all())
         self.assertTrue((y__.data == y_expects).all())
 
+    def test_get_slice_starts_stops(self):
+        """Check area slice end-points."""
+        from pyresample import utils
+        area_id = 'orig'
+        area_name = 'Test area'
+        proj_id = 'test'
+        x_size = 3712
+        y_size = 3712
+        area_extent = (-5570248.477339745, -5561247.267842293, 5567248.074173927, 5570248.477339745)
+        proj_dict = {'a': 6378169.0, 'b': 6356583.8, 'h': 35785831.0,
+                     'lon_0': 0.0, 'proj': 'geos', 'units': 'm'}
+        target_area = utils.get_area_def(area_id,
+                                         area_name,
+                                         proj_id,
+                                         proj_dict,
+                                         x_size, y_size,
+                                         area_extent)
+
+        # Expected result is the same for all cases
+        expected = (3, 3709, 3, 3709)
+
+        # Source and target have the same orientation
+        area_extent = (-5580248.477339745, -5571247.267842293, 5577248.074173927, 5580248.477339745)
+        source_area = utils.get_area_def(area_id,
+                                         area_name,
+                                         proj_id,
+                                         proj_dict,
+                                         x_size, y_size,
+                                         area_extent)
+        res = source_area._get_slice_starts_stops(target_area)
+        assert res == expected
+
+        # Source is flipped in X direction
+        area_extent = (5577248.074173927, -5571247.267842293, -5580248.477339745, 5580248.477339745)
+        source_area = utils.get_area_def(area_id,
+                                         area_name,
+                                         proj_id,
+                                         proj_dict,
+                                         x_size, y_size,
+                                         area_extent)
+        res = source_area._get_slice_starts_stops(target_area)
+        assert res == expected
+
+        # Source is flipped in Y direction
+        area_extent = (-5580248.477339745, 5580248.477339745, 5577248.074173927, -5571247.267842293)
+        source_area = utils.get_area_def(area_id,
+                                         area_name,
+                                         proj_id,
+                                         proj_dict,
+                                         x_size, y_size,
+                                         area_extent)
+        res = source_area._get_slice_starts_stops(target_area)
+        assert res == expected
+
+        # Source is flipped in both X and Y directions
+        area_extent = (5577248.074173927, 5580248.477339745, -5580248.477339745, -5571247.267842293)
+        source_area = utils.get_area_def(area_id,
+                                         area_name,
+                                         proj_id,
+                                         proj_dict,
+                                         x_size, y_size,
+                                         area_extent)
+        res = source_area._get_slice_starts_stops(target_area)
+        assert res == expected
+
     def test_get_area_slices(self):
         """Check area slicing."""
         from pyresample import utils