Convert AreaDefinitions to odc geoboxes

continuous_integration/environment.yaml

@@ -27,3 +27,4 @@ dependencies:
   - pytest-lazy-fixture
   - importlib-metadata
   - sphinx-reredirects
+  - odc-geo

pyresample/geometry.py

@@ -57,6 +57,11 @@
 except ImportError:
     da = None
 
+try:
+    import odc.geo as odc_geo
+except ModuleNotFoundError:
+    odc_geo = None
+
 from pyproj import CRS
 from pyproj.enums import TransformDirection
 
@@ -354,7 +359,7 @@
         new_dim1_sides = []
         new_dim2_sides = []
         for dim1_side, dim2_side in zip(dim1_sides, dim2_sides):
             # FIXME: ~(~np.isfinite(dim1_side) | ~np.isfinite(dim1_side))
             is_valid_mask = ~(np.isnan(dim1_side) | np.isnan(dim2_side))
             if not is_valid_mask.any():
                 raise ValueError("Can't compute boundary coordinates. At least one side is completely invalid.")
@@ -363,7 +368,7 @@
         return new_dim1_sides, new_dim2_sides
 
     def _get_bbox_slices(self, vertices_per_side):
         # FIXME: This currently replicate values if heigh/width < row_num/col_num !
         height, width = self.shape
         if vertices_per_side is None:
             row_num = height
@@ -630,7 +635,7 @@
                 (i.e. full disc geostationary area, Robinson projection, polar projections, ...)
                 by default only 50 points are selected.
         """
         # FIXME: Add logic for out-of-earth disk projections
         if self.is_geostationary:
             return self._get_geostationary_boundary_sides(vertices_per_side=vertices_per_side, coordinates="geographic")
         sides_lons, sides_lats = self._get_sides(coord_fun=self.get_lonlats, vertices_per_side=vertices_per_side)
@@ -1998,6 +2003,18 @@
             return "EPSG:{}".format(self.crs.to_epsg())
         return self.crs.to_proj4()
 
+    def to_odc_geobox(self):
+        """Convert AreaDefinition to ODC GeoBox.
+
+        See: https://odc-geo.readthedocs.io/en/latest/
+        """
+        if odc_geo is None:
+            raise ModuleNotFoundError("Please install 'odc-geo' to use this method.")
+
+        return odc_geo.geobox.GeoBox.from_bbox(bbox=self.area_extent, crs=self.crs,
+                                               resolution=odc_geo.Resolution(x=self.pixel_size_x, y=-self.pixel_size_y),
+                                               tight=True)
+
     def create_areas_def(self):
         """Generate YAML formatted representation of this area.
 
@@ -2723,7 +2740,7 @@
             Each list element is a numpy array representing a specific side of the geometry.
             The order of the sides are [top", "right", "bottom", "left"]
         """
         # FIXME: Add logic for out-of-earth-disk
         if self.is_geostationary:
             return self._get_geostationary_boundary_sides(vertices_per_side=vertices_per_side,
                                                           coordinates="projection")

pyresample/test/test_geometry/test_area.py

@@ -404,6 +404,45 @@ def test_cartopy_crs_latlon_bounds(self, create_test_area):
         latlong_crs = area_def.to_cartopy_crs()
         np.testing.assert_allclose(latlong_crs.bounds, [-180, 180, -90, 90])
 
+    def test_to_odc_geobox_odc_missing(self, monkeypatch, stere_area):
+        """Test odc-geo not installed."""
+        area = stere_area
+
+        with monkeypatch.context() as m:
+            m.setattr(pyresample.geometry, "odc_geo", None)
+
+            with pytest.raises(ModuleNotFoundError):
+                area.to_odc_geobox()
+
+    def test_to_odc_geobox(self, stere_area, create_test_area):
+        """Test conversion from area definition to odc GeoBox."""
+        from odc.geo.geobox import GeoBox
+
+        europe = stere_area
+        seviri = create_test_area(
+            {
+                'proj': 'geos',
+                'lon_0': 0.0,
+                'a': 6378169.00,
+                'b': 6356583.80,
+                'h': 35785831.00,
+                'units': 'm'},
+            123,
+            123,
+            (-5500000, -5500000, 5500000, 5500000))
+
+        for area_def in [europe, seviri]:
+            geobox = area_def.to_odc_geobox()
+
+            assert isinstance(geobox, GeoBox)
+
+            # Affine coefficiants
+            af = geobox.affine
+            assert af.a == area_def.pixel_size_x
+            assert af.e == -area_def.pixel_size_y
+            assert af.xoff == area_def.area_extent[0]
+            assert af.yoff == area_def.area_extent[3]
+
     def test_dump(self, create_test_area):
         """Test exporting area defs."""
         from io import StringIO

setup.py

@@ -32,16 +32,22 @@
     requirements.append('importlib_metadata')
 
 test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr',
-                 'pytest-lazy-fixtures', 'shapely']
+                 'pytest-lazy-fixtures', 'shapely', 'odc-geo']
 extras_require = {'numexpr': ['numexpr'],
                   'quicklook': ['matplotlib', 'cartopy>=0.20.0', 'pillow'],
                   'rasterio': ['rasterio'],
                   'dask': ['dask>=0.16.1'],
                   'cf': ['xarray'],
                   'gradient_search': ['shapely'],
                   'xarray_bilinear': ['xarray', 'dask', 'zarr'],
+                  'odc-geo': ['odc-geo'],
                   'tests': test_requires}
 
+all_extras = []
+for extra_deps in extras_require.values():
+    all_extras.extend(extra_deps)
+extras_require['all'] = list(set(all_extras))
+
 if sys.platform.startswith("win"):
     extra_compile_args = []
 else: