Merge pull request #306 from gerritholl/area-join-utility

Add a function to join/enclose areas.

pyresample/geometry.py

@@ -2,12 +2,7 @@
 # -*- coding: utf-8 -*-
 # pyresample, Resampling of remote sensing image data in python
 #
-# Copyright (C) 2010-2016
-#
-# Authors:
-# Esben S. Nielsen
-# Thomas Lavergne
-# Adam Dybbroe
+# Copyright (C) 2010-2020 Pyresample developers
 #
 # This program is free software: you can redistribute it and/or modify it under
 # the terms of the GNU Lesser General Public License as published by the Free
@@ -2439,3 +2434,43 @@ def _dict_representer(dumper, data):
 
  OrderedDumper.add_representer(OrderedDict, _dict_representer)
  return yaml.dump(data, stream, OrderedDumper, **kwds)
+
+
+def enclose_areas(*areas, area_id="joint-area"):
+ """Return the smallest areadefinition enclosing one or more others.
+
+ From one or more AreaDefinition objects (most usefully at least
+ two), which shall differ only in extent, calculate the smallest
+ AreaDefinition that encloses all. Touches only the ``area_extent``;
+ projection and units must be identical in all input areas and will be
+ unchanged in the resulting area. When the input areas :math:`i=1..n`
+ have extent :math:`(a_i, b_i, c_i, d_i)`, the resulting area will have
+ extent :math:`(\\min_i{a_i}, \\min_i{b_i}, \\max_i{c_i}, \\max_i{d_i})`.
+
+ Args:
+ *areas (AreaDefinition): AreaDefinition objects to enclose.
+ area_id (Optional[str]): Name of joint area, defaults to "joint-area".
+ """
+
+ first = None
+ if not areas:
+ raise TypeError("Must pass at least one area, found zero.")
+ for area in areas:
+ if first is None:
+ first = area
+ largest_extent = list(area.area_extent)
+ else:
+ if not area.crs == first.crs:
+ raise ValueError("Inconsistent projections between areas")
+ if not np.isclose(area.resolution, first.resolution).all():
+ raise ValueError("Inconsistent resolution between areas")
+ largest_extent[0] = min(largest_extent[0], area.area_extent[0])
+ largest_extent[1] = min(largest_extent[1], area.area_extent[1])
+ largest_extent[2] = max(largest_extent[2], area.area_extent[2])
+ largest_extent[3] = max(largest_extent[3], area.area_extent[3])
+
+ return create_area_def(
+ area_id=area_id,
+ projection=first.crs,
+ area_extent=largest_extent,
+ resolution=first.resolution)

pyresample/test/test_geometry.py

@@ -2,12 +2,7 @@
 # -*- coding: utf-8 -*-
 # pyresample, Resampling of remote sensing image data in python
 #
-# Copyright (C) 2010-2016
-#
-# Authors:
-# Esben S. Nielsen
-# Thomas Lavergne
-# Adam Dybbroe
+# Copyright (C) 2010-2020 Pyresample developers
 #
 # This program is free software: you can redistribute it and/or modify it under
 # the terms of the GNU Lesser General Public License as published by the Free
@@ -24,6 +19,7 @@
 """Test the geometry objects."""
 import random
 import sys
+import pytest
 
 import numpy as np
 
@@ -36,7 +32,6 @@
 from unittest.mock import MagicMock, patch
 import unittest
 import pyproj
-import pytest
 
 try:
  from pyproj import CRS
@@ -2506,3 +2501,56 @@ def test_aggregate(self):
  np.testing.assert_allclose(res.area_extent, area.area_extent)
  self.assertEqual(res.shape[0], area.shape[0] / 2)
  self.assertEqual(res.shape[1], area.shape[1] / 4)
+
+
+def test_enclose_areas():
+ from pyresample.geometry import (enclose_areas, create_area_def)
+ proj_dict = {'proj': 'geos', 'sweep': 'x', 'lon_0': 0, 'h': 35786023,
+ 'x_0': 0, 'y_0': 0, 'ellps': 'GRS80', 'units': 'm',
+ 'no_defs': None, 'type': 'crs'}
+ proj_dict_alt = {'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0,
+ 'units': 'm'}
+
+ ar1 = create_area_def(
+ "test-area",
+ projection=proj_dict,
+ units="m",
+ area_extent=[0, 20, 100, 120],
+ shape=(10, 10))
+
+ ar2 = create_area_def(
+ "test-area",
+ projection=proj_dict,
+ units="m",
+ area_extent=[20, 40, 120, 140],
+ shape=(10, 10))
+
+ ar3 = create_area_def(
+ "test-area",
+ projection=proj_dict,
+ units="m",
+ area_extent=[20, 0, 120, 100],
+ shape=(10, 10))
+
+ ar4 = create_area_def(
+ "test-area",
+ projection=proj_dict_alt,
+ units="m",
+ area_extent=[20, 0, 120, 100],
+ shape=(10, 10))
+
+ ar5 = create_area_def(
+ "test-area",
+ projection=proj_dict,
+ units="m",
+ area_extent=[-50, -50, 50, 50],
+ shape=(100, 100))
+
+ ar_joined = enclose_areas(ar1, ar2, ar3)
+ np.testing.assert_allclose(ar_joined.area_extent, [0, 0, 120, 140])
+ with pytest.raises(ValueError):
+ enclose_areas(ar3, ar4)
+ with pytest.raises(ValueError):
+ enclose_areas(ar3, ar5)
+ with pytest.raises(TypeError):
+ enclose_areas()