diff --git a/pyresample/area_config.py b/pyresample/area_config.py index 58b9eafbc..c4e6fefbf 100644 --- a/pyresample/area_config.py +++ b/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: diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 8cb919172..460796d77 100644 --- a/pyresample/geometry.py +++ b/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: diff --git a/pyresample/test/test_files/areas.yaml b/pyresample/test/test_files/areas.yaml index 4859d6e8c..40762ace2 100644 --- a/pyresample/test/test_files/areas.yaml +++ b/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: diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py index e5cec8a4b..caee8f8d7 100644 --- a/pyresample/test/test_geometry.py +++ b/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]) diff --git a/pyresample/test/test_utils.py b/pyresample/test/test_utils.py index 854792b79..d99f342f2 100644 --- a/pyresample/test/test_utils.py +++ b/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