4

I'm quite new to GIS and PostGIS and, probably, this question has been already answered multiple times, but unfortunately I was not able to find an answer.

I have few polygons and I want to find those that are intersecting with each other using the query like:

select * from polygons where ST_Intersects(ST_GeogFromText('POLYGON ((-128 50, -65 50, -65 22, -128 22, -128 50))'), geography_column)

But sometimes this query returns wrong results, i.e. returns polygons that are not really intersecting, however when I use the ST_Intersection to find the intersecting parts PostGIS tells me that these polygons don't have any intersections.

For example, I have 2 polygons that are not intersecting 'POLYGON ((-128 50, -65 50, -65 22, -128 22, -128 50))' and 'POLYGON ((-114.77891171734832 56.331842688203174, -118.11874428779515 52.76935461305989, -112.99766767977668 51.43342158488116, -108.76721309054403 52.101388098970524, -109.21252409993694 55.2185651647209, -111.66173465159795 56.99980920229254, -114.77891171734832 56.331842688203174))'

If I try to get intersections between them as follows:

select St_AsText(
    St_Intersection(
        St_GeogFromText('POLYGON ((-128 50, -65 50, -65 22, -128 22, -128 50))'),
        St_GeogFromText('POLYGON ((-114.77891171734832 56.331842688203174, -118.11874428779515 52.76935461305989, -112.99766767977668 51.43342158488116, -108.76721309054403 52.101388098970524, -109.21252409993694 55.2185651647209, -111.66173465159795 56.99980920229254, -114.77891171734832 56.331842688203174))')
        )
    );

select St_AsText(
    St_Intersection(
        St_GeomFromText('POLYGON ((-128 50, -65 50, -65 22, -128 22, -128 50))'),
        St_GeomFromText('POLYGON ((-114.77891171734832 56.331842688203174, -118.11874428779515 52.76935461305989, -112.99766767977668 51.43342158488116, -108.76721309054403 52.101388098970524, -109.21252409993694 55.2185651647209, -111.66173465159795 56.99980920229254, -114.77891171734832 56.331842688203174))')
        )
    );

I get the GEOMETRYCOLLECTION EMPTY which seems to be correct.

But if I try to check intersections with ST_Intersects like this:

select St_Intersects(
        St_GeogFromText('POLYGON ((-128 50, -65 50, -65 22, -128 22, -128 50))'),
        St_GeogFromText('POLYGON ((-114.77891171734832 56.331842688203174, -118.11874428779515 52.76935461305989, -112.99766767977668 51.43342158488116, -108.76721309054403 52.101388098970524, -109.21252409993694 55.2185651647209, -111.66173465159795 56.99980920229254, -114.77891171734832 56.331842688203174))')
        );

I get true, however both of the following returns false as expected:

select St_Intersects(
        St_GeomFromText('POLYGON ((-128 50, -65 50, -65 22, -128 22, -128 50))'),
        St_GeomFromText('POLYGON ((-114.77891171734832 56.331842688203174, -118.11874428779515 52.76935461305989, -112.99766767977668 51.43342158488116, -108.76721309054403 52.101388098970524, -109.21252409993694 55.2185651647209, -111.66173465159795 56.99980920229254, -114.77891171734832 56.331842688203174))')
        );

select St_Intersects(
        St_GeomFromText('POLYGON ((-128 50, -65 50, -65 22, -128 22, -128 50))', 4326),
        St_GeomFromText('POLYGON ((-114.77891171734832 56.331842688203174, -118.11874428779515 52.76935461305989, -112.99766767977668 51.43342158488116, -108.76721309054403 52.101388098970524, -109.21252409993694 55.2185651647209, -111.66173465159795 56.99980920229254, -114.77891171734832 56.331842688203174))', 4326)
        );

Can someone explain this behavior and suggest what would be the correct way for finding intersecting polygons in my case?

PostgreSQL: 9.4.5

PostGIS: 2.1.8 r13780

  • When you play with geography the shortest path from one corner into another does not look like a straight line as it is does when you use geometry. Perhaps that is the reason. – user30184 Feb 12 '16 at 14:54
  • @user30184 In this case shouldn't both ST_Intersects and ST_Intersection return the same results for the same inputs, i.e. if there is no intersecting areas then ST_Intersects returns false? Also, shouldn't the geometry intersection with provided SRID give the same result as geography intersection? – Vader Feb 12 '16 at 15:06
7

The question has indeed been asked before, but it is still a good question :) The manual already gives you a hint towards the problem.

The main problem is that you are trying to use geograpy types (that work with real world coordinates) while the function st_intersection is not completely ready for geography.

In the real (spherical) world, your geometries would overlap, since the top of the big square you give would follow the circle of the earth (and show as convex arc on your 2D map). ST_Intersects knows this for geography and hence gives 'true' ST_Intersection on the other hand is less smart and uses a 2D projection to calculate where your square will remain a square. It chooses the 2D projection itself but doesn't always choose right and defaults to mercator (especially not for larger areas). There is a nice explanation overhere

Your solution will be to first transform the geography to some appropriate geometry (don't know by heart) and then do the intersection.

  • Thanks for the answer. Wasn't aware of these ST_Intersection internals. Just one more thing to clarify, in my understanding geometry with spatial reference id set to 4326 should be the same to geography, i.e. should have the same projection to 2D, is this correct? Also, could you suggest software (preferably free) that can visualize geography data? I tried OpenJUMP, but it doesn't seem to work with geography. – Vader Feb 12 '16 at 17:16
  • Maybe you can put these questions in different posts? It keeps the overview better and helps other people to find answers. Just add a comment here with the links to them when you did. – tilt Feb 13 '16 at 18:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.