Skip to content

Commit

Permalink
Merge pull request #264 from djhoese/bugfix-crs-threadsafe
Browse files Browse the repository at this point in the history 
Switch to storing CRS WKT in AreaDefinitions instead of the CRS object
  • Loading branch information
@mraspaud
mraspaud authored May 4, 2020
2 parents d8b1b88 + f3f2765 commit 066da3b
Show file tree
Hide file tree
Showing 9 changed files with 246 additions and 39 deletions.
1 change: 1 addition & 0 deletions docs/source/conf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
}
16 changes: 12 additions & 4 deletions docs/source/geometry_utils.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand Down Expand Up @@ -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``:
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down
43 changes: 38 additions & 5 deletions pyresample/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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."""
Expand Down Expand Up @@ -1877,7 +1910,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
Expand Down
3 changes: 1 addition & 2 deletions pyresample/test/test_files/areas.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
128 changes: 104 additions & 24 deletions pyresample/test/test_geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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."""
Expand Down Expand Up @@ -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'},
Expand All @@ -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'}
Expand Down Expand Up @@ -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():
Expand Down Expand Up @@ -1276,6 +1298,55 @@ def test_area_def_geocentric_resolution(self):
geo_res = area_def.geocentric_resolution()
np.testing.assert_allclose(248.594116, geo_res)

@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."""
Expand Down Expand Up @@ -1663,6 +1734,7 @@ def test_swath_def_geocentric_resolution(self):
sd = SwathDefinition(xlons, xlats)
self.assertRaises(RuntimeError, sd.geocentric_resolution)


class TestStackedAreaDefinition(unittest.TestCase):
"""Test the StackedAreaDefition."""

Expand Down Expand Up @@ -2063,28 +2135,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))

Expand Down
4 changes: 3 additions & 1 deletion pyresample/test/test_grid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()))
Loading

0 comments on commit 066da3b

Please sign in to comment.