Merge pull request #125 from djhoese/bugfix-equal-area-slices

Fix area slices not working for non-geos projections

.travis.yml

@@ -29,7 +29,7 @@ script:
 - coverage run --source=pyresample setup.py test && cd docs && mkdir doctest && sphinx-build
   -E -n -b doctest ./source ./doctest && cd ..
 after_success:
-- if [[ $TRAVIS_PYTHON_VERSION == 3.6 ]]; then coveralls; codecov; fi
+- if [[ $PYTHON_VERSION == 3.6 ]]; then coveralls; codecov; fi
 deploy:
   - provider: pypi
     user: dhoese

pyresample/boundary.py

@@ -106,7 +106,7 @@ class AreaDefBoundary(AreaBoundary):
     def __init__(self, area, frequency=1):
         lons, lats = area.get_bbox_lonlats()
         AreaBoundary.__init__(self,
-                              *zip(*(lons, lats)))
+                              *zip(lons, lats))
 
         if frequency != 1:
             self.decimate(frequency)

pyresample/geometry.py

@@ -556,7 +556,6 @@ def _compute_omerc_parameters(self, ellipsoid):
                             'no_rot': True
                             }
 
-        # return proj_dict2points
         # 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)
@@ -1342,9 +1341,6 @@ def get_area_slices(self, area_to_cover):
         if not isinstance(area_to_cover, AreaDefinition):
             raise NotImplementedError('Only AreaDefinitions can be used')
 
-        if self.proj_dict.get('proj') != 'geos':
-            raise NotImplementedError('Only geos supported')
-
         # Intersection only required for two different projections
         if area_to_cover.proj_str == self.proj_str:
             logger.debug('Projections for data and slice areas are'
@@ -1360,6 +1356,11 @@ def get_area_slices(self, area_to_cover):
 
             return slice(xstart, xstop), slice(ystart, ystop)
 
+        if self.proj_dict.get('proj') != 'geos':
+            raise NotImplementedError("Source projection must be 'geos' if "
+                                      "source/target projections are not "
+                                      "equal.")
+
         data_boundary = Boundary(*get_geostationary_bounding_box(self))
         if area_to_cover.proj_dict['proj'] == 'geos':
             area_boundary = Boundary(

pyresample/test/test_geometry.py

@@ -806,6 +806,41 @@ def test_get_area_slices(self):
         self.assertEqual(slice_x, slice(1610, 2343))
         self.assertEqual(slice_y, slice(158, 515, None))
 
+    def test_get_area_slices_nongeos(self):
+        """Check area slicing for non-geos projections."""
+        from pyresample import utils
+
+        # The area of our source data
+        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, 'lat_1': 25.,
+                     'lat_2': 25., 'lon_0': 0.0, 'proj': 'lcc', 'units': 'm'}
+        area_def = utils.get_area_def(area_id,
+                                      area_name,
+                                      proj_id,
+                                      proj_dict,
+                                      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(3, 3709, None), slice_x)
+        self.assertEqual(slice(3, 3709, None), slice_y)
+
     def test_proj_str(self):
         from collections import OrderedDict
         proj_dict = OrderedDict()