Fix support for pyproj-2 EPSG syntax #186

sfinkens · 2019-04-15T15:09:48Z

This PR fixes the creation of area defs using pyproj-2 EPSG syntax: AreaDefinition(projection='EPSG:XXXX', ...). Additionally it adds (limited) support for EPSG codes in AreaDefinition.to_cartopy_crs().

Syntax

Valid EPSG definitions are now:

AreaDefinition(projection='+init=EPSG:XXXX', ...)
AreaDefinition(projection={'init': 'EPSG:XXXX'}, ...)
AreaDefinition(projection='EPSG:XXXX', ...)

The former two cases in conjunction with pyproj 2+ will trigger a warning. For the implementation I came up with two options:

Convert EPSG codes to their formal definition:

>>> pyproj.Proj('EPSG:4326').definition_string()
'proj=longlat datum=WGS84 no_defs ellps=WGS84 towgs84=0,0,0'

Then convert that string to a proj_dict.

Add extra logic for EPSG:XXXX syntax

Option 1. would have been the cleaner solution, because we could have used a uniform proj_dict everywhere. However, I found out that the conversion from EPSG code to proj4 string is a potentially lossy transformation (see pyproj4/pyproj#207 (comment)). Furthermore I had a performance issue with pyproj-2 and definition strings. The following snippet took > 10 minutes (I cancelled after 10 minutes):

In [1]: import pyresample.geometry as geometry 
   ...: import pyproj 
   ...:  
   ...: area_def = geometry.AreaDefinition('geos', 
   ...:                                    'Geostationary test area', 
   ...:                                    'geos', 
   ...:                                    {'a': 6378169.0, 'b': 6356583.8, 'h': 35785831.0,  
   ...:                                     'lon_0': 0.0, 'proj': 'geos', 'units': 'm'}, 
   ...:                                    3712, 3712, 
   ...:                                    (-5570248.477339745, -5561247.267842293, 5567248.074173927, 5570248.477339745)) 
   ...:  
   ...: area_to_cover = geometry.AreaDefinition('epsg4326',  
   ...:                                         'Global equal latitude/longitude grid for global sphere', 
   ...:                                         'epsg4326', 
   ...:                                         pyproj.Proj('EPSG:4326').definition_string(), 
   ...:                                         8192, 
   ...:                                         4096, 
   ...:                                         [-180.0, -90.0, 180.0, 90.0])                                                                                                                                                                

In [2]: area_to_cover                                                                                                                                                                                                                        
Out[2]: 
Area ID: epsg4326
Description: Global equal latitude/longitude grid for global sphere
Projection ID: epsg4326
Projection: {'datum': 'WGS84', 'ellps': 'WGS84', 'no_defs': 'True', 'proj': 'longlat', 'towgs84': '0,0,0'}
Number of columns: 8192
Number of rows: 4096
Area extent: (-180.0, -90.0, 180.0, 90.0)

In [3]: slice_x, slice_y = area_def.get_area_slices(area_to_cover)

Without the extra towgs84 parameter in the proj_dict it finishes within a few seconds.

Cartopy CRS

Additionally, this PR adds limited support for EPSG codes in AreaDefinition.to_cartopy_crs(). Again, there were two options:

Use cartopy.crs.epsg() (https://github.com/SciTools/cartopy/blob/v0.16.0/lib/cartopy/crs.py#L1993)
Convert EPSG code to proj4 string and issue a warning that this transformation is possibly lossy

Unfortunately cartopy.crs.epsg() requires an active internet connection and doesn't work for geodetic coordinate systems such as EPSG:4326. That is why I chose option 2.

Tests added
Tests passed
Passes git diff origin/master **/*py | flake8 --diff
Fully documented

Usage: AreaDefinition(projection='EPSG:XXXX', ...)

pyresample/utils/_proj4.py

…o fix-epsg

pyresample/_spatial_mp.py

coveralls · 2019-04-15T15:53:25Z

Coverage decreased (-0.04%) to 86.585% when pulling bc861f8 on sfinkens:fix-epsg into 18807c9 on pytroll:master.

codecov · 2019-04-15T15:53:34Z

Codecov Report

Merging #186 into master will decrease coverage by 0.03%.
The diff coverage is 85.18%.

@@            Coverage Diff             @@
##           master     #186      +/-   ##
==========================================
- Coverage   86.65%   86.61%   -0.04%     
==========================================
  Files          36       37       +1     
  Lines        6653     6776     +123     
==========================================
+ Hits         5765     5869     +104     
- Misses        888      907      +19

Impacted Files	Coverage Δ
pyresample/test/__init__.py	`95.45% <100%> (+0.21%)`	⬆️
pyresample/utils/__init__.py	`84.61% <100%> (+0.61%)`	⬆️
pyresample/test/test_utils.py	`99.62% <100%> (+0.02%)`	⬆️
pyresample/utils/_proj4.py	`93.33% <100%> (+1.84%)`	⬆️
pyresample/_cartopy.py	`94.33% <100%> (+0.58%)`	⬆️
pyresample/_spatial_mp.py	`62.11% <50%> (-1.53%)`	⬇️
pyresample/test/test_spatial_mp.py	`66.66% <66.66%> (ø)`
pyresample/test/test_geometry.py	`97.06% <88.88%> (-0.68%)`	⬇️
pyresample/geometry.py	`84.43% <0%> (+0.11%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 18807c9...bc861f8. Read the comment docs.

djhoese

I'm not happy how much work has had to go in to supporting EPSG codes, but I suppose we don't have much of a choice. Another option that I'm now realizing, could/should we just pass whatever the user provided to Proj and let their deprecation warnings be the "source of truth"? I suppose with our use of the proj dicts that isn't possible.

As for your performance issues, could you maybe get this down to a simple Proj(...)(x, y) use case? Maybe file a bug with pyproj if one does not exist already?

mraspaud · 2019-04-15T18:13:35Z

I had an example that failed reading from a yaml file. I would like this to be tested if possible ?

sfinkens · 2019-04-16T07:53:02Z

@djhoese Yeah, I'm not happy either. Your suggestion to pass whatever the user provided to Proj is the cleanest solution. However that would be lots of internal changes. IIRC you mentioned moving to CRSes at some point. I think this could solve most of the problems. Maybe by replacing proj_dicts with CRSes. Many properties are accessible from the CRS object (crs.ellipsoid, crs.prime_meridian). Only something like crs.projection is currently missing.

pyresample/test/test_utils.py

pyresample/test/test_geometry.py

djhoese · 2019-04-16T13:00:53Z

pyresample/test/test_geometry.py

+                    area_def.proj_dict['a'] = area_def.proj_dict.pop('R')
+                for key in ['x_0', 'y_0', 'no_defs', 'b', 'init']:
+                    area_def.proj_dict.pop(key, None)
+


What is the above actually testing? Is this showing that area definitions with EPSG codes aren't comparable out of the box?

The generated area defs are compared to base_def whos proj_dict looks like {'proj': 'laea', 'lat_0': -90, ...}. However, the EPSG area defs have proj_dict={'EPSG': 1234} so that they are not directly comparable to base_def. The workaround is to use their definition_string. But that contains some additional (redundant) parameters like x_0=0, y_0=0. The ellipsoid parameters in the definition string are also different depending on the proj version (a and b in pyproj >= 2 vs R in pyproj < 2)

djhoese · 2019-04-16T13:07:02Z

pyresample/utils/_proj4.py

+        try:
+            code = int(proj4_str.split(':', 1)[1])
+        except ValueError as err:
+            six.raise_from(ValueError("Invalid EPSG code '{}': {}".format(proj4_str, err)),


Am I correct in thinking that this won't be needed when python 2.7 is dropped?

Yes, correct! Then we could use raise ... from None

pyresample/utils/_proj4.py

If EPSG code is present, ignore other parameters

sfinkens · 2019-04-18T09:35:53Z

@djhoese Regarding the performance issue in AreaDefinition.get_area_slices(): The root cause is that proj.crs.is_geographic = False if the projection is created using the EPSG:4326 definition_string

In [3]: pyproj.Proj('EPSG:4326').definition_string()                                                                                                                                                                                         
Out[3]: 'proj=longlat datum=WGS84 no_defs ellps=WGS84 towgs84=0,0,0'

In [4]: pyproj.Proj('proj=longlat datum=WGS84 no_defs ellps=WGS84 towgs84=0,0,0').crs.is_geographic                                                                                                                                          
Out[4]: False

See pyproj4/pyproj#274. As a consequence, the shortcut

if self.is_latlong():
    return data1, data2

in spatial_mp.Proj.__call__ doesn't work. This in turn leads to invalid corner coordinates in AreaDefinition.__init__ because the projection coordinates are expected to be radians, not degrees. And those invalid corner coordinates seem to slow down the intersection computation in AreaDefinition.get_area_slices()... 😅

Conversion proj4 string <-> dict can be tested for all pyproj versions

sfinkens and others added 3 commits April 15, 2019 14:43

Add support for pyproj-2 EPSG syntax

a7d0027

Usage: AreaDefinition(projection='EPSG:XXXX', ...)

Add EPSG support in AreaDefinition.to_cartopy_crs

54c6f3f

Fixing style errors.

f899a47

stickler-ci reviewed Apr 15, 2019

View reviewed changes

pyresample/utils/_proj4.py Outdated Show resolved Hide resolved

sfinkens added 2 commits April 15, 2019 15:11

Fix unused import

f7a8c65

Merge branch 'fix-epsg' of https://github.com/sfinkens/pyresample int…

5e2b196

…o fix-epsg

sfinkens changed the title ~~Add support for pyproj-2 EPSG syntax~~ Fix support for pyproj-2 EPSG syntax Apr 15, 2019

djhoese reviewed Apr 15, 2019

View reviewed changes

pyresample/_spatial_mp.py Outdated Show resolved Hide resolved

sfinkens added 2 commits April 15, 2019 15:25

Remove unnecessary copy of keyword arguments

3ffa4cb

Fix mock import

17a14ad

wroberts4 reviewed Apr 15, 2019

View reviewed changes

pyresample/_spatial_mp.py Show resolved Hide resolved

djhoese approved these changes Apr 15, 2019

View reviewed changes

sfinkens added 4 commits April 16, 2019 09:20

Make EPSG code an integer

6b5839f

Add more tests for proj4 str <-> dict conversions

8eae103

Test YAML parsing of EPSG area defs

8891385

Add is_pyproj2 utility

26631de

stickler-ci reviewed Apr 16, 2019

View reviewed changes

pyresample/test/test_utils.py Outdated Show resolved Hide resolved

pyresample/test/test_utils.py Outdated Show resolved Hide resolved

sfinkens added 2 commits April 16, 2019 09:43

Fix unused imports

28ea4fd

Use yaml.safe_load in YAML tests

4e25dab

djhoese reviewed Apr 16, 2019

View reviewed changes

pyresample/test/test_geometry.py Outdated Show resolved Hide resolved

djhoese reviewed Apr 16, 2019

View reviewed changes

pyresample/utils/_proj4.py Outdated Show resolved Hide resolved

sfinkens added 3 commits April 17, 2019 07:21

Update proj4_dict_to_str

e7704d8

If EPSG code is present, ignore other parameters

Test +init=EPSG, regardless of the pyproj version

b733d5e

Update documentation

e27fb30

sfinkens added 3 commits April 23, 2019 08:07

Avoid some test misses

5971b3e

Conversion proj4 string <-> dict can be tested for all pyproj versions

Add EPSG YAML example to the docs

3ddfbe1

Use new location of load_area in the docs

bc861f8

djhoese merged commit b10854b into pytroll:master Apr 23, 2019

sfinkens deleted the fix-epsg branch April 23, 2019 13:28

djhoese added the bug label Apr 25, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix support for pyproj-2 EPSG syntax #186

Fix support for pyproj-2 EPSG syntax #186

sfinkens commented Apr 15, 2019 •

edited

Loading

coveralls commented Apr 15, 2019 •

edited

Loading

codecov bot commented Apr 15, 2019 •

edited

Loading

djhoese left a comment

mraspaud commented Apr 15, 2019

sfinkens commented Apr 16, 2019

djhoese Apr 16, 2019

sfinkens Apr 16, 2019

djhoese Apr 16, 2019

sfinkens Apr 16, 2019 •

edited

Loading

sfinkens commented Apr 18, 2019

Fix support for pyproj-2 EPSG syntax #186

Fix support for pyproj-2 EPSG syntax #186

Conversation

sfinkens commented Apr 15, 2019 • edited Loading

Syntax

Cartopy CRS

coveralls commented Apr 15, 2019 • edited Loading

codecov bot commented Apr 15, 2019 • edited Loading

Codecov Report

djhoese left a comment

Choose a reason for hiding this comment

mraspaud commented Apr 15, 2019

sfinkens commented Apr 16, 2019

djhoese Apr 16, 2019

Choose a reason for hiding this comment

sfinkens Apr 16, 2019

Choose a reason for hiding this comment

djhoese Apr 16, 2019

Choose a reason for hiding this comment

sfinkens Apr 16, 2019 • edited Loading

Choose a reason for hiding this comment

sfinkens commented Apr 18, 2019

sfinkens commented Apr 15, 2019 •

edited

Loading

coveralls commented Apr 15, 2019 •

edited

Loading

codecov bot commented Apr 15, 2019 •

edited

Loading

sfinkens Apr 16, 2019 •

edited

Loading