diff --git a/docs/source/conf.py b/docs/source/conf.py index abb61e2e7..34d7e672a 100644 --- a/docs/source/conf.py +++ b/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), } diff --git a/docs/source/geometry_utils.rst b/docs/source/geometry_utils.rst index 7bb5cc182..6f3486a7a 100644 --- a/docs/source/geometry_utils.rst +++ b/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 ` 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,7 +111,7 @@ 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', @@ -111,7 +119,7 @@ can be obtained as follows: >>> 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) diff --git a/pyresample/bucket/__init__.py b/pyresample/bucket/__init__.py index 9e87678fc..9b9768610 100644 --- a/pyresample/bucket/__init__.py +++ b/pyresample/bucket/__init__.py @@ -90,8 +90,8 @@ 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 ---------- @@ -99,15 +99,8 @@ def _get_proj_coordinates(self, lons, lats, x_res, y_res): 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) & diff --git a/pyresample/geometry.py b/pyresample/geometry.py index a0e5a27cd..6e6e63487 100644 --- a/pyresample/geometry.py +++ b/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 diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py index bd72ca37c..5f3e11a81 100644 --- a/pyresample/kd_tree.py +++ b/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: diff --git a/pyresample/test/test_bucket.py b/pyresample/test/test_bucket.py index bf1583944..c936bbfc3 100644 --- a/pyresample/test/test_bucket.py +++ b/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.""" diff --git a/pyresample/test/test_files/areas.yaml b/pyresample/test/test_files/areas.yaml index fa9e3255c..4859d6e8c 100644 --- a/pyresample/test/test_files/areas.yaml +++ b/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 diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py index acbe2e899..aea2a7a93 100644 --- a/pyresample/test/test_geometry.py +++ b/pyresample/test/test_geometry.py @@ -36,6 +36,11 @@ from unittest.mock import MagicMock, patch import unittest +try: + from pyproj import CRS +except ImportError: + CRS = None + class Test(unittest.TestCase): """Unit testing the geometry and geo_filter modules.""" @@ -166,7 +171,7 @@ def test_create_areas_def(self): 'areaD', {'a': '6378144.0', 'b': '6356759.0', - 'lat_0': '50.00', + 'lat_0': '90.00', 'lat_ts': '50.00', 'lon_0': '8.00', 'proj': 'stere'}, @@ -179,13 +184,25 @@ def test_create_areas_def(self): res = yaml.safe_load(area_def.create_areas_def()) expected = yaml.safe_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' + ' lat_0: 90.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) + self.assertEqual(set(res.keys()), set(expected.keys())) + res = res['areaD'] + expected = expected['areaD'] + self.assertEqual(set(res.keys()), set(expected.keys())) + self.assertEqual(res['description'], expected['description']) + self.assertEqual(res['shape'], expected['shape']) + self.assertEqual(res['area_extent']['lower_left_xy'], + expected['area_extent']['lower_left_xy']) + # pyproj versions may effect how the PROJ is formatted + for proj_key in ['a', 'lat_0', 'lon_0', 'proj', 'lat_ts']: + self.assertEqual(res['projection'][proj_key], + expected['projection'][proj_key]) + # EPSG projections = {'+init=epsg:3006': 'init: epsg:3006'} @@ -1112,34 +1129,39 @@ def test_proj_str(self): """Test the 'proj_str' property of AreaDefinition.""" from collections import OrderedDict from pyresample import utils + from pyresample.test.utils import friendly_crs_equal # pyproj 2.0+ adds a +type=crs parameter - extra_params = ' +type=crs' if utils.is_pyproj2() else '' proj_dict = OrderedDict() proj_dict['proj'] = 'stere' proj_dict['a'] = 6378144.0 proj_dict['b'] = 6356759.0 - proj_dict['lat_0'] = 50.00 + proj_dict['lat_0'] = 90.00 proj_dict['lat_ts'] = 50.00 proj_dict['lon_0'] = 8.00 area = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', proj_dict, 10, 10, [-1370912.72, -909968.64, 1029087.28, 1490031.36]) - 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' + extra_params) + assert friendly_crs_equal( + '+a=6378144.0 +b=6356759.0 +lat_0=90.0 +lat_ts=50.0 ' + '+lon_0=8.0 +proj=stere', + area + ) # try a omerc projection and no_rot parameters proj_dict['proj'] = 'omerc' + proj_dict['lat_0'] = 50.0 proj_dict['alpha'] = proj_dict.pop('lat_ts') 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.assertEqual(area.proj_str, - '+a=6378144.0 +alpha=50.0 +b=6356759.0 +lat_0=50.0 ' - '+lon_0=8.0 +no_rot +proj=omerc' + extra_params) + assert friendly_crs_equal( + '+a=6378144.0 +alpha=50.0 +b=6356759.0 +lat_0=50.0 ' + '+lon_0=8.0 +no_rot +proj=omerc', + area + ) # EPSG if utils.is_pyproj2(): @@ -1291,6 +1313,55 @@ def test_from_epsg(self): (181896.3291, 6101648.0705, 1086312.942376, 7689478.3056)) + @unittest.skipIf(CRS is None, "pyproj 2.0+ required") + def test_area_def_init_projection(self): + """Test AreaDefinition with different projection definitions.""" + proj_dict = { + 'a': '6378144.0', + 'b': '6356759.0', + 'lat_0': '90.00', + 'lat_ts': '50.00', + 'lon_0': '8.00', + 'proj': 'stere' + } + crs = CRS(CRS.from_dict(proj_dict).to_wkt()) + # pass CRS object directly + area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', + crs, + 800, 800, + [-1370912.72, -909968.64000000001, + 1029087.28, 1490031.3600000001]) + self.assertEqual(crs, area_def.crs) + # PROJ dictionary + area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', + crs.to_dict(), + 800, 800, + [-1370912.72, -909968.64000000001, + 1029087.28, 1490031.3600000001]) + self.assertEqual(crs, area_def.crs) + # PROJ string + area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', + crs.to_string(), + 800, 800, + [-1370912.72, -909968.64000000001, + 1029087.28, 1490031.3600000001]) + self.assertEqual(crs, area_def.crs) + # WKT2 + area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', + crs.to_wkt(), + 800, 800, + [-1370912.72, -909968.64000000001, + 1029087.28, 1490031.3600000001]) + self.assertEqual(crs, area_def.crs) + # WKT1_ESRI + area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', + crs.to_wkt(version='WKT1_ESRI'), + 800, 800, + [-1370912.72, -909968.64000000001, + 1029087.28, 1490031.3600000001]) + # WKT1 to WKT2 has some different naming of things so this fails + # self.assertEqual(crs, area_def.crs) + class TestMakeSliceDivisible(unittest.TestCase): """Test the _make_slice_divisible.""" @@ -2079,28 +2150,36 @@ def test_get_geostationary_bbox(self): def test_get_geostationary_angle_extent(self): """Get max geostationary angles.""" geos_area = MagicMock() - geos_area.proj_dict = {'a': 6378169.00, - 'b': 6356583.80, - 'h': 35785831.00} + del geos_area.crs + geos_area.proj_dict = { + 'proj': 'geos', + 'sweep': 'x', + 'lon_0': -89.5, + 'a': 6378169.00, + 'b': 6356583.80, + 'h': 35785831.00, + 'units': 'm'} expected = (0.15185342867090912, 0.15133555510297725) - np.testing.assert_allclose(expected, geometry.get_geostationary_angle_extent(geos_area)) - geos_area.proj_dict = {'ellps': 'GRS80', - 'h': 35785831.00} - expected = (0.15185277703584374, 0.15133971368991794) + geos_area.proj_dict['a'] = 1000.0 + geos_area.proj_dict['b'] = 1000.0 + geos_area.proj_dict['h'] = np.sqrt(2) * 1000.0 - 1000.0 + expected = (np.deg2rad(45), np.deg2rad(45)) np.testing.assert_allclose(expected, geometry.get_geostationary_angle_extent(geos_area)) - geos_area.proj_dict = {'a': 1000.0, - 'b': 1000.0, - 'h': np.sqrt(2) * 1000.0 - 1000.0} - - expected = (np.deg2rad(45), np.deg2rad(45)) - + geos_area.proj_dict = { + 'proj': 'geos', + 'sweep': 'x', + 'lon_0': -89.5, + 'ellps': 'GRS80', + 'h': 35785831.00, + 'units': 'm'} + expected = (0.15185277703584374, 0.15133971368991794) np.testing.assert_allclose(expected, geometry.get_geostationary_angle_extent(geos_area)) diff --git a/pyresample/test/test_grid.py b/pyresample/test/test_grid.py index b6ecb7325..4c75e4409 100644 --- a/pyresample/test/test_grid.py +++ b/pyresample/test/test_grid.py @@ -205,6 +205,8 @@ def test_proj4_string(self): proj4_string = self.area_def.proj_str expected_string = '+a=6378144.0 +b=6356759.0 +lat_ts=50.0 +lon_0=8.0 +proj=stere +lat_0=50.0' if is_pyproj2(): - expected_string += ' +type=crs' + expected_string = '+a=6378144 +k=1 +lat_0=50 +lon_0=8 ' \ + '+no_defs +proj=stere +rf=298.253168108487 ' \ + '+type=crs +units=m +x_0=0 +y_0=0' self.assertEqual( frozenset(proj4_string.split()), frozenset(expected_string.split())) diff --git a/pyresample/test/test_kd_tree.py b/pyresample/test/test_kd_tree.py index 1e9f4e6cf..51d37a4b4 100644 --- a/pyresample/test/test_kd_tree.py +++ b/pyresample/test/test_kd_tree.py @@ -621,6 +621,30 @@ def test_nearest_from_sample(self): expected = 15874591.0 self.assertEqual(cross_sum, expected) + def test_nearest_from_sample_np_dtypes(self): + lons = np.fromfunction(lambda y, x: 3 + x, (50, 10)) + lats = np.fromfunction(lambda y, x: 75 - y, (50, 10)) + swath_def = geometry.SwathDefinition(lons=lons, lats=lats) + valid_input_index, valid_output_index, index_array, distance_array = \ + kd_tree.get_neighbour_info(swath_def, + self.area_def, + 50000, neighbours=1, segments=1) + + for dtype in [np.uint16, np.float32]: + with self.subTest(dtype): + data = np.fromfunction(lambda y, x: y * x, (50, 10)).astype(dtype) + fill_value = dtype(0.0) + res = \ + kd_tree.get_sample_from_neighbour_info('nn', (800, 800), + data.ravel(), + valid_input_index, + valid_output_index, + index_array, + fill_value=fill_value) + cross_sum = res.sum() + expected = 15874591.0 + self.assertEqual(cross_sum, expected) + def test_custom_multi_from_sample(self): def wf1(dist): return 1 - dist / 100000.0 diff --git a/pyresample/test/test_utils.py b/pyresample/test/test_utils.py index cf5dbb3d8..e7675da6d 100644 --- a/pyresample/test/test_utils.py +++ b/pyresample/test/test_utils.py @@ -19,6 +19,7 @@ import os import unittest +from unittest import mock import io import pathlib from tempfile import NamedTemporaryFile @@ -173,8 +174,8 @@ def test_area_parser_yaml(self): if is_pyproj2(): # pyproj 2.0+ adds some extra parameters - projection = ("{'ellps': 'WGS84', 'lat_0': '27.12', " - "'lon_0': '-81.36', 'proj': 'longlat', " + projection = ("{'ellps': 'WGS84', 'no_defs': 'None', " + "'pm': '-81.36', 'proj': 'longlat', " "'type': 'crs'}") else: projection = ("{'ellps': 'WGS84', 'lat_0': '27.12', " @@ -353,6 +354,14 @@ def test_proj4_radius_parameters_provided(self): np.testing.assert_almost_equal(a, 6378273) np.testing.assert_almost_equal(b, 6356889.44891) + # test again but force pyproj <2 behavior + with mock.patch.object(utils._proj4, 'CRS', None): + a, b = utils._proj4.proj4_radius_parameters( + '+proj=stere +a=6378273 +b=6356889.44891', + ) + np.testing.assert_almost_equal(a, 6378273) + np.testing.assert_almost_equal(b, 6356889.44891) + def test_proj4_radius_parameters_ellps(self): """Test proj4_radius_parameters with ellps.""" from pyresample import utils @@ -362,16 +371,33 @@ def test_proj4_radius_parameters_ellps(self): np.testing.assert_almost_equal(a, 6378137.) np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) + # test again but force pyproj <2 behavior + with mock.patch.object(utils._proj4, 'CRS', None): + a, b = utils._proj4.proj4_radius_parameters( + '+proj=stere +ellps=WGS84', + ) + np.testing.assert_almost_equal(a, 6378137.) + np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) + def test_proj4_radius_parameters_default(self): """Test proj4_radius_parameters with default parameters.""" from pyresample import utils a, b = utils._proj4.proj4_radius_parameters( - '+proj=lcc', + '+proj=lcc +lat_0=10 +lat_1=10', ) # WGS84 np.testing.assert_almost_equal(a, 6378137.) np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) + # test again but force pyproj <2 behavior + with mock.patch.object(utils._proj4, 'CRS', None): + a, b = utils._proj4.proj4_radius_parameters( + '+proj=lcc +lat_0=10 +lat_1=10', + ) + # WGS84 + np.testing.assert_almost_equal(a, 6378137.) + np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) + def test_proj4_radius_parameters_spherical(self): """Test proj4_radius_parameters in case of a spherical earth.""" from pyresample import utils @@ -381,6 +407,14 @@ def test_proj4_radius_parameters_spherical(self): np.testing.assert_almost_equal(a, 6378273.) np.testing.assert_almost_equal(b, 6378273.) + # test again but force pyproj <2 behavior + with mock.patch.object(utils._proj4, 'CRS', None): + a, b = utils._proj4.proj4_radius_parameters( + '+proj=stere +R=6378273', + ) + np.testing.assert_almost_equal(a, 6378273.) + np.testing.assert_almost_equal(b, 6378273.) + def test_convert_proj_floats(self): from collections import OrderedDict import pyresample.utils as utils diff --git a/pyresample/test/utils.py b/pyresample/test/utils.py index f7f676f23..7c34fab0e 100644 --- a/pyresample/test/utils.py +++ b/pyresample/test/utils.py @@ -26,6 +26,11 @@ import numpy as np +try: + from pyproj import CRS +except ImportError: + CRS = None + _deprecations_as_exceptions = False _include_astropy_deprecations = False AstropyDeprecationWarning = None @@ -194,3 +199,38 @@ def __call__(self, dsk, keys, **kwargs): raise RuntimeError("Too many dask computations were scheduled: " "{}".format(self.total_computes)) return dask.get(dsk, keys, **kwargs) + + +def friendly_crs_equal(expected, actual, keys=None, use_obj=True, use_wkt=True): + """Test if two projection definitions are equal. + + The main purpose of this function is to help manage differences + between pyproj versions. Depending on the version installed and used + pyresample may provide a different `proj_dict` or other similar + CRS definition. + + Args: + expected (dict, str, pyproj.crs.CRS): Expected CRS definition as + a PROJ dictionary or string or CRS object. + actual (dict, str, pyproj.crs.CRS): Actual CRS definition + keys (list): Specific PROJ parameters to look for. Only takes effect + if `use_obj` is `False`. + use_obj (bool): Use pyproj's CRS object to test equivalence. Default + is True. + use_wkt (bool): Increase likely hood of making CRS objects equal by + converting WellKnownText before converting to the final CRS + object. Requires `use_obj`. Defaults to True. + + """ + if CRS is not None and use_obj: + if hasattr(expected, 'crs'): + expected = expected.crs + if hasattr(actual, 'crs'): + actual = actual.crs + expected_crs = CRS(expected) + actual_crs = CRS(actual) + if use_wkt: + expected_crs = CRS(expected_crs.to_wkt()) + actual_crs = CRS(actual_crs.to_wkt()) + return expected_crs == actual_crs + raise NotImplementedError("""TODO""") diff --git a/pyresample/utils/_proj4.py b/pyresample/utils/_proj4.py index 9e07695b0..8e8be8456 100644 --- a/pyresample/utils/_proj4.py +++ b/pyresample/utils/_proj4.py @@ -90,6 +90,16 @@ def proj4_radius_parameters(proj4_dict): Returns: a (float), b (float): equatorial and polar radius """ + if CRS is not None: + import math + crs = CRS(proj4_dict) + a = crs.ellipsoid.semi_major_metre + b = crs.ellipsoid.semi_minor_metre + if not math.isnan(b): + return a, b + # older versions of pyproj didn't always have a valid minor radius + proj4_dict = crs.to_dict() + if isinstance(proj4_dict, str): new_info = proj4_str_to_dict(proj4_dict) else: