[postgis-users] ST_DWithin with a rectangular shape

zwieberl at lavabit.com zwieberl at lavabit.com
Wed Dec 19 04:00:13 PST 2012


Hello,

I'm using Postgis 2.0.2 with Geography-Points only
(satellite-measurements). I have two tables for two different satellites
and need to find all Points from table A that are in the vicinity of the
points in table B.

Now ST_DWithin does this job amazingly well (200 s for 14.216.117 hits),
but since earth is rotating I need a asymmetrical search-area (bigger
longitudes than latitudes).
Basically I need an ST_DWithin-function for rectangles, where I can give
two distances (one in a direction parallel to the equator and one
perpendicular to it) around a given point.

Right now I was trying to solve the problem the following way:

I create a POLYGON on-the-fly during the SELECT-statement out of points
from table B using ST_Project and use than ST_Covers() on the created
Polygon:

select A.id, B.id from A, B WHERE ST_Covers(
ST_GeogFromText(
'POLYGON(('
|| trim(trailing ')' from trim(leading 'POINT(' from
ST_AsText(ST_Project(B.geo, 1500000, 0)))) || ','
|| trim(trailing ')' from trim(leading 'POINT(' from
ST_AsText(ST_Project(B.geo, 1500000, 1.57)))) || ','
|| trim(trailing ')' from trim(leading 'POINT(' from
ST_AsText(ST_Project(B.geo, 1500000, 3.1415)))) || ','
|| trim(trailing ')' from trim(leading 'POINT(' from
ST_AsText(ST_Project(B.geo, 1500000, 4.71)))) || ','
|| trim(trailing ')' from trim(leading 'POINT(' from
ST_AsText(ST_Project(B.geo, 1500000, 0)))) ||
'))'
), A.geo);

Note 1: This Polygon is wrong! It is turned by 45°. I would need to do a
ST_Project twice for every point, but for shortness I omitted it here. (I
would need to do something like ST_Project(ST_Project(B.geo, 1500000,
1.57), 800000, 0) for every corner of the rectangle).

Note 2: I know I can use ST_X(point::geometry) to extract the coordinates,
but its basically the same with respect to performance, especially since I
need to project the point then four times instead of two (twice for ST_X
and twice for ST_Y)

This solution is bad. Really, really bad! As stated above ST_DWithin does
this in 200 s for 14 million hits. I had to cancel my solution after 4
hours for the same tables.

Does anybody now a solution to this problem?

Any help would be really appreciated!

Thank you!






More information about the postgis-users mailing list