cs2cs calculations fail for certain projections

[corinnebosley]

I have noticed that coordinate transforms between certain projections are failing upon calling the transform from proj4. This does not happen with all coordinates, and not with all projections.

The projections I have been having this problem with are: Orthographic Transverse Mercator OSGB Geostationary Gnomonic

To illustrate the problem I am experiencing with some coordinates but not others, here is an example:

import cartopy.crs as ccrs
import numpy as np

good_xy = np.array([[292.49993896, -77.1010721],
                    [296.24993896, -75.62961341],
                    [297.518888, -77.49995422]])

bad_xy = np.array([[142.49996948, -72.48229747],
                   [164.99995422, -72.9603603],
                   [112.99996948, -77.406934]])

for points in (good_xy, bad_xy):
    print('xy input coordinates:')
    print(points)
    result = ccrs.Orthographic().transform_points(ccrs.PlateCarree(),
                                                  points[:, 0],
                                                  points[:, 1])[:, 0:2]
    print('result of coordinate transform:')
    print(result)
    print

If this is not a bug, but it is in fact a problem with the way that we are invoking proj4's transformations, I would appreciate a chat about it so that we can configure our software to catch and correct these problems before going anywhere near the actual transforms.

Thank you, and merry christmas and things.

[hobu]

Please demonstrate this issue with only the cs2cs command and not with cartopy or numpy or anything else that can potentially confound the issue.

[corinnebosley]

Good coords:

cs2cs +over +proj=latlon +ellps=WGS84 +to +proj=ortho +ellps=WGS84 
-r <<EOF
>77.1S 68W
>EOF
-1320235.47 -6217160.44 0.00

Bad coords:

cs2cs +over +proj=latlon +ellps=WGS84 +to +proj=ortho +ellps=WGS84 
-r <<EOF
>72.5S 142.5E
>EOF
Rel. 4.9.1, 04 March 2015
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): tolerance condition error
*   * 0.00

[kbevers]

Confirmed with most recent release:

$ cs2cs +over +proj=latlon +ellps=WGS84 +to +proj=ortho +ellps=WGS84 -r <<EOF
> 72.5S 142.5E
> EOF
Rel. 4.9.3, 15 August 2016
<cs2cs>: while processing file: <stdin>, line 1
pj_transform(): tolerance condition error
*   * 0.00

Same results with cs2cs from master.

[micahcochran]

With the ortho projection, that point is in Antarctica. From my perspective it seems as if it is on the other side of the globe. Some type of out of bounds type of error might be reasonable response from proj.4.

[kbevers]

So I had a proper look at this, instead of just confirming the results with the most recent version of Proj.4, and there is actually two things at play here:

1. The inputs are reversed. Here they are latititude, longitude and they should be longitude, latitude. Can be corrected by using the -r flag.
2. As @micahcochran states above, the point is outside the defined area of the projection. The Orthographic projection only maps parts of the globe (inside +/- 90 deg longitude of the longitude origin if I am not mistaken). This can be corrected by changing the longitude origin +lon_0=140E

With these adjustments the cs2cs call now looks something like:

> echo 72.5S 142.5E | cs2cs -r +proj=latlon +to +proj=ortho +lon_0=140E
83659.49        -6082937.37 0.00

Note that +over is not needed in this case.

For the rest of the projections you mention, I bet you run into something similar.

[corinnebosley]

@kbevers @micahcochran Thanks for this, I will add some bits to our software to catch this scenario. Appreciate the help.