Merge pull request #424 from mraspaud/dynamic-resolution

Fix DynamicAreaDefinition resolution handling for incomplete projection definitions

pyresample/area_config.py

@@ -477,7 +477,8 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No
         p = Proj(crs, preserve_units=True)
     except (RuntimeError, CRSError):
         # Assume that an invalid projection will be "fixed" by a dynamic area definition later
-        return _make_area(area_id, description, proj_id, projection, shape, area_extent, **kwargs)
+        return _make_area(area_id, description, proj_id, projection, shape, area_extent,
+                          resolution=resolution, **kwargs)
 
     # If no units are provided, try to get units used in proj_dict. If still none are provided, use meters.
     if units is None:

pyresample/geometry.py

@@ -827,7 +827,7 @@ def compute_bb_proj_params(self, proj_dict):
             new_proj.update(proj_dict)
             return new_proj
 
-    def _compute_uniform_shape(self):
+    def _compute_uniform_shape(self, resolution=None):
         """Compute the height and width of a domain to have uniform resolution across dimensions."""
         g = Geod(ellps='WGS84')
 
@@ -861,14 +861,17 @@ def notnull(arr):
         az1, az2, width2 = g.inv(leftlons[-1], leftlats[-1], rightlons[-1], rightlats[-1])
         az1, az2, height = g.inv(middlelons[0], middlelats[0], middlelons[-1], middlelats[-1])
         width = min(width1, width2)
-        vresolution = height * 1.0 / self.lons.shape[0]
-        hresolution = width * 1.0 / self.lons.shape[1]
-        resolution = min(vresolution, hresolution)
+        if resolution is None:
+            vresolution = height * 1.0 / self.lons.shape[0]
+            hresolution = width * 1.0 / self.lons.shape[1]
+            resolution = (vresolution, hresolution)
+        if isinstance(resolution, (tuple, list)):
+            resolution = min(*resolution)
         width = int(width * 1.1 / resolution)
         height = int(height * 1.1 / resolution)
         return height, width
 
-    def compute_optimal_bb_area(self, proj_dict=None):
+    def compute_optimal_bb_area(self, proj_dict=None, resolution=None):
         """Compute the "best" bounding box area for this swath with `proj_dict`.
 
         By default, the projection is Oblique Mercator (`omerc` in proj.4), in
@@ -884,7 +887,7 @@ def compute_optimal_bb_area(self, proj_dict=None):
         projection = proj_dict.setdefault('proj', 'omerc')
         area_id = projection + '_otf'
         description = 'On-the-fly ' + projection + ' area'
-        height, width = self._compute_uniform_shape()
+        height, width = self._compute_uniform_shape(resolution)
         proj_dict = self.compute_bb_proj_params(proj_dict)
 
         area = DynamicAreaDefinition(area_id, description, proj_dict)
@@ -1024,7 +1027,7 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None):
         proj_info:
           complementing parameters to the projection info.
 
-        Resolution and shape parameters are ignored if the instance is created
+        Shape parameters are ignored if the instance is created
         with the `optimize_projection` flag set to True.
         """
         proj_dict = self._get_proj_dict()
@@ -1035,7 +1038,7 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None):
             projection = proj_dict
 
         if self.optimize_projection:
-            return lonslats.compute_optimal_bb_area(proj_dict)
+            return lonslats.compute_optimal_bb_area(proj_dict, resolution=resolution or self.resolution)
         if resolution is None:
             resolution = self.resolution
         if shape is None:

pyresample/test/test_files/areas.yaml

@@ -219,6 +219,14 @@ test_latlong:
     lower_left_xy: [-0.08115781021773638, 0.4038691889114878]
     upper_right_xy: [0.08115781021773638, 0.5427973973702365]
 
+omerc_bb_1000:
+  description: Oblique Mercator Bounding Box for Polar Overpasses
+  projection:
+    ellps: sphere
+    proj: omerc
+  optimize_projection: True
+  resolution: 1000
+
 test_dynamic_resolution:
   description: Dynamic with resolution specified in meters
   projection:

pyresample/test/test_geometry.py

@@ -1871,6 +1871,25 @@ def test_compute_optimal_bb(self):
         assert_np_dict_allclose(res.proj_dict, proj_dict)
         self.assertEqual(res.shape, (6, 3))
 
+    def test_compute_optimal_bb_with_resolution(self):
+        """Test computing the bb area while passing in the resolution."""
+        import xarray as xr
+        nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+                           [80.84000396728516, 60.74200439453125, 34.08500289916992],
+                           [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+        lats = xr.DataArray(nplats)
+        nplons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+                           [79.11000061035156, 7.284000396728516, -5.107000350952148],
+                           [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+        lons = xr.DataArray(nplons)
+
+        area = geometry.SwathDefinition(lons, lats)
+
+        res1000 = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'}, resolution=1000)
+        res10000 = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'}, resolution=10000)
+        assert res1000.shape[0] // 10 == res10000.shape[0]
+        assert res1000.shape[1] // 10 == res10000.shape[1]
+
     def test_aggregation(self):
         """Test aggregation on SwathDefinitions."""
         import dask.array as da
@@ -2343,6 +2362,49 @@ def test_freeze(self):
         assert result.width == 395
         assert result.height == 539
 
+    def test_freeze_when_area_is_optimized_and_has_a_resolution(self):
+        """Test freezing an optimized area with a resolution."""
+        import xarray as xr
+        nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+                           [80.84000396728516, 60.74200439453125, 34.08500289916992],
+                           [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+        lats = xr.DataArray(nplats)
+        nplons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+                           [79.11000061035156, 7.284000396728516, -5.107000350952148],
+                           [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+        lons = xr.DataArray(nplons)
+
+        swath = geometry.SwathDefinition(lons, lats)
+
+        area10km = geometry.DynamicAreaDefinition('test_area', 'A test area',
+                                                  {'ellps': 'WGS84', 'proj': 'omerc'},
+                                                  resolution=10000,
+                                                  optimize_projection=True)
+
+        result10km = area10km.freeze(swath)
+        assert result10km.shape == (679, 330)
+
+    def test_freeze_when_area_is_optimized_and_a_resolution_is_provided(self):
+        """Test freezing an optimized area when provided a resolution."""
+        import xarray as xr
+        nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+                           [80.84000396728516, 60.74200439453125, 34.08500289916992],
+                           [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+        lats = xr.DataArray(nplats)
+        nplons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+                           [79.11000061035156, 7.284000396728516, -5.107000350952148],
+                           [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+        lons = xr.DataArray(nplons)
+
+        swath = geometry.SwathDefinition(lons, lats)
+
+        area10km = geometry.DynamicAreaDefinition('test_area', 'A test area',
+                                                  {'ellps': 'WGS84', 'proj': 'omerc'},
+                                                  optimize_projection=True)
+
+        result10km = area10km.freeze(swath, 10000)
+        assert result10km.shape == (679, 330)
+
     @pytest.mark.parametrize(
         ('lats',),
         [
@@ -2399,7 +2461,7 @@ def test_freeze_with_bb(self):
                 [50, 51, 52, 53]]
         import xarray as xr
         sdef = geometry.SwathDefinition(xr.DataArray(lons), xr.DataArray(lats))
-        result = area.freeze(sdef, resolution=1000)
+        result = area.freeze(sdef)
         np.testing.assert_allclose(result.area_extent,
                                    [-335439.956533, 5502125.451125,
                                     191991.313351, 7737532.343683])

pyresample/test/test_utils.py

@@ -191,6 +191,13 @@ def test_dynamic_area_parser_yaml(self):
         self.assertTrue(hasattr(test_area, 'resolution'))
         np.testing.assert_allclose(test_area.resolution, (1.0, 1.0))
 
+    def test_dynamic_area_parser_passes_resolution(self):
+        """Test that the resolution from the file is passed to a dynamic area."""
+        from pyresample import parse_area_file
+        test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+        test_area = parse_area_file(test_area_file, 'omerc_bb_1000')[0]
+        assert test_area.resolution == (1000, 1000)
+
     def test_multiple_file_content(self):
         from pyresample import parse_area_file
         from pyresample.area_config import load_area_from_string