Merge branch 'master' into feature-from_epsg

docs/source/conf.py

@@ -236,5 +236,6 @@ def __getattr__(cls, name):
     'pyresample': ('https://pyresample.readthedocs.io/en/stable', None),
     'trollsift': ('https://trollsift.readthedocs.io/en/stable', None),
     'trollimage': ('https://trollimage.readthedocs.io/en/stable', None),
+    'pyproj': ('https://pyproj4.github.io/pyproj/dev/', None),
     'proj4': ('https://proj.org', None),
 }

docs/source/geometry_utils.rst

@@ -22,7 +22,8 @@ necessary to get an area's ``shape`` and ``area_extent``.
 The ``create_area_def`` function has the following required arguments:
 
 * **area_id**: ID of area
-* **projection**: Projection parameters as a proj4_dict or proj4_string
+* **projection**: Projection parameters as a dictionary or string of PROJ
+  parameters.
 
 and optional arguments:
 
@@ -61,6 +62,13 @@ and optional arguments:
  Number of rows: 425
  Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
 
+.. note::
+
+    Projection (CRS) information is stored internally using the pyproj
+    library's :class:`CRS <pyproj.crs.CRS>` object. To meet certain standards
+    for representing CRS information, pyproj may rename parameters or use
+    completely different parameters from what you provide.
+
 The ``create_area_def`` function accepts some parameters in multiple forms
 to make it as easy as possible. For example, the **resolution** and **radius**
 keyword arguments can be specified with one value if ``dx == dy``:
@@ -92,7 +100,7 @@ the mercator projection with radius and resolution defined in degrees.
  >>> print(area_def)
  Area ID: ease_sh
  Description: Antarctic EASE grid
