6

I want to search objects that have at least one point inside of a rectangle. Based on the size of the rectangle, PostGIS sometimes returns the contained object, sometimes not.

create table t1 (id integer not null primary key, coords geography(linestring));
insert into t1 values (1, 'linestring(9.4 48.5, 9.2 48.6, 9.0 48.8)')

So the table contains a linestring, that consists of 3 coordinates (somewhere in the area around Stuttgart).

Now I want to search all objects, that have at least one point inside of a given rectangle. I use a polygon for this boundary. The following query returns the linestring from above:

select t.id, st_astext(t.coords) from t1 t
where ST_Intersects(
    ST_GeogFromText('POLYGON((-110.0 49.0, 10.0 49.0, 10.0 47.0, -110.0 47.0, -110.0 49.0))'),
    t.coords)

> 1, LINESTRING(9.4 48.5,9.2 48.6,9 48.8)

In the following query, I just modified the right (eastern) border of the rectangle and moved it from lon=10.0 to lon=20.0. This query does not return any result:

select t.id, st_astext(t.coords) from t1 t
where ST_Intersects(
    ST_GeogFromText('POLYGON((-110.0 49.0, 20.0 49.0, 20.0 47.0, -110.0 47.0, -110.0 49.0))'),
    t.coords)

What am I doing wrong?

Some info about my environment:

SELECT version() || ' ' || postgis_full_version();
PostgreSQL 12.3 (Debian 12.3-1.pgdg100+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 8.3.0-6) 8.3.0, 64-bit POSTGIS="3.0.1 ec2a9aa" [EXTENSION] PGSQL="120" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 5.2.0, September 15th, 2018" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1" WAGYU="0.4.3 (Internal)"
4
  • If your data is in Germany, why does your polygon span the Atlantic and most of the United States (110W longitude)?
    – Vince
    Jul 11 '20 at 15:09
  • When casting to geometry e.g select t.id, st_astext(t.coords) from t1 t where ST_Intersects(ST_GeogFromText('POLYGON((-110.0 49.0, 10.0 49.0, 10.0 47.0, -110.0 47.0, -110.0 49.0))')::geometry, t.coords::geometry) no more problem. No more clue than this
    – ThomasG77
    Jul 11 '20 at 15:37
  • Casting to geometry solves my problem, thank you. @JGH explaines below (and in the linked post) the source of the problem.
    – mre
    Jul 11 '20 at 16:14
  • BTW: The "rectangle" (polygon) is the visible area of a leaflet map. That's the reason, why the rectangle could be so large. The linestring is the geography object stored in the database. If at least one single point of the linestring is inside the visible rectangle of the leaflet, then I have to fetch it from the database in order to draw it and display its details in a tooltop/popup.
    – mre
    Jul 11 '20 at 16:39
6

Geography uses great circles instead of straight lines to connect two points. A great circle arc is the line intersecting the earth surface and a plan going through the 2 points and the center of the earth.

Consequently, a coordinate slightly north to the coordinate of one of the 2 connected point will be "above" (north of) the great circle arc near the 2 points, but will be "below" (south of) it further away.

By moving the intersection box further to the east, you also moved the distance between your line and the box edge and it is now outside of it.

This similar post contains some illustration to depict the observed behavior.

Let's quantify the difference. The 1st step is to segmentize the polygon (i.e. to add more vertices so its great-circle shape is preserved when being projected), then let's evaluate where are the intersections between a north-south line going though longitude 9.2

  1. using the original polygons:
WITH bb (geog) as (VALUES (ST_Segmentize(ST_GeogFromText('POLYGON((-110.0 49.0, 10.0 49.0, 10.0 47.0, -110.0 47.0, -110.0 49.0))'),1000))),
    blade(geog) as (VALUES (ST_GeogFromText('LINESTRING(9.2 0, 9.2 80)')))
SELECT st_asText(st_intersection(bb.geog,blade.geog))
FROM bb, blade;
                       st_astext
-------------------------------------------------------
 LINESTRING(9.2 47.6795376273658,9.2 49.6740119023496)
  1. using the extended polygon:
WITH bb (geog) as (VALUES (ST_Segmentize(ST_GeogFromText('POLYGON((-110.0 49.0, 20.0 49.0, 20.0 47.0, -110.0 47.0, -110.0 49.0))'),1000))),
    blade(geog) as (VALUES (ST_GeogFromText('LINESTRING(9.2 0, 9.2 80)')))
SELECT st_asText(st_intersection(bb.geog,blade.geog))
FROM bb, blade;
                       st_astext
-------------------------------------------------------
 LINESTRING(9.2 56.0310479315236,9.2 57.8695612625077)

we see that the intersection between the two is located much more to the north

You can find below a visual illustration

Original vs extended polygon

5
  • Thank you for this explanation (and the link). When I cast to "geometry", all my queries work as expected. select t.id, st_astext(t.coords) from t1 t where ST_Intersects(ST_GeomFromText('POLYGON((-110.0 49.0, 20.0 49.0, 20.0 47.0, -110.0 47.0, -110.0 49.0))'), t.coords::geometry)
    – mre
    Jul 11 '20 at 16:15
  • 2
    Well, it "works" but using straight lines on such large distances is not necessarily the proper way to go...
    – JGH
    Jul 11 '20 at 16:18
  • 1
    @JGH Thanks for your answer. Added a visual illustration to it as I feel the need for it (not edited any other part of your answer)
    – ThomasG77
    Jul 11 '20 at 18:42
  • Very enlightening. Thank you.
    – mre
    Jul 11 '20 at 19:19
  • 1
    @ThomasG77 thanks! a picture is indeed worth a thousand line of code, hum, I mean words
    – JGH
    Jul 12 '20 at 15:29

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.