pyshp-2.2.0 test failure

[manisandro]

I'm building proj9 for fedora, and rebuilding cartopy against it yielded the following test failure:

=================================== FAILURES ===================================
___________________________ TestLakes.test_geometry ____________________________
[gw0] linux -- Python 3.10.2 /usr/bin/python3

self = <cartopy.tests.test_shapereader.TestLakes object at 0x7f7df66ec190>

    def test_geometry(self):
        lake_geometry = self.test_lake_geometry
        assert lake_geometry.type == 'Polygon'
    
        polygon = lake_geometry
    
        expected = np.array([(-84.85548682324658, 11.147898667846633),
                            (-85.29013729525353, 11.176165676310276),
                            (-85.79132117383625, 11.509737046754324),
                            (-85.8851655748783, 11.900100816287136),
                            (-85.5653401354239, 11.940330918826362),
                            (-85.03684526237491, 11.5216484643976),
                            (-84.85548682324658, 11.147898667846633),
                            (-84.85548682324658, 11.147898667846633)])
    
>       assert_array_almost_equal(expected, polygon.exterior.coords)
E       AssertionError: 
E       Arrays are not almost equal to 6 decimals
E       
E       Mismatched elements: 12 / 16 (75%)
E       Max absolute difference: 0.75447591
E       Max relative difference: 0.00887234
E        x: array([[-84.855487,  11.147899],
E              [-85.290137,  11.176166],
E              [-85.791321,  11.509737],...
E        y: array([[-84.855487,  11.147899],
E              [-84.855487,  11.147899],
E              [-85.036845,  11.521648],...

../../BUILDROOT/python-cartopy-0.20.1-1.fc37.x86_64/usr/lib64/python3.10/site-packages/cartopy/tests/test_shapereader.py:45: AssertionError
---------------------------- Captured stdout setup -----------------------------

[manisandro]

Actually this looks like a bug in pyshp-2.2.0. With pyshp-2.1.3 the test passes.

[manisandro]

So upon further investigation, pyshp-2.2.0 enforces the positive orientation of the outer polygon ring, while pyshp-2.1.3 does not. The lake_geometry in the test is negatively oriented, and the ring of expected points is also negatively orientated. Since pyshp-2.2.0 reverses the orientation of the ring, assert_array_almost_equal(expected, polygon.exterior.coords) will fail. See also GeospatialPython/pyshp@298d981

[manisandro]

I just noticed there's also the issue that depending on whether pyshp or fiona are used to read the shapefiles, the returned ring orientations will differ, as fiona won't reverse the orientation. Not sure what convention cartopy expects, but it might also something that needs to be solved by shapely (which deals with transforming the geometry from pyshp into a shapely.geometry which is what cartopy ends up using AFAICS.

CC @QuLogic as Fedora cartopy maintainer.

[QuLogic]

Yes, I commented on that commit; I'm still confused as to whether this change is intentional.

[karimbahgat]

The switch to enforcing ring orientation in pyshp 2.2 was indeed intentional and in an effort to conform to the new GeoJSON RFC7946 specification, where it is one of the requirements. In theory the switch should not have caused any types of issues, since the original 2008 GeoJSON specification did not say anything about ring orientation. The exception as evidenced in this issue is for test cases where expected coordinate sequences have to be made explicit.

However, after reading more up on the GeoJSON RFC7946 spec, I have come to conclude that the next version pyshp will probably revert back to not enforcing ring orientation as before (instead keeping the rings in whichever orientation they were stored in the shapefile). I had not realized this before, but RFC7946 also requires that all coordinates must be in WGS84 coordinates, and some other rules about geometry splitting around the antimeridian. Supporting these requirements would then mean that pyshp would have to be dependent on pyproj in order to do projection transformations for non-WGS84 shapefiles and possibly shapely to do clipping where necessary. This would also make things much slower.

Indeed, the reason fiona does not enforce ring orientation is that it relies on GDAL to generate the geojson, and the GDAL GeoJSON driver still uses the 2008 specification by default. Given the backwards incompatible changes, the consensus in the GDAL community appears to have been that the switch to RFC7946 should be manually set and should be "all-or-nothing". Since pyshp cannot reasonably be expected to support "all" of RFC7946, I think it's better to only support the 2008 GeoJSON spec and just make this clear in the docs. After all, the purpose of pyshp is not to generate perfectly valid geojson, but rather to have quick access to the geometry contents in a human-readable standard - for these purposes 2008 GeoJSON is sufficient enough. If anyone wants to generate RFC7946-valid GeoJSON for file-based storage or sending across the internet, the burden should then fall on the user or some third-party library to do the necessary coordinate projection, clipping, etc.

While this means that cartopy's test error will disappear when I revert this in the next pyshp version, I do suspect this issue might come up again soon. As the only official standard, it's likely that RFC7946 will only become more widely used, and if GDAL makes RFC7946 the default in future versions then fiona will produce geojson which enforces ring orientation as pyshp 2.2 does, and the cartopy test will fail once again. Given this, maybe a more robust fix here is to enforce a particular ring orientation in the test itself? This might be as simple as sorting the rings before comparing them: assert_array_almost_equal(sorted(expected), sorted(polygon.exterior.coords))?

[QuLogic]

Thanks for clarifying things @karimbahgat.