.travis.yml

@@ -4,12 +4,13 @@ env:
   - PYTHON_VERSION=$PYTHON_VERSION
   - NUMPY_VERSION=stable
   - MAIN_CMD='python setup.py'
-  - CONDA_DEPENDENCIES='xarray dask toolz Cython pykdtree sphinx cartopy pillow matplotlib
+  - CONDA_DEPENDENCIES='xarray dask toolz Cython pykdtree sphinx cartopy rasterio pillow matplotlib
     pyyaml pyproj coveralls configobj coverage codecov'
   - SETUP_XVFB=False
   - EVENT_TYPE='push pull_request'
   - SETUP_CMD='test'
   - CONDA_CHANNELS='conda-forge'
+  - CONDA_CHANNEL_PRIORITY='True'
 matrix:
   include:
   - env: PYTHON_VERSION=2.7

CHANGELOG.md

+## Version 1.10.3 (2018/11/23)
+
+### Issues Closed
+
+* [Issue 92](https://github.com/pytroll/pyresample/issues/92) - Add utility function for converting geotiffs to area definitions ([PR 143](https://github.com/pytroll/pyresample/pull/143))
+
+In this release 1 issue was closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 147](https://github.com/pytroll/pyresample/pull/147) - Fix dtype preservation for kdtree resampling
+* [PR 144](https://github.com/pytroll/pyresample/pull/144) - Non-contiguous area definitions are now not concatenable ([491](https://github.com/pytroll/satpy/issues/491))
+
+#### Features added
+
+* [PR 143](https://github.com/pytroll/pyresample/pull/143) - get_area_def_from_raster ([92](https://github.com/pytroll/pyresample/issues/92))
+* [PR 142](https://github.com/pytroll/pyresample/pull/142) - Add converter from def to yaml
+
+In this release 4 pull requests were closed.
+
+
 ## Version 1.10.2 (2018/10/01)
 
 ### Issues Closed

appveyor.yml

@@ -3,8 +3,9 @@ environment:
     PYTHON: "C:\\conda"
     MINICONDA_VERSION: "latest"
     CMD_IN_ENV: "cmd /E:ON /V:ON /C .\\ci-helpers\\appveyor\\windows_sdk.cmd"
-    CONDA_DEPENDENCIES: "xarray dask toolz Cython pykdtree sphinx cartopy pillow matplotlib pyyaml pyproj coveralls configobj coverage"
+    CONDA_DEPENDENCIES: "xarray dask toolz Cython pykdtree sphinx cartopy rasterio pillow matplotlib pyyaml pyproj coveralls configobj coverage"
     CONDA_CHANNELS: "conda-forge"
+    CONDA_CHANNEL_PRIORITY: "True"
 
   matrix:
     - PYTHON: "C:\\Python27_32"

debian/changelog

+pyresample (1.10.3-1) unstable; urgency=medium
+
+  * New upstream release.
+  * debian/control
+    - add dependency form rasterio
+  * debian/patches
+    - refresh all patches
+
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Sat, 01 Dec 2018 20:10:19 +0000
+
 pyresample (1.10.2-1) unstable; urgency=medium
 
   * New upstream release

debian/control

@@ -33,6 +33,8 @@ Build-Depends: debhelper (>= 11.0.0),
                python3-six,
                python-scipy,
                python3-scipy,
+	       python-rasterio,
+	       python3-rasterio,
                cython,
                cython3
 Standards-Version: 4.2.1
@@ -51,6 +53,7 @@ Depends: ${misc:Depends},
          python-yaml,
          python-pykdtree (>= 1.3.1),
          python-six,
+	 python-rasterio,
          python-pyresample-test
 Recommends: python-numexpr,
             python-pil,
@@ -80,6 +83,7 @@ Depends: ${misc:Depends},
          python3-yaml,
          python3-pykdtree (>= 1.3.1),
          python3-six,
+	 python3-rasterio,
          python-pyresample-test
 Recommends: python3-numexpr,
             python3-pil,

debian/patches/0001-fix-doc-build.patch

@@ -7,55 +7,55 @@ Subject: fix doc build
  1 file changed, 9 insertions(+), 9 deletions(-)
 
 diff --git a/docs/source/API.rst b/docs/source/API.rst
-index 1731916..d73c6f3 100644
+index 5282938..8d20c17 100644
 --- a/docs/source/API.rst
 +++ b/docs/source/API.rst
-@@ -3,47 +3,47 @@ pyresample API
+@@ -3,45 +3,45 @@ pyresample API
  
  pyresample.geometry
- ---------------------------------
+ -------------------
 -.. automodule:: geometry
 +.. automodule:: pyresample.geometry
      :members:
  
  pyresample.image
- ---------------------------------
+ ----------------
 -.. automodule:: image
 +.. automodule:: pyresample.image
      :members:
  
  pyresample.grid
- ---------------------------------
+ ---------------
 -.. automodule:: grid
 +.. automodule:: pyresample.grid
      :members:
  
  pyresample.kd_tree
- ---------------------------------
+ ------------------
 -.. automodule:: kd_tree
 +.. automodule:: pyresample.kd_tree
      :members:
  
  pyresample.bilinear
- ---------------------------------
+ -------------------
 -.. automodule:: bilinear
 +.. automodule:: pyresample.bilinear
      :members:
  
  pyresample.utils
- ---------------------------------
+ ----------------
 -.. automodule:: utils
 +.. automodule:: pyresample.utils
      :members:
  
  pyresample.data_reduce
- ---------------------------------
+ ----------------------
 -.. automodule:: data_reduce
 +.. automodule:: pyresample.data_reduce
      :members:
  
  pyresample.plot
- ---------------------------------
+ ---------------
 -.. automodule:: plot
 +.. automodule:: pyresample.plot
      :members:
@@ -65,5 +65,3 @@ index 1731916..d73c6f3 100644
 -.. automodule:: ewa
 +.. automodule:: pyresample.ewa
      :members:
- 
- 	

debian/patches/0002-fix-proj4-initialization.patch

@@ -8,10 +8,10 @@ Subject: fix proj4 initialization
  2 files changed, 5 insertions(+), 5 deletions(-)
 
 diff --git a/pyresample/test/test_files/areas.cfg b/pyresample/test/test_files/areas.cfg
-index 3c6ef3c..5346cb4 100644
+index 063264d..620e665 100644
 --- a/pyresample/test/test_files/areas.cfg
 +++ b/pyresample/test/test_files/areas.cfg
-@@ -19,7 +19,7 @@ REGION: ease_nh {
+@@ -28,7 +28,7 @@ REGION: ease_nh {
  REGION: pc_world {
    NAME:    Plate Carree world map
    PCS_ID:  pc_world
@@ -21,10 +21,10 @@ index 3c6ef3c..5346cb4 100644
    YSIZE: 480
    AREA_EXTENT:  (-20037508.342789244, -10018754.171394622, 20037508.342789244, 10018754.171394622)
 diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
-index 8c6aa9a..11ec7cd 100644
+index 19520b9..2f52e8c 100644
 --- a/pyresample/test/test_geometry.py
 +++ b/pyresample/test/test_geometry.py
-@@ -401,7 +401,7 @@ class Test(unittest.TestCase):
+@@ -428,7 +428,7 @@ class Test(unittest.TestCase):
          swath_def = geometry.SwathDefinition(lons, lats)
          filter_area = geometry.AreaDefinition('test', 'test', 'test',
                                                {'proj': 'eqc', 'lon_0': 0.0,
@@ -33,7 +33,7 @@ index 8c6aa9a..11ec7cd 100644
                                                8, 8,
                                                (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
          filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0],
-@@ -426,7 +426,7 @@ class Test(unittest.TestCase):
+@@ -453,7 +453,7 @@ class Test(unittest.TestCase):
          data = np.array([1, 2, 3, 4])
          filter_area = geometry.AreaDefinition('test', 'test', 'test',
                                                {'proj': 'eqc', 'lon_0': 0.0,
@@ -42,7 +42,7 @@ index 8c6aa9a..11ec7cd 100644
                                                8, 8,
                                                (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
          filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0],
-@@ -461,7 +461,7 @@ class Test(unittest.TestCase):
+@@ -488,7 +488,7 @@ class Test(unittest.TestCase):
          data = np.dstack((data1, data2, data3))
          filter_area = geometry.AreaDefinition('test', 'test', 'test',
                                                {'proj': 'eqc', 'lon_0': 0.0,
@@ -51,7 +51,7 @@ index 8c6aa9a..11ec7cd 100644
                                                8, 8,
                                                (-20037508.34, -10018754.17, 20037508.34, 10018754.17))
          filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0],
-@@ -530,7 +530,7 @@ class Test(unittest.TestCase):
+@@ -557,7 +557,7 @@ class Test(unittest.TestCase):
  
      def test_latlong_area(self):
          area_def = geometry.AreaDefinition('', '', '',

debian/patches/0003-Skip-TestXArrayResamplerNN-if-dask-is-not-available.patch

@@ -7,7 +7,7 @@ Subject: Skip TestXArrayResamplerNN if dask is not available
  1 file changed, 6 insertions(+)
 
 diff --git a/pyresample/test/test_kd_tree.py b/pyresample/test/test_kd_tree.py
-index dcb3943..fbf91ef 100644
+index ccaea19..65007cc 100644
 --- a/pyresample/test/test_kd_tree.py
 +++ b/pyresample/test/test_kd_tree.py
 @@ -14,6 +14,11 @@ if sys.version_info < (2, 7):

docs/source/API.rst

 pyresample API
-======================
+==============
 
 pyresample.geometry
----------------------------------
+-------------------
 .. automodule:: geometry
     :members:
 
 pyresample.image
----------------------------------
+----------------
 .. automodule:: image
     :members:
 
 pyresample.grid
----------------------------------
+---------------
 .. automodule:: grid
     :members:
 
 pyresample.kd_tree
----------------------------------
+------------------
 .. automodule:: kd_tree
     :members:
 
 pyresample.bilinear
----------------------------------
+-------------------
 .. automodule:: bilinear
     :members:
 
 pyresample.utils
----------------------------------
+----------------
 .. automodule:: utils
     :members:
 
 pyresample.data_reduce
----------------------------------
+----------------------
 .. automodule:: data_reduce
     :members:
 
 pyresample.plot
----------------------------------
+---------------
 .. automodule:: plot
     :members:
 
@@ -45,5 +45,3 @@ pyresample.ewa
 --------------
 .. automodule:: ewa
     :members:
-
-	

pyresample/geometry.py

@@ -920,21 +920,18 @@ class AreaDefinition(BaseDefinition):
         return crs
 
     def create_areas_def(self):
-        to_dump = OrderedDict()
-        res = OrderedDict()
-        to_dump[self.area_id] = res
-
-        res['description'] = self.name
-        res['shape'] = OrderedDict([('height', self.y_size),
-                                    ('width', self.x_size)])
-        res['area_extent'] = OrderedDict([('lower_left_xy',
-                                           list(self.area_extent[:2])),
-                                          ('upper_right_xy',
-                                           list(self.area_extent[2:])),
-                                          ('units', 'm')
-                                          ])
-
-        return ordered_dump(to_dump)
+
+        res = OrderedDict(description=self.name,
+                          projection=OrderedDict(self.proj_dict),
+                          shape=OrderedDict([('height', self.y_size), ('width', self.x_size)]))
+        units = res['projection'].pop('units', None)
+        extent = OrderedDict([('lower_left_xy', list(self.area_extent[:2])),
+                              ('upper_right_xy', list(self.area_extent[2:]))])
+        if units is not None:
+            extent['units'] = units
+        res['area_extent'] = extent
+
+        return ordered_dump(OrderedDict([(self.area_id, res)]))
 
     def create_areas_def_legacy(self):
         proj_dict = self.proj_dict
@@ -1452,13 +1449,17 @@ def get_geostationary_bounding_box(geos_area, nb_points=50):
 
 def combine_area_extents_vertical(area1, area2):
     """Combine the area extents of areas 1 and 2."""
-    if (area1.area_extent[0] == area2.area_extent[0] and
-            area1.area_extent[2] == area2.area_extent[2]):
+    if (area1.area_extent[0] == area2.area_extent[0]
+            and area1.area_extent[2] == area2.area_extent[2]):
         current_extent = list(area1.area_extent)
         if np.isclose(area1.area_extent[1], area2.area_extent[3]):
             current_extent[1] = area2.area_extent[1]
         elif np.isclose(area1.area_extent[3], area2.area_extent[1]):
             current_extent[3] = area2.area_extent[3]
+        else:
+            raise IncompatibleAreas(
+                "Can't concatenate non-contiguous area definitions: "
+                "{0} and {1}".format(area1, area2))
     else:
         raise IncompatibleAreas(
             "Can't concatenate area definitions with "

pyresample/kd_tree.py

@@ -1077,12 +1077,16 @@ class XArrayResamplerNN(object):
         Args:
             data (dask.array.Array): Source data pixels to sample
             fill_value (float): Output fill value when no source data is
-                                near the target pixel. If the input data
-                                is a integer array then the minimum value
-                                for that integer type is used. Otherwise,
-                                NaN is used and can be detected in the result
-                                with ``res.isnull()``.
-
+                near the target pixel. When omitted, if the input data is an
+                integer array then the maximum value for that integer type is
+                used, but otherwise, NaN is used and can be detected in the
+                result with ``res.isnull()``.
+
+        Returns:
+            dask.array.Array: The resampled array. The dtype of the array will
+                be the same as the input data. Pixels with no matching data from
+                the input array will be filled (see the `fill_value` parameter
+                description above).
         """
         if fill_value is not None and np.isnan(fill_value) and \
                 np.issubdtype(data.dtype, np.integer):
@@ -1200,10 +1204,7 @@ class XArrayResamplerNN(object):
                       dtype=new_data.dtype, concatenate=True)
         res = DataArray(res, dims=dst_dims, coords=coords,
                         attrs=data.attrs.copy())
-        res.attrs['_FillValue'] = fill_value
-        # if fill_value isn't NaN then we have to tell xarray what null is
-        if not np.isnan(fill_value):
-            res = res.where(res != fill_value)
+
         return res
 
 

pyresample/test/test_files/areas.cfg

@@ -16,6 +16,15 @@ REGION: ease_nh {
         AREA_EXTENT:    (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
 };
 
+#REGION: commented {
+#         NAME:           Arctic EASE grid
+#         PCS_ID:         ease_nh
+#         PCS_DEF:        proj=laea, lat_0=90, lon_0=0, a=6371228.0, units=m
+#         XSIZE:          425
+#         YSIZE:          425
+#         AREA_EXTENT:    (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
+# };
+
 REGION: pc_world {
   NAME:    Plate Carree world map
   PCS_ID:  pc_world

pyresample/test/test_geometry.py

@@ -89,6 +89,33 @@ class Test(unittest.TestCase):
                           area_def.area_extent[1],
                           area_def.area_extent[3]))
 
+    def test_create_areas_def(self):
+        area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)',
+                                           'areaD',
+                                           {'a': '6378144.0',
+                                            'b': '6356759.0',
+                                            'lat_0': '50.00',
+                                            'lat_ts': '50.00',
+                                            'lon_0': '8.00',
+                                            'proj': 'stere'},
+                                           800,
+                                           800,
+                                           [-1370912.72,
+                                            -909968.64000000001,
+                                            1029087.28,
+                                            1490031.3600000001])
+        import yaml
+        res = yaml.load(area_def.create_areas_def())
+        expected = yaml.load(('areaD:\n  description: Europe (3km, HRV, VTC)\n'
+                              '  projection:\n    a: 6378144.0\n    b: 6356759.0\n'
+                              '    lat_0: 50.0\n    lat_ts: 50.0\n    lon_0: 8.0\n'
+                              '    proj: stere\n  shape:\n    height: 800\n'
+                              '    width: 800\n  area_extent:\n'
+                              '    lower_left_xy: [-1370912.72, -909968.64]\n'
+                              '    upper_right_xy: [1029087.28, 1490031.36]\n'))
+
+        self.assertDictEqual(res, expected)
+
     def test_base_type(self):
         lons1 = np.arange(-135., +135, 50.)
         lats = np.ones_like(lons1) * 70.
@@ -866,14 +893,14 @@ class Test(unittest.TestCase):
                                        proj_dict, 10, 10,
                                        [-1370912.72, -909968.64, 1029087.28,
                                         1490031.36])
-        self.assertEquals(area.proj_str,
+        self.assertEqual(area.proj_str,
                          '+a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0 +proj=stere')
         proj_dict['no_rot'] = ''
         area = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
                                        proj_dict, 10, 10,
                                        [-1370912.72, -909968.64, 1029087.28,
                                         1490031.36])
-        self.assertEquals(area.proj_str,
+        self.assertEqual(area.proj_str,
                          '+a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0 +no_rot +proj=stere')
 
 
@@ -1239,6 +1266,14 @@ class TestStackedAreaDefinition(unittest.TestCase):
         res = combine_area_extents_vertical(area1, area2)
         self.assertListEqual(res, [1, 2, 3, 6])
 
+        # Non contiguous area extends shouldn't be combinable
+        area1 = MagicMock()
+        area1.area_extent = (1, 2, 3, 4)
+        area2 = MagicMock()
+        area2.area_extent = (1, 5, 3, 7)
+        self.assertRaises(IncompatibleAreas,
+                          combine_area_extents_vertical, area1, area2)
+
     def test_append_area_defs_fail(self):
         """Fail appending areas."""
         area1 = MagicMock()

pyresample/test/test_kd_tree.py

@@ -797,6 +797,34 @@ class TestXArrayResamplerNN(unittest.TestCase):
         ])
         np.testing.assert_allclose(actual, expected)
 
+    def test_nearest_type_preserve(self):
+        """Test 1D swath definition to 2D grid definition; 1 neighbor."""
+        from pyresample.kd_tree import XArrayResamplerNN
+        import xarray as xr
+        import dask.array as da
+        resampler = XArrayResamplerNN(self.tswath_1d, self.tgrid,
+                                      radius_of_influence=100000,
+                                      neighbours=1)
+        data = self.tdata_1d
+        data = xr.DataArray(da.from_array(np.array([1, 2, 3]),
+                                          chunks=5),
+                            dims=('my_dim1',))
+        ninfo = resampler.get_neighbour_info()
+        for val in ninfo[:3]:
+            # vii, ia, voi
+            self.assertIsInstance(val, da.Array)
+        res = resampler.get_sample_from_neighbour_info(data, fill_value=255)
+        self.assertIsInstance(res, xr.DataArray)
+        self.assertIsInstance(res.data, da.Array)
+        actual = res.values
+        expected = np.array([
+            [1, 2, 2],
+            [1, 2, 2],
+            [1, 255, 2],
+            [1, 2, 2],
+        ])
+        np.testing.assert_equal(actual, expected)
+
     def test_nearest_swath_2d_mask_to_area_1n(self):
         """Test 2D swath definition to 2D area definition; 1 neighbor."""
         from pyresample.kd_tree import XArrayResamplerNN

pyresample/test/test_utils.py

@@ -2,6 +2,7 @@ import os
 import unittest
 
 import numpy as np
+import uuid
 
 from pyresample.test.utils import create_test_longitude, create_test_latitude
 
@@ -11,6 +12,16 @@ def tmp(f):
     return f
 
 
+def tmptiff(width=100, height=100, transform=None, crs=None, dtype=np.uint8):
+    import rasterio
+    array = np.ones((width, height)).astype(dtype)
+    fname = '/vsimem/%s' % uuid.uuid4()
+    with rasterio.open(fname, 'w', driver='GTiff', count=1, transform=transform,
+                       width=width, height=height, crs=crs, dtype=dtype) as dst:
+        dst.write(array, 1)
+    return fname
+
+
 class TestLegacyAreaParser(unittest.TestCase):
     def test_area_parser_legacy(self):
         """Test legacy area parser."""
@@ -26,7 +37,7 @@ Projection: {'a': '6371228.0', 'lat_0': '90.0', 'lon_0': '0.0', 'proj': 'laea',
 Number of columns: 425
 Number of rows: 425
 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
-        self.assertEquals(ease_nh.__str__(), nh_str)
+        self.assertEqual(ease_nh.__str__(), nh_str)
         self.assertIsInstance(ease_nh.proj_dict['lat_0'], float)
 
         sh_str = """Area ID: ease_sh
@@ -36,7 +47,7 @@ Projection: {'a': '6371228.0', 'lat_0': '-90.0', 'lon_0': '0.0', 'proj': 'laea',
 Number of columns: 425
 Number of rows: 425
 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
-        self.assertEquals(ease_sh.__str__(), sh_str)
+        self.assertEqual(ease_sh.__str__(), sh_str)
         self.assertIsInstance(ease_sh.proj_dict['lat_0'], float)
 
     def test_load_area(self):
@@ -51,7 +62,7 @@ Projection: {'a': '6371228.0', 'lat_0': '90.0', 'lon_0': '0.0', 'proj': 'laea',
 Number of columns: 425
 Number of rows: 425
 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
-        self.assertEquals(nh_str, ease_nh.__str__())
+        self.assertEqual(nh_str, ease_nh.__str__())
 
     def test_not_found_exception(self):
         from pyresample import utils
@@ -60,6 +71,12 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
                               os.path.dirname(__file__), 'test_files', 'areas.cfg'),
                           'no_area')
 
+    def test_commented(self):
+        from pyresample import utils
+        areas = utils.parse_area_file(os.path.join(os.path.dirname(__file__),
+                                                   'test_files', 'areas.cfg'))
+        self.assertNotIn('commented', [area.name for area in areas])
+
 
 class TestYAMLAreaParser(unittest.TestCase):
     def test_area_parser_yaml(self):
@@ -76,7 +93,7 @@ Projection: {'a': '6371228.0', 'lat_0': '90.0', 'lon_0': '0.0', 'proj': 'laea',
 Number of columns: 425
 Number of rows: 425
 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
-        self.assertEquals(ease_nh.__str__(), nh_str)
+        self.assertEqual(ease_nh.__str__(), nh_str)
 
         sh_str = """Area ID: ease_sh
 Description: Antarctic EASE grid
@@ -84,7 +101,7 @@ Projection: {'a': '6371228.0', 'lat_0': '-90.0', 'lon_0': '0.0', 'proj': 'laea',
 Number of columns: 425
 Number of rows: 425
 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
-        self.assertEquals(ease_sh.__str__(), sh_str)
+        self.assertEqual(ease_sh.__str__(), sh_str)
 
     def test_multiple_file_content(self):
         from pyresample import utils
@@ -121,7 +138,7 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
     units: m
 """]
         results = utils.parse_area_file(area_list)
-        self.assertEquals(len(results), 2)
+        self.assertEqual(len(results), 2)
         self.assertIn(results[0].area_id, ('ease_sh', 'ease_sh2'))
         self.assertIn(results[1].area_id, ('ease_sh', 'ease_sh2'))
 
@@ -251,6 +268,79 @@ class TestMisc(unittest.TestCase):
         self.assertIsInstance(proj_dict['lon_0'], float)
         self.assertIsInstance(proj_dict2['lon_0'], float)
 
+    def test_def2yaml_converter(self):
+        from pyresample import utils
+        import tempfile
+        def_file = os.path.join(os.path.dirname(__file__), 'test_files',
+                                'areas.cfg')
+        filehandle, yaml_file = tempfile.mkstemp()
+        os.close(filehandle)
+        try:
+            utils.convert_def_to_yaml(def_file, yaml_file)
+            areas_new = set(utils.parse_area_file(yaml_file))
+            areas_old = set(utils.parse_area_file(def_file))
+
+            self.assertEqual(areas_new, areas_old)
+        finally:
+            os.remove(yaml_file)
+
+    def test_get_area_def_from_raster(self):
+        from rasterio.crs import CRS
+        from affine import Affine
+        from pyresample import utils
+        x_size = 791
+        y_size = 718
+        transform = Affine(300.0379266750948, 0.0, 101985.0,
+                           0.0, -300.041782729805, 2826915.0)
+        crs = CRS(init='epsg:3857')
+        source = tmptiff(x_size, y_size, transform, crs=crs)
+        area_id = 'area_id'
+        proj_id = 'proj_id'
+        name = 'name'
+        area_def = utils.get_area_def_from_raster(source, area_id=area_id, name=name, proj_id=proj_id)
+        self.assertEqual(area_def.area_id, area_id)
+        self.assertEqual(area_def.proj_id, proj_id)
+        self.assertEqual(area_def.name, name)
+        self.assertEqual(area_def.x_size, x_size)
+        self.assertEqual(area_def.y_size, y_size)
+        self.assertDictEqual(crs.to_dict(), area_def.proj_dict)
+        self.assertTupleEqual(area_def.area_extent, (transform.c, transform.f + transform.e * y_size,
+                                                     transform.c + transform.a * x_size, transform.f))
+
+    def test_get_area_def_from_raster_extracts_proj_id(self):
+        from rasterio.crs import CRS
+        from pyresample import utils
+        crs = CRS(init='epsg:3857')
+        source = tmptiff(crs=crs)
+        area_def = utils.get_area_def_from_raster(source)
+        self.assertEqual(area_def.proj_id, 'WGS 84 / Pseudo-Mercator')
+
+    def test_get_area_def_from_raster_rotated_value_err(self):
+        from pyresample import utils
+        from affine import Affine
+        transform = Affine(300.0379266750948, 0.1, 101985.0,
+                           0.0, -300.041782729805, 2826915.0)
+        source = tmptiff(transform=transform)
+        self.assertRaises(ValueError, utils.get_area_def_from_raster, source)
+
+    def test_get_area_def_from_raster_non_georef_value_err(self):
+        from pyresample import utils
+        from affine import Affine
+        transform = Affine(300.0379266750948, 0.0, 101985.0,
+                           0.0, -300.041782729805, 2826915.0)
+        source = tmptiff(transform=transform)
+        self.assertRaises(ValueError, utils.get_area_def_from_raster, source)
+
+    def test_get_area_def_from_raster_non_georef_respects_proj_dict(self):
+        from pyresample import utils
+        from affine import Affine
+        transform = Affine(300.0379266750948, 0.0, 101985.0,
+                           0.0, -300.041782729805, 2826915.0)
+        source = tmptiff(transform=transform)
+        proj_dict = {'init': 'epsg:3857'}
+        area_def = utils.get_area_def_from_raster(source, proj_dict=proj_dict)
+        self.assertDictEqual(area_def.proj_dict, proj_dict)
+
 
 def suite():
     """The test suite.
@@ -263,3 +353,7 @@ def suite():
     mysuite.addTest(loader.loadTestsFromTestCase(TestMisc))
 
     return mysuite
+
+
+if __name__ == '__main__':
+    unittest.main()

pyresample/utils.py

@@ -28,7 +28,7 @@ import numpy as np
 import six
 import yaml
 from configobj import ConfigObj
-from collections import Mapping
+from collections import Mapping, OrderedDict
 
 
 class AreaNotFound(KeyError):
@@ -149,6 +149,10 @@ def _parse_yaml_area_file(area_file_name, *regions):
                            params['area_extent']['upper_right_xy'])
         except KeyError:
             area_extent = None
+        try:
+            projection['units'] = params['area_extent']['units']
+        except KeyError:
+            pass
         try:
             rotation = params['rotation']
         except KeyError:
@@ -202,7 +206,7 @@ def _parse_legacy_area_file(area_file_name, *regions):
     in_area = False
     for line in area_file:
         if not in_area:
-            if 'REGION' in line:
+            if 'REGION' in line and not line.strip().startswith('#'):
                 area_id = line.replace('REGION:', ''). \
                     replace('{', '').strip()
                 if area_id in area_list or select_all_areas:
@@ -210,11 +214,14 @@ def _parse_legacy_area_file(area_file_name, *regions):
                     area_content = ''
         elif '};' in line:
             in_area = False
+            try:
                 if select_all_areas:
                     area_defs.append(_create_area(area_id, area_content))
                 else:
                     area_defs[area_list.index(area_id)] = _create_area(area_id,
                                                                        area_content)
+            except KeyError:
+                raise ValueError('Invalid area definition: %s, %s' % (area_id, area_content))
         else:
             area_content += line
 
@@ -298,6 +305,110 @@ def get_area_def(area_id, area_name, proj_id, proj4_args, x_size, y_size,
                           x_size, y_size, area_extent)
 
 
+def _get_area_def_from_gdal(dataset, area_id=None, name=None, proj_id=None, proj_dict=None):
+    from pyresample.geometry import AreaDefinition
+
+    # a: width of a pixel
+    # b: row rotation (typically zero)
+    # c: x-coordinate of the upper-left corner of the upper-left pixel
+    # d: column rotation (typically zero)
+    # e: height of a pixel (typically negative)
+    # f: y-coordinate of the of the upper-left corner of the upper-left pixel
+    c, a, b, f, d, e = dataset.GetGeoTransform()
+    if not (b == d == 0):
+        raise ValueError('Rotated rasters are not supported at this time.')
+    area_extent = (c, f + e * dataset.RasterYSize, c + a * dataset.RasterXSize, f)
+
+    if proj_dict is None:
+        from osgeo import osr
+        proj = dataset.GetProjection()
+        if proj != '':
+            sref = osr.SpatialReference(wkt=proj)
+            proj_dict = proj4_str_to_dict(sref.ExportToProj4())
+        else:
+            raise ValueError('The source raster is not gereferenced, please provide the value of proj_dict')
+
+        if proj_id is None:
+            proj_id = proj.split('"')[1]
+
+    area_def = AreaDefinition(area_id, name, proj_id, proj_dict,
+                              dataset.RasterXSize, dataset.RasterYSize, area_extent)
+    return area_def
+
+
+def _get_area_def_from_rasterio(dataset, area_id, name, proj_id=None, proj_dict=None):
+    from pyresample.geometry import AreaDefinition
+
+    a, b, c, d, e, f, _, _, _ = dataset.transform
+    if not (b == d == 0):
+        raise ValueError('Rotated rasters are not supported at this time.')
+
+    if proj_dict is None:
+        crs = dataset.crs
+        if crs is not None:
+            proj_dict = dataset.crs.to_dict()
+        else:
+            raise ValueError('The source raster is not gereferenced, please provide the value of proj_dict')
+
+        if proj_id is None:
+            proj_id = crs.wkt.split('"')[1]
+
+    area_def = AreaDefinition(area_id, name, proj_id, proj_dict,
+                              dataset.width, dataset.height, dataset.bounds)
+    return area_def
+
+
+def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, proj_dict=None):
+    """Construct AreaDefinition object from raster
+
+    Parameters
+    ----------
+    source : str, Dataset, DatasetReader or DatasetWriter
+        A file name. Also it can be ``osgeo.gdal.Dataset``,
+        ``rasterio.io.DatasetReader`` or ``rasterio.io.DatasetWriter``
+    area_id : str, optional
+        ID of area
+    name : str, optional
+        Name of area
+    proj_id : str, optional
+        ID of projection
+    proj_dict : dict, optional
+        PROJ.4 parameters
+
+    Returns
+    -------
+    area_def : object
+        AreaDefinition object
+    """
+    try:
+        import rasterio
+    except ImportError:
+        rasterio = None
+        try:
+            from osgeo import gdal
+        except ImportError:
+            raise ImportError('Either rasterio or gdal must be available')
+
+    cleanup_gdal = cleanup_rasterio = None
+    if isinstance(source, (str, six.text_type)):
+        if rasterio is not None:
+            source = rasterio.open(source)
+            cleanup_rasterio = True
+        else:
+            source = gdal.Open(source)
+            cleanup_gdal = True
+
+    try:
+        if rasterio is not None and isinstance(source, (rasterio.io.DatasetReader, rasterio.io.DatasetWriter)):
+            return _get_area_def_from_rasterio(source, area_id, name, proj_id, proj_dict)
+        return _get_area_def_from_gdal(source, area_id, name, proj_id, proj_dict)
+    finally:
+        if cleanup_rasterio:
+            source.close()
+        elif cleanup_gdal:
+            source = None
+
+
 def generate_quick_linesample_arrays(source_area_def, target_area_def,
                                      nprocs=1):
     """Generate linesample arrays for quick grid resampling
@@ -410,7 +521,7 @@ def fwhm2sigma(fwhm):
 
 def convert_proj_floats(proj_pairs):
     """Convert PROJ.4 parameters to floats if possible."""
-    proj_dict = {}
+    proj_dict = OrderedDict()
     for x in proj_pairs:
         if len(x) == 1 or x[1] is True:
             proj_dict[x[0]] = True
@@ -429,10 +540,10 @@ def _get_proj4_args(proj4_args):
     """
 
     if isinstance(proj4_args, (str, six.text_type)):
-        proj_config = ConfigObj(str(proj4_args).replace('+', '').split())
+        proj_config = proj4_str_to_dict(str(proj4_args))
     else:
         proj_config = ConfigObj(proj4_args)
-    return convert_proj_floats(proj_config.dict().items())
+    return convert_proj_floats(proj_config.items())
 
 
 def proj4_str_to_dict(proj4_str):
@@ -570,3 +681,14 @@ def recursive_dict_update(d, u):
         else:
             d[k] = u[k]
     return d
+
+
+def convert_def_to_yaml(def_area_file, yaml_area_file):
+    """Convert a legacy area def file to the yaml counter partself.
+
+    *yaml_area_file* will be overwritten by the operation.
+    """
+    areas = _parse_legacy_area_file(def_area_file)
+    with open(yaml_area_file, 'w') as yaml_file:
+        for area in areas:
+            yaml_file.write(area.create_areas_def())

pyresample/version.py

@@ -16,4 +16,4 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-__version__ = '1.10.2'
+__version__ = '1.10.3'

setup.py

@@ -31,9 +31,10 @@ requirements = ['setuptools>=3.2', 'pyproj>=1.9.5.1', 'numpy>=1.10.0', 'configob
 extras_require = {'pykdtree': ['pykdtree>=1.1.1'],
                   'numexpr': ['numexpr'],
                   'quicklook': ['matplotlib', 'cartopy', 'pillow'],
+                  'rasterio': ['rasterio'],
                   'dask': ['dask>=0.16.1']}
 
-test_requires = []
+test_requires = ['rasterio']
 if sys.version_info < (3, 3):
     test_requires.append('mock')
 if sys.version_info < (2, 6):