North Pole Intersection fails with "Coordinate values out of range"

[jgilman]

The following sql correctly results in "POINT(0 90)". select ST_AsText( ST_Intersection( ST_GeomFromText( 'POINT(0 90)'), ST_GeomFromText( 'POINT(0 90)' ) ) );

If I change the geometry to geography it fails: select ST_Intersection( ST_GeogFromText( 'POINT(0 90)'), ST_GeogFromText( 'POINT(0 90)' ) ) ; ERROR: Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type CONTEXT: SQL function "st_intersection" statement 1

[jgilman]

I tested this with a fresh build of postgis 2.0.0 on Ubuntu and I still got the same error.

[chodgson]

The following query shows that it is possible to represent the pole correctly as a geography.

select st_astext(ST_GeogFromText( 'POINT(0 90)'));

However, the actual intersection calculation is done in a projected, planar coordinate system, and then converted back to spherical coordinates - I suspect that something goes amiss in the conversion, and the result is possibly 90.00000000000001 or something like that.

The current implementation of intersection for geographies is only an approximation using this technique of re-projecting to a "good enough" local planar projection and doing the intersection calculation in the land of planes before projecting back.

[bnordgren]

Does it choose the "local planar projection" dynamically or is the projection fixed? What are the properties of the selected projection near the poles?

It may be that there is no single projection which will work for intersecting with all geographies.

[pramsey]

You can see the logic in _best_srid, however the nutshell version is: if it's above a certain latitude, it gets a polar stereographic. If it's below, it gets an appropriate UTM, and if it's too big, it gets a mercator. We could move to a fully dynamic system, with arbitrary best projections fit to the working areas, but would fall out of the SRID cache at that point, so performance would drop. And we'd still have maximum sizes of objects we could work with.

[pramsey]

Yes, this is it, as taking apart just the transformation step shows:

select st_y(st_transform(st_transform(ST_GeomFromText( 'SRID=4326;POINT(0 90)'), _st_bestsrid(ST_GeogFromText( 'POINT(0 90)'))),4326)) - 90.0;

      ?column?       
---------------------
 2.8421709430404e-14

It shouldn't, but proj is creating the returned point every so slightly above 90.

[pramsey]

And the example works in 2.0, because of the "bump into range" code added there that nudges out of range data into range.

[pramsey]

Backported the nudger from 2.0 at r9665