Opened 8 years ago

Closed 7 years ago

#1170 closed defect (fixed)

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

Reported by: jgilman Owned by: pramsey
Priority: critical Milestone: PostGIS 1.5.4
Component: postgis Version: 1.5.X
Keywords: Cc:

Description

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

Change History (7)

comment:1 Changed 8 years ago by jgilman

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

comment:2 Changed 8 years ago by 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.

comment:3 Changed 8 years ago by 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.

comment:4 Changed 8 years ago by 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.

comment:5 Changed 7 years ago by 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.

comment:6 Changed 7 years ago by 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.

comment:7 Changed 7 years ago by pramsey

Resolution: fixed
Status: newclosed

Backported the nudger from 2.0 at r9665

Note: See TracTickets for help on using tickets.