[GIS] Find out whether geography point lies within geography bounds

postgispostgis-2.0

I have three points:

Point A
Point SW
Point NE

I want to know if Point A resides within the bounds determined by points SW (for south-west coordinates) and NE (for north-east coordinates).

I've so far tried both ST_Intersects and && to determine this, which works fine for small regions, but fails whenever the bounding area determined by SW and NE points covers more than half the globe. In this case it seems to instead use a different, smaller, bounding area using the same four points.

As an example, if I'd want to cover everything from the SW point of South Africa to the NE point of Ireland it would instead catch points where the the NW point is in Ireland and the the SE point is in Madagascar. That is, China would not considered to be within the bounds.

Here are some SQL queries to illustrate this:

select 
ST_GeographyFromText('SRID=4326;POINT(0 0)') && ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  -80 -80, -80 80, 180 80, 180 -80))') as expectTrue_1, 
ST_Intersects(ST_GeographyFromText('SRID=4326;POINT(0 0)'), ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  -80 -80, -80 80, 180 80, 180 -80))')) as expectTrue_2, 
ST_GeographyFromText('SRID=4326;POINT(0 0)') && ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  180 80, -80 80, -80 -80, 180 -80))') as expectTrue_3, 
ST_Intersects(ST_GeographyFromText('SRID=4326;POINT(0 0)'), ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  180 80, -80 80, -80 -80, 180 -80))')) as expectTrue_5, 
ST_GeographyFromText('SRID=4326;POINT(-90 0)') && ST_GeographyFromText('SRID=4326;POLYGON((180 -80,  -80 -80, -80 80, 180 80, 180 -80))') as expectFalse_1, 
ST_GeographyFromText('SRID=4326;POINT(0 0)') && ST_GeographyFromText('SRID=4326;POLYGON((1 -1,  -1 -1, -1 1, 1 1, 1 -1))') as expectTrue_5;

I'd want for this to return t t t t f t but it returns f f f f t t.

I've also tried ST_Intersects(ST_MakeBox2D()) as per the following SO post:
https://stackoverflow.com/questions/24212355/query-by-coordinates-takes-too-long-options-to-optimize/24218861?iemail=1&noredirect=1#24218861

I'm using PostGIS 2.1.3 with the Geography data type.

What functions or approaches should I use to have PostGIS treat my more-than-half-globe-spanning polygons as more-than-half-globe-spanning ones instead of 'optimizing' them to be less?

I understand I could just split the bounds in, say, two parts, and check for each, but there must be a cleaner solution.

Best Answer

I'm guessing that you probably figure the connectors between your points are "straight lines" and would run parallel to the edges of a mercator map of the world. If you are using the geography type, that won't happen, the edges will be great circles and will run in directions you'd never guess.

http://boundlessgeo.com/2012/04/the-earth-is-not-flat-volume-2/

In any event, no edge in geography can be longer than 180d, since any two points define a great circle, the system has to decide what portion of the circle the points bound, and chooses the shortest one (the one less than 180d). If you want longer edges you have to add intermediate points that fall in the places you'd expect.

But remember, you aren't going to get "straight lines" out of this, you're going to get great circles.

If you want a straight line solution, work in geometry, and recognize that you're going to have to put in your own logic to handle dateline splitting.

(Dateline splitting is not hard to do, if your input box is not defined as "two points" but rather as a "lower left point" and "upper right point". Then you can easily test for cases in which the "lower left" point appears to be to the right of the "upper right point" and say "aha! this is actually a dateline crossing case!" and construct two query polygons with the dateline as the bounding edge.)

Related Solutions

Python – Minimum Bounding Circle of Geometry Crossing the 180th Meridian

This is how I solved the problem:

import psycopg2
import re

def replace(matchObj):
    # Group 0 is all groups, group 1 is the first match contained in
    # parentheses. Group 2 is the whole float value.
    value = float(matchObj.group(2))
    value = value - 360
    return matchObj.group(1) + str(value) + " "

def main():
    # connect to postgresql
    conn = psycopg2.connect("conn info")
    cur = conn.cursor()
    cur.execute("select st_astext(the_geom) from tl_2010_us_state10 \
                WHERE name10='Alaska'")

    # Get query results
    rows = cur.fetchall()
    alaskaText = rows[0][0]

    # The regular expression finds a positive float between 170 and 179.9, 
    # inclusive. Matches any character not a negative sign to ensure that the
    # longitude value is positive, then it matches the rest of the number.  
    # To ensure that the number matched is longitude, we check that it is
    # followed by a space (the format being "longitude latitude,".
    matchText = re.compile('([^-])(1[67][0-9]\.[0-9]*) ')

    # Substitute longitude values
    alaskaText = matchText.sub(replace, alaskaText)

    # insert into dabase and commit transaction
    cur.execute("insert into tl_2010_us_state10 (the_geom, name10) \
                VALUES (st_geomfromtext('%s', 4269), 'Alaska2')" %(alaskaText))
    conn.commit()
    conn.close()

if __name__ == "__main__":
    main()

[GIS] Problem PostGIS Spatial Query ST_Contains ST_Disjoint ST_Intersects

You have a set of polygon, they cover an area of the map.

You have a set of points, some are inside the polygons and some are outside.

The st_intersect will compare all points against all polygons and only return points that are "on top" one or more polygons.

The st_disjoint will return a point if a point is not "on top" of a polygon. So for a set of points against a set of polygons, so unless the polygons overlap, each point is going to be "outside" of one or more polygons.

You can have a nested query of the intersect, and another inverted nested query of the intersect to get your t/f result

Adding a solution

There's a fiddle that goes with this here, since it was a PostgreSQL fiddle rather than PostGIS I've change the fiddle test to be a simple integer equality test.

So the first half of the sql is your intersect test:

    SELECT distinct points.id
    FROM points, multipolygon
    WHERE st_intersects(points.geom, multipolygon.geom)

I add a TRUE column to match your sample, by adding a wrapping SQL select with a hardcoded value of TRUE:

SELECT
   id,
   TRUE AS intersects_at_least_one_polygon
FROM 
(
    SELECT distinct points.id
    FROM points, multipolygon
    WHERE st_intersects(points.geom, multipolygon.geom)
) intersecting ;

To find the points that are NOT intersecting the polygons, well that is the total list of point id, minus the above points. So I select all the points, the outer select below, and then say "WHERE id NOT IN" and just repeat your SQL, plus I add a hardcode a FALSE:

SELECT
   points.id,
   FALSE AS intersects_at_least_one_polygon
FROM
   points WHERE id NOT IN
   (
        SELECT distinct points.id
        FROM points, multipolygon
        WHERE st_intersects(points.geom, multipolygon.geom)
   )
   order by id
   ;


SELECT
   id,
   TRUE AS intersects_at_least_one_polygon
FROM 
(
    SELECT distinct points.id
    FROM points, multipolygon
    WHERE st_intersects(points.geom, multipolygon.geom)
) intersecting
UNION
SELECT
   points.id,
   FALSE AS intersects_at_least_one_polygon
FROM
   points WHERE id NOT IN
   (
        SELECT distinct points.id
        FROM points, multipolygon
        WHERE st_intersects(points.geom, multipolygon.geom)
   )
   order by id
   ;

I'm certain that there are more elegant, efficient, solutions to this, but hope you're happy with a simple working version.

Best Answer

Related Solutions

Python – Minimum Bounding Circle of Geometry Crossing the 180th Meridian

[GIS] Problem PostGIS Spatial Query ST_Contains ST_Disjoint ST_Intersects

Related Question