cartopy.crs.Orthographic.project_geometry method returns an empty geometry

lmmarsano commented 5 years ago

Description

geoplot gives a ValueError: Failed to determine the required bounds in projection coordinates. message whose root cause traces to cartopy.crs.Orthographic.project_geometry returning an empty geometry. The unit test added in my branch exemplifies the problem. Is project_geometry correct to return an empty geometry in such circumstances?

Code to reproduce

Refer to lmmarsano/cartopy@9b133e3

Traceback

============================= test session starts ==============================
platform linux -- Python 3.7.1, pytest-3.9.2, py-1.7.0, pluggy-0.8.0
rootdir: /home/luism/project/cartopy, inifile: setup.cfg
collected 1 item                                                               

cartopy/tests/test_ortho_proj.py F                                       [100%]

=================================== FAILURES ===================================
____________________________ test_project_geometry _____________________________

    def test_project_geometry():
        """Test project_geometry to Orthographic. Fails in geoplot.kdeplot"""
        x2, y2 = 90, 1
        x1, y1 = -x2, -y2
>       assert not cp.crs.Orthographic().project_geometry(
            sgeom.polygon.LineString([[x1, y1], [x2, y1],
                                     [x2, y2], [x1, y2],
                                     [x1, y1]])
        ).is_empty
E    AssertionError: assert not True
E     +  where True = <shapely.geometry.multilinestring.MultiLineString object at 0x7fb0e993c748>.is_empty
E     +    where <shapely.geometry.multilinestring.MultiLineString object at 0x7fb0e993c748> = <bound method Projection.project_geometry of <cartopy.crs.Orthographic object at 0x7fb0ea3e2a98>>(<shapely.geometry.linestring.LineString object at 0x7fb0e993c710>)
E     +      where <bound method Projection.project_geometry of <cartopy.crs.Orthographic object at 0x7fb0ea3e2a98>> = <cartopy.crs.Orthographic object at 0x7fb0ea3e2a98>.project_geometry
E     +        where <cartopy.crs.Orthographic object at 0x7fb0ea3e2a98> = <class 'cartopy.crs.Orthographic'>()
E     +          where <class 'cartopy.crs.Orthographic'> = <module 'cartopy.crs' from '/home/luism/project/cartopy/build/lib.linux-x86_64-3.7/cartopy/crs.py'>.Orthographic
E     +            where <module 'cartopy.crs' from '/home/luism/project/cartopy/build/lib.linux-x86_64-3.7/cartopy/crs.py'> = cp.crs
E     +      and   <shapely.geometry.linestring.LineString object at 0x7fb0e993c710> = <class 'shapely.geometry.linestring.LineString'>([[-90, -1], [90, -1], [90, 1], [-90, 1], [-90, -1]])
E     +        where <class 'shapely.geometry.linestring.LineString'> = <module 'shapely.geometry.polygon' from '/home/luism/.pyenv/versions/3.7.1/envs/cartopy/lib/python3.7/site-packages/shapely/geometry/polygon.py'>.LineString
E     +          where <module 'shapely.geometry.polygon' from '/home/luism/.pyenv/versions/3.7.1/envs/cartopy/lib/python3.7/site-packages/shapely/geometry/polygon.py'> = sgeom.polygon

cartopy/tests/test_ortho_proj.py:7: AssertionError
=========================== 1 failed in 1.24 seconds ===========================

Full environment definition

### Operating system WSL (Ubuntu 18.10) ### Cartopy version v0.16.0 ### pip list ``` -r requirements/default.txt -r requirements/tests.txt Cython ```

QuLogic commented 5 years ago

Since you don't specify a src_crs, it is using a Geodetic projection as source. In Geodetic, the path between the first two points is the closest path in spherical terms, which is one going through the poles. We can see that (as opposed to PlateCarree) in this plot: figure_1-6 So, since your input is essentially a loop around the poles, for a default Orthographic with central longitude of 0, the result is right around the boundary.

Due to #589, the boundary is slightly smaller than your line string, so it is clipped. There have been a few fixes to projections, so maybe that could be reverted, but that shrunk boundary fixed a lot of weird polygon issues. I don't know if reverting that change is possible or not.

lmmarsano commented 5 years ago

Thanks. The original issue actually occurred when transforming from a PlateCarree CRS (I oversimplified, sorry). I've updated my branch to include that test case and fail from both source CRSs. Please refer to lmmarsano/cartopy@8bf9e18

Should this other test case also fail?

QuLogic commented 5 years ago

Well, since you changed +/-90 to +/-91, it's definitely true that the Geodetic should fail since it's completely outside the default Orthographic limits: figure_1-7 I'm not sure about PlateCarree; it might be related to #116.

QuLogic commented 5 years ago

Due to #589, the boundary is slightly smaller than your line string, so it is clipped. There have been a few fixes to projections, so maybe that could be reverted, but that shrunk boundary fixed a lot of weird polygon issues. I don't know if reverting that change is possible or not.

So @pelson, do you think this can be reverted? It seems pretty stable without the 0.99999, but I'm not sure how you concluded that would be the right fix.

pelson commented 5 years ago

So @pelson, do you think this can be reverted? It seems pretty stable without the 0.99999, but I'm not sure how you concluded that would be the right fix.

I'm not against it. The geometry transformation pipeline has undergone a significant amount of refinement over the last 3 years, so it is plausible that it is no longer necessary.

SciTools / cartopy