- Projection: {'a': '6371228.0', 'lat_0': '0', 'lon_0': '0', 'proj': 'merc', 'type': 'crs', 'units': 'm'}
+ Projection: {'R': '6371228', 'k': '1', 'lon_0': '0', 'no_defs': 'None', 'proj': 'merc', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
  Number of columns: 425
  Number of rows: 425
  Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
@@ -103,15 +111,15 @@ can be obtained as follows:
 .. doctest::
 
  >>> area_def = create_area_def('my_area',
- ...                            {'proj': 'latlong', 'lon_0': 0},
+ ...                            {'proj': 'longlat', 'datum': 'WGS84'},
  ...                            area_extent=[-180, -90, 180, 90],
  ...                            resolution=1,
  ...                            units='degrees',
  ...                            description='Global 1x1 degree lat-lon grid')
  >>> print(area_def)
  Area ID: my_area
  Description: Global 1x1 degree lat-lon grid
- Projection: {'lon_0': '0', 'proj': 'latlong', 'type': 'crs'}
+ Projection: {'datum': 'WGS84', 'no_defs': 'None', 'proj': 'longlat', 'type': 'crs'}
  Number of columns: 360
  Number of rows: 180
  Area extent: (-180.0, -90.0, 180.0, 90.0)

pyresample/bucket/__init__.py

@@ -90,24 +90,17 @@ def __init__(self, target_area, source_lons, source_lats):
         self._get_indices()
         self.counts = None
 
-    def _get_proj_coordinates(self, lons, lats, x_res, y_res):
-        """Calculate projection coordinates and round to resolution unit.
+    def _get_proj_coordinates(self, lons, lats):
+        """Calculate projection coordinates.
 
         Parameters
         ----------
         lons : Numpy or Dask array
             Longitude coordinates
         lats : Numpy or Dask array
             Latitude coordinates
-        x_res : float
-            Resolution of the output in X direction
-        y_res : float
-            Resolution of the output in Y direction
         """
         proj_x, proj_y = self.prj(lons, lats)
-        proj_x = round_to_resolution(proj_x, x_res)
-        proj_y = round_to_resolution(proj_y, y_res)
-
         return np.stack((proj_x, proj_y))
 
     def _get_indices(self):
@@ -121,30 +114,20 @@ def _get_indices(self):
             Y indices of the target grid where the data are put
         """
         LOG.info("Determine bucket resampling indices")
-        adef = self.target_area
 
+        # Transform source lons/lats to target projection coordinates x/y
         lons = self.source_lons.ravel()
         lats = self.source_lats.ravel()
-
-        # Create output grid coordinates in projection units
-        x_res = (adef.area_extent[2] - adef.area_extent[0]) / adef.width
-        y_res = (adef.area_extent[3] - adef.area_extent[1]) / adef.height
-        x_vect = da.arange(adef.area_extent[0] + x_res / 2.,
-                           adef.area_extent[2] - x_res / 2., x_res)
-        # Orient so that 0-meridian is pointing down
-        y_vect = da.arange(adef.area_extent[3] - y_res / 2.,
-                           adef.area_extent[1] + y_res / 2.,
-                           -y_res)
-
-        result = da.map_blocks(self._get_proj_coordinates, lons,
-                               lats, x_res, y_res,
+        result = da.map_blocks(self._get_proj_coordinates, lons, lats,
                                new_axis=0, chunks=(2,) + lons.chunks)
         proj_x = result[0, :]
         proj_y = result[1, :]
 
-        # Calculate array indices
-        x_idxs = ((proj_x - np.min(x_vect)) / x_res).astype(np.int)
-        y_idxs = ((np.max(y_vect) - proj_y) / y_res).astype(np.int)
+        # Calculate array indices. Orient so that 0-meridian is pointing down.
+        adef = self.target_area
+        x_res, y_res = adef.resolution
+        x_idxs = da.floor((proj_x - adef.area_extent[0]) / x_res).astype(np.int)
+        y_idxs = da.floor((adef.area_extent[3] - proj_y) / y_res).astype(np.int)
 
         # Get valid index locations
         mask = ((x_idxs >= 0) & (x_idxs < adef.width) &

pyresample/geometry.py

@@ -1024,8 +1024,9 @@ class AreaDefinition(BaseDefinition):
         Human-readable description of the area
     proj_id : str
         ID of projection
-    projection: dict or str
-        Dictionary or string of Proj.4 parameters
+    projection: dict or str or pyproj.crs.CRS
+        Dictionary of PROJ parameters or string of PROJ or WKT parameters.
+        Can also be a :class:`pyproj.crs.CRS` object.
     width : int
         x dimension in number of pixels, aka number of grid columns
     height : int
@@ -1077,8 +1078,22 @@ class AreaDefinition(BaseDefinition):
     pixel_offset_y : float
         y offset between projection center and upper left corner of upper
         left pixel in units of pixels..
+    crs : pyproj.crs.CRS
+        Coordinate reference system object similar to the PROJ parameters in
+        `proj_dict` and `proj_str`. This is the preferred attribute to use
+        when working with the `pyproj` library. Note, however, that this
+        object is not thread-safe and should not be passed between threads.
+    crs_wkt : str
+        WellKnownText version of the CRS object. This is the preferred
+        way of describing CRS information as a string.
+    proj_dict : dict
+        Projection defined as a dictionary of PROJ parameters. This is no
+        longer the preferred way of describing CRS information. Switch to
+        the `crs` or `crs_wkt` properties for the most flexibility.
     proj_str : str
-        Projection defined as Proj.4 string
+        Projection defined as a PROJ string. This is no
+        longer the preferred way of describing CRS information. Switch to
+        the `crs` or `crs_wkt` properties for the most flexibility.
     cartesian_coords : object
         Grid cartesian coordinates
     projection_x_coords : object
@@ -1113,8 +1128,9 @@ def __init__(self, area_id, description, proj_id, projection, width, height,
         self.pixel_size_y = (area_extent[3] - area_extent[1]) / float(height)
         self.area_extent = tuple(area_extent)
         if CRS is not None:
-            self.crs = CRS(projection)
+            self.crs_wkt = CRS(projection).to_wkt()
             self._proj_dict = None
+            self.crs = self._crs  # see _crs property for details
         else:
             if isinstance(projection, str):
                 proj_dict = proj4_str_to_dict(projection)
@@ -1149,6 +1165,23 @@ def __init__(self, area_id, description, proj_id, projection, width, height,
 
         self.dtype = dtype
 
+    @property
+    def _crs(self):
+        """Helper property for the `crs` property.
+
+        The :class:`pyproj.crs.CRS` object is not thread-safe. To avoid
+        accidentally passing it between threads, we only create it when it
+        is requested (the `self.crs` property). The alternative of storing it
+        as a normal instance attribute could cause issues between threads.
+
+        For backwards compatibility, we only create the `.crs` property if
+        pyproj 2.0+ is installed. Users can then check
+        `hasattr(area_def, 'crs')` to easily support older versions of
+        pyresample and pyproj.
+
+        """
+        return CRS.from_wkt(self.crs_wkt)
+
     @property
     def proj_dict(self):
         """Return the projection dictionary."""
@@ -1896,7 +1929,7 @@ def get_lonlats(self, nprocs=None, data_slice=None, cache=False, dtype=None, chu
             raise ValueError("Can't specify 'nprocs' and 'chunks' at the same time")
 
         # Proj.4 definition of target area projection
-        proj_def = self.crs if hasattr(self, 'crs') else self.proj_dict
+        proj_def = self.crs_wkt if hasattr(self, 'crs_wkt') else self.proj_dict
         if hasattr(target_x, 'chunks'):
             # we are using dask arrays, map blocks to th
             from dask.array import map_blocks

pyresample/kd_tree.py

@@ -79,7 +79,7 @@ def resample_nearest(source_geo_def,
     epsilon : float, optional
         Allowed uncertainty in meters. Increasing uncertainty
         reduces execution time
-    fill_value : int or None, optional
+    fill_value : int, float, numpy floating, numpy integer or None, optional
             Set undetermined pixels to this value.
             If fill_value is None a masked array is returned
             with undetermined pixels masked
@@ -617,7 +617,7 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
         If only one channel is resampled weight_funcs is
         a single function object.
         Must be supplied when using 'custom' resample type
-    fill_value : int or None, optional
+    fill_value : int, float, numpy floating, numpy integer or None, optional
         Set undetermined pixels to this value.
         If fill_value is None a masked array is returned
         with undetermined pixels masked
@@ -678,9 +678,6 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
         raise ValueError('weight_funcs must be supplied when using '
                          'custom resampling')
 
-    if not isinstance(fill_value, (long, int, float)) and fill_value is not None:
-        raise TypeError('fill_value must be number or None')
-
     if index_array.ndim == 1:
         neighbours = 1
     else:

pyresample/test/test_bucket.py

@@ -5,6 +5,7 @@
 import xarray as xr
 from unittest.mock import MagicMock, patch
 
+from pyresample import create_area_def
 from pyresample.geometry import AreaDefinition
 from pyresample import bucket
 from pyresample.test.utils import CustomScheduler
@@ -19,7 +20,6 @@ class Test(unittest.TestCase):
                            'lon_0': '0.0',
                            'proj': 'stere'}, 2560, 2048,
                           (-3780000.0, -7644000.0, 3900000.0, -1500000.0))
-
     chunks = 2
     lons = da.from_array(np.array([[25., 25.], [25., 25.]]),
                          chunks=chunks)
@@ -70,14 +70,12 @@ def test_get_proj_coordinates(self):
         prj.return_value = ([3.1, 3.1, 3.1], [4.8, 4.8, 4.8])
         lons = [1., 1., 1.]
         lats = [2., 2., 2.]
-        x_res, y_res = 0.5, 0.5
         self.resampler.prj = prj
-        result = self.resampler._get_proj_coordinates(lons, lats, x_res, y_res)
+        result = self.resampler._get_proj_coordinates(lons, lats)
         prj.assert_called_once_with(lons, lats)
         self.assertTrue(isinstance(result, np.ndarray))
-        self.assertEqual(result.shape, (2, 3))
-        self.assertTrue(np.all(result == np.array([[3., 3., 3.],
-                                                   [5., 5., 5.]])))
+        np.testing.assert_equal(result, np.array([[3.1, 3.1, 3.1],
+                                                  [4.8, 4.8, 4.8]]))
 
     def test_get_bucket_indices(self):
         """Test calculation of array indices."""
@@ -86,8 +84,28 @@ def test_get_bucket_indices(self):
             self.resampler._get_indices()
         x_idxs, y_idxs = da.compute(self.resampler.x_idxs,
                                     self.resampler.y_idxs)
-        self.assertTrue(np.all(x_idxs == np.array([1709, 1709, 1706, 1705])))
-        self.assertTrue(np.all(y_idxs == np.array([465, 465, 458, 455])))
+        np.testing.assert_equal(x_idxs, np.array([1710, 1710, 1707, 1705]))
+        np.testing.assert_equal(y_idxs, np.array([465, 465, 459, 455]))
+
+        # Additional small test case
+        adef = create_area_def(
+            area_id='test',
+            projection={'proj': 'latlong'},
+            width=2, height=2,
+            center=(0, 0),
+            resolution=10)
+        lons = da.from_array(
+            np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, -10.1, 0]),
+            chunks=2)
+        lats = da.from_array(
+            np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, 0, 10.1]),
+            chunks=2)
+        resampler = bucket.BucketResampler(source_lats=lats,
+                                           source_lons=lons,
+                                           target_area=adef)
+        resampler._get_indices()
+        np.testing.assert_equal(resampler.x_idxs, np.array([-1, 0, 0, 1, 1, 1, -1, -1, -1]))
+        np.testing.assert_equal(resampler.y_idxs, np.array([-1, 1, 1, 1, 0, 0, -1, -1, -1]))
 
     def test_get_sum(self):
         """Test drop-in-a-bucket sum."""

pyresample/test/test_files/areas.yaml

@@ -210,8 +210,7 @@ test_latlong:
   description: Basic latlong grid
   projection:
     proj: longlat
-    lat_0: 27.12
-    lon_0: -81.36
+    pm: -81.36
     ellps: WGS84
   shape:
     height: 4058