Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

OSGB tests fail without datum transformation grids available #1880

Closed
QuLogic opened this issue Sep 18, 2021 · 8 comments · Fixed by #1884
Closed

OSGB tests fail without datum transformation grids available #1880

QuLogic opened this issue Sep 18, 2021 · 8 comments · Fixed by #1884

Comments

@QuLogic
Copy link
Member

QuLogic commented Sep 18, 2021

Description

Currently, tests use conda-forge, which is on 8.0.1, but Fedora Rawhide is on 8.1.1. With that version of Proj, a few tests fail now.

Traceback

___________________________ TestCRS.test_osgb[True] ____________________________

self = <cartopy.tests.test_crs.TestCRS object at 0x7f0fd2cd1d50>, approx = True

    @pytest.mark.parametrize('approx', [True, False])
    def test_osgb(self, approx):
>       self._check_osgb(ccrs.OSGB(approx=approx))

../../BUILDROOT/python-cartopy-0.20.0-1.fc36.x86_64/usr/lib64/python3.10/site-packages/cartopy/tests/test_crs.py:73: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <cartopy.tests.test_crs.TestCRS object at 0x7f0fd2cd1d50>
osgb = <Projected CRS: +proj=tmerc +datum=OSGB36 +ellps=airy +lon_0=-2 +l ...>
Name: unknown
Axis Info [cartesian]:
- E[east]...ion:
- name: unknown
- method: Transverse Mercator
Datum: OSGB 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwich


    def _check_osgb(self, osgb):
        ll = ccrs.Geodetic()
    
        # results obtained by streetmap.co.uk.
        lat, lon = np.array([50.462023, -3.478831], dtype=np.double)
        east, north = np.array([295132.1,  63512.6], dtype=np.double)
    
        # note the handling of precision here...
>       assert_arr_almost_eq(np.array(osgb.transform_point(lon, lat, ll)),
                             np.array([east, north]),
                             1)
E       AssertionError: 
E       Arrays are not almost equal to 1 decimals
E       
E       Mismatched elements: 2 / 2 (100%)
E       Max absolute difference: 1.62307515
E       Max relative difference: 2.55551679e-05
E        x: array([295131.,  63511.])
E        y: array([295132.1,  63512.6])

../../BUILDROOT/python-cartopy-0.20.0-1.fc36.x86_64/usr/lib64/python3.10/site-packages/cartopy/tests/test_crs.py:56: AssertionError
___________________________ TestCRS.test_osgb[False] ___________________________

self = <cartopy.tests.test_crs.TestCRS object at 0x7f0fd14d8550>, approx = False

    @pytest.mark.parametrize('approx', [True, False])
    def test_osgb(self, approx):
>       self._check_osgb(ccrs.OSGB(approx=approx))

../../BUILDROOT/python-cartopy-0.20.0-1.fc36.x86_64/usr/lib64/python3.10/site-packages/cartopy/tests/test_crs.py:73: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <cartopy.tests.test_crs.TestCRS object at 0x7f0fd14d8550>
osgb = <Projected CRS: +proj=tmerc +datum=OSGB36 +ellps=airy +lon_0=-2 +l ...>
Name: unknown
Axis Info [cartesian]:
- E[east]...ion:
- name: unknown
- method: Transverse Mercator
Datum: OSGB 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwich


    def _check_osgb(self, osgb):
        ll = ccrs.Geodetic()
    
        # results obtained by streetmap.co.uk.
        lat, lon = np.array([50.462023, -3.478831], dtype=np.double)
        east, north = np.array([295132.1,  63512.6], dtype=np.double)
    
        # note the handling of precision here...
>       assert_arr_almost_eq(np.array(osgb.transform_point(lon, lat, ll)),
                             np.array([east, north]),
                             1)
E       AssertionError: 
E       Arrays are not almost equal to 1 decimals
E       
E       Mismatched elements: 2 / 2 (100%)
E       Max absolute difference: 1.62307537
E       Max relative difference: 2.55551713e-05
E        x: array([295131.,  63511.])
E        y: array([295132.1,  63512.6])

../../BUILDROOT/python-cartopy-0.20.0-1.fc36.x86_64/usr/lib64/python3.10/site-packages/cartopy/tests/test_crs.py:56: AssertionError
______________________________ TestCRS.test_epsg _______________________________

self = <cartopy.tests.test_crs.TestCRS object at 0x7f0fd2c35ae0>

    def test_epsg(self):
        uk = ccrs.epsg(27700)
        assert uk.epsg_code == 27700
        assert_almost_equal(uk.x_limits, (-104009.357,  688806.007), decimal=3)
        assert_almost_equal(uk.y_limits, (-8908.37, 1256558.45), decimal=2)
        assert_almost_equal(uk.threshold, 7928.15, decimal=2)
>       self._check_osgb(uk)

../../BUILDROOT/python-cartopy-0.20.0-1.fc36.x86_64/usr/lib64/python3.10/site-packages/cartopy/tests/test_crs.py:81: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = <cartopy.tests.test_crs.TestCRS object at 0x7f0fd2c35ae0>
osgb = _EPSGProjection(27700)

    def _check_osgb(self, osgb):
        ll = ccrs.Geodetic()
    
        # results obtained by streetmap.co.uk.
        lat, lon = np.array([50.462023, -3.478831], dtype=np.double)
        east, north = np.array([295132.1,  63512.6], dtype=np.double)
    
        # note the handling of precision here...
>       assert_arr_almost_eq(np.array(osgb.transform_point(lon, lat, ll)),
                             np.array([east, north]),
                             1)
E       AssertionError: 
E       Arrays are not almost equal to 1 decimals
E       
E       Mismatched elements: 2 / 2 (100%)
E       Max absolute difference: 1.62307537
E       Max relative difference: 2.55551713e-05
E        x: array([295131.,  63511.])
E        y: array([295132.1,  63512.6])

../../BUILDROOT/python-cartopy-0.20.0-1.fc36.x86_64/usr/lib64/python3.10/site-packages/cartopy/tests/test_crs.py:56: AssertionError

The differences are rather small, but I did not see anything obvious that might have been the cause in Proj.

Full environment definition

Operating system

Fedora Rawhide

Cartopy version

0.20.0

@QuLogic
Copy link
Member Author

QuLogic commented Sep 21, 2021

Hmm, so I tried downgrading to 8.0.0 and it still fails, so now I'm not sure what the real cause is.

@QuLogic
Copy link
Member Author

QuLogic commented Sep 21, 2021

Additionally, the expected projected values were changed in 726ebac#diff-535d7d5e8b3b49d4b6a3bee21502ba1c3835941b30f5b629e5f0ef6784faaef0R53 for the pyproj conversion, but I am getting back the values from before the change.

I can replicate the "new" values using a conda environment, but not in any online converter (i.e., the streetmap.co.uk mentioned in the test comment, or britishnationalgrid.uk).

Any ideas @snowman2?

@snowman2
Copy link
Contributor

I wonder if the datum/transformation grids were updated ref?

@QuLogic
Copy link
Member Author

QuLogic commented Sep 21, 2021

Ah, that's a good thought. The test currently transforms from/to Geodetic, which is WGS84, while OSGB is using OSGB36/Airy. If I change the test to use osgb.as_geodetic() instead and copy the OSGB36 lat/lon from streetmap.co.uk, then it passes both in conda and in Fedora builds.

There's no network available on the Fedora build, and only the main package data is installed. There are many data subpackages available, so maybe I can install one of those to get the missing transformation grid.

@QuLogic QuLogic changed the title Some CRS tests fail with Proj 8.1.1 Some CRS tests fail without datum transformation grids available Sep 21, 2021
@QuLogic QuLogic changed the title Some CRS tests fail without datum transformation grids available OSGB tests fail without datum transformation grids available Sep 21, 2021
@QuLogic
Copy link
Member Author

QuLogic commented Sep 21, 2021

So that's it; I installed proj-data-uk, which adds:

/usr/share/proj/uk_os_OSGM15_Belfast.tif
/usr/share/proj/uk_os_OSGM15_GB.tif
/usr/share/proj/uk_os_OSGM15_Malin.tif
/usr/share/proj/uk_os_OSTN15_NTv2_OSGBtoETRS.tif
/usr/share/proj/uk_os_README.txt

and now the tests pass without any patching.

Is there a way to determine if the necessary grids are available? Should we just tweak the test to work either way?

@snowman2
Copy link
Contributor

@QuLogic
Copy link
Member Author

QuLogic commented Sep 21, 2021

That seems good; I'll put together a PR.

For that pyproj check, shouldn't you check if files exist before checking the network flag?

@snowman2
Copy link
Contributor

shouldn't you check if files exist before checking the network flag?

When PROJ_NETWORK is enabled, it behaves as if you have all of the files downloaded. So, there is no need to check for the file. This assumes that your network works ...

QuLogic added a commit to QuLogic/cartopy that referenced this issue Sep 21, 2021
@QuLogic QuLogic added this to the 0.20.1 milestone Sep 22, 2021
QuLogic added a commit to QuLogic/cartopy that referenced this issue Oct 7, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants