4

I'm a newbie to POSTGIS and I wonder if someone could help me out.

I have loaded all postal codes with their geographical coordinates (geography) and created an index on this coordinate.

When I run the following query it works just fine (within a few millisecond over a few million instances)

SELECT
  postcode, huisnummer, huisletter, huisnummertoevoeging,
  ST_X(geolocatie :: geometry) as latitude,
  ST_Y(geolocatie :: geometry) as longitude,
  st_distance(geolocatie, ST_SetSRID(ST_Point(4.3387478, 51.9808117), 4326)) as afstand
FROM inspire
WHERE ST_DWithin(geolocatie, ST_SetSRID(ST_Point(4.3387478, 51.9808117), 4326), 1000)
ORDER BY afstand
LIMIT 5;

geolocatie is of type geography and I have defined the following index on geolocatie

CREATE INDEX inspire_geolocatie_index
    ON public.inspire USING gist
    (geolocatie)
    TABLESPACE pg_default;

This works fine as long as there are neighbors, but the query seems to have difficulties for the use case that there are no neighbors. Any suggestions please?

  • If you do "LIMIT 5" then it stops after finding the first five matches. That might be the first five rows of your table. It will be quick. If your condition has no matches, PostGIS has to try every row. Is that the problem? – Spacedman Sep 24 '18 at 19:39
  • @Spacedman I think in this case, PG would always execute the ORDER BY (that you can try and rewrite to geolocatie <-> ST_SetSRID(ST_Point(4.3387478, 51.9808117), 4326) for a possible boost...and maybe place the point creation into a CTE), and should not interfere with the no-neighbor case...I also cannot reproduce this; could you add your table structure in more detail and the EXPLAIN ANALYZE <your_query> output for both cases? PG/PostGIS versions might be nice, too. – ThingumaBob Sep 24 '18 at 20:19
  • 1
    btw., it is advisable to use explicit casts, if only for clarity here, e.g. ST_SetSRID(ST_Point(4.3387478, 51.9808117), 4326)::geography; most functions will have multiple signatures where a cast can make a huge difference (and many don't have an implicit cast like here). – ThingumaBob Sep 24 '18 at 20:24
  • My mistake yes, I didn't notice the "order by". it will have to compute all the distances. – Spacedman Sep 24 '18 at 22:18
  • You are so helpful. Thank you so much! For me it's a real adventure :). BTW I looked at your suggestions and using the explicit cast ::geography did the trick I now do get a very fast response without any answers. And that's exactly what I want. – user2767718 Sep 26 '18 at 21:49
2

You do have a solid query right there, and I can´t seem to reproduce your issue (as per the title), assuming a correctly created and indexed geography column.

If I were to improve things here, I´d use

SELECT
  postcode, huisnummer, huisletter, huisnummertoevoeging,
  ST_X(geolocatie::geometry) AS latitude,
  ST_Y(geolocatie::geometry) AS longitude,
  ST_Distance(geolocatie, ST_SetSRID(ST_MakePoint(4.3387478, 51.9808117), 4326)::geography) AS afstand
FROM inspire
ORDER BY geolocatie <-> ST_SetSRID(ST_MakePoint(4.3387478, 51.9808117), 4326)::geography
LIMIT 5;

The KNN operator <-> as the ORDER BY parameter is ultimate in performance when used properly; I would even forgo a limiting filter, especially as ST_DWithin, in my experience, adds overhead and denies the planner the much desired index (only) scans more often than not.


I wouldn´t recommend this without proper transformation (as, e.g. @Michael proposed), but an index forcing alternative filter option would be

SELECT
  postcode, huisnummer, huisletter, huisnummertoevoeging,
  ST_X(geolocatie::geometry) AS latitude,
  ST_Y(geolocatie::geometry) AS longitude,
  ST_Distance(geolocatie, ST_SetSRID(ST_MakePoint(4.3387478, 51.9808117), 4326)::geography) AS afstand
FROM inspire
WHERE geolocatie::geography && ST_Expand(ST_SetSRID(ST_MakePoint(4.3387478, 51.9808117), 4326), 0.01)
ORDER BY geolocatie <-> ST_SetSRID(ST_MakePoint(4.3387478, 51.9808117), 4326)::geography
LIMIT 5;

i.e. using the bbox intersection operator && directly; major drawback here is the dependency on the geometry type (and thus the CRS units) for ST_Expand, that leaves you with a dangerous 0.01 degrees (~1.1km, but only at the equator) to be used as distance parameter without a meter based projection.


Without your adding details to the question, issue tracking will be difficult; to start with, the EXPLAIN ANALYZE <your_query> for both plans is almost mandatory for query performance related questions, as well as your PostgreSQl/PostGIS versions.

  • To be sure to have an geography index: CREATE INDEX idx_inspire_geolocatie ON inspire USING GIST (cast(geolocatie as geography)); – robin loche Sep 26 '18 at 9:02
  • Thank you ThingumaBob. Your answer was very helpful. Using the explicit cast, as you suggested, solved the problem. – user2767718 Sep 26 '18 at 21:57
  • Your suggested queries also perform very well, viz.1 millisecond. Really good!. The original modified query takes 0.057 seconds using your first suggestion to use an explicit cast, which is not bad either. I still have a lot to learn.... – user2767718 Sep 26 '18 at 22:18
2

You mention your query has difficulties for the case that are no neighbours - this is due to the 'where' statement, which excludes observations beyond 1000 units of your coordinate system. In the documentation of st_dwithin it states the following concerning distance: units are in those of spatial reference.

This brings me to a second point: 1000 units in WGS84 would be very far. You should use a projected coordinate system with st_transform. For the Netherlands this would be EPSG 28992.

To know which points are within 1000 meters you can execute following query:

SELECT
  postcode, huisnummer, huisletter, huisnummertoevoeging,
  ST_X(ST_Transform(geolocatie,28992) :: geometry) as latitude, -- assuming your geometry is in 4326
  ST_Y(ST_Transform(geolocatie,28992)  :: geometry) as longitude,
  ST_Distance(
        ST_Transform(geolocatie,28992),
        ST_Transform(ST_SetSRID(ST_Point(4.3387478, 51.9808117),4326),28992)) as afstand , 
  ST_DWithin(
        ST_Transform(geolocatie,28992) , 
        ST_Transform(ST_SetSRID(ST_Point(4.3387478, 51.9808117),4326),28992),1000) as within_distance 
FROM inspire 
ORDER BY afstand

Performance wise, this can be improved by adding a geometry to the table in advance with the project geometry and adding an index.

-- Add the column
ALTER TABLE inspire ADD COLUMN geom_28992 geometry;
-- Set the column
UPDATE inspire SET geom_28992 = st_transform(geom,28992);
-- Create an index
CREATE INDEX idx_geom_28992_inspire ON inspire USING GIST(geom_28992);

For st_dwithin to work efficiently both geometries need to be indexed: the index of reference geometry needs to be in a table as well, with index.

-- Reference geom
CREATE TABLE reference_geom AS
SELECT ST_Transform(ST_SetSRID(ST_Point(4.3387478, 51.9808117),4326),28992) geom_28992;
-- Create an index
CREATE INDEX idx_geom_28992_reference_geom ON reference_geom USING GIST(geom_28992);

Finally, you can execute the query

SELECT postcode, huisnummer, huisletter, huisnummertoevoeging, 
st_distance(i.geom_28992,g.geom_28992) distance,
true::boolean within_distance
FROM inspire i
INNER JOIN reference_geom  g ON st_dwithin(i.geom_28992,g.geom_28992,1000)
UNION ALL
SELECT postcode, huisnummer, huisletter, huisnummertoevoeging, 
null AS distance,
false::boolean within_distance
FROM inspire i,reference_geom
WHERE st_dwithin(i.geom_28992,g.geom_28992,1000) = FALSE

For more interesting applications have a look at nearest neighbour problems in Postgis, for instance here and here.

  • there's hidden advice here, but some misconceptions, too; 1st: for point-on-point proximity searches, an index is not necessarily beneficial (but seeing to it being in place doesn´t hurt either), since bboxes of points are in fact the points themselves; 2nd: for each single geometry, only the index of the filtered table will make a difference in performance. 3rd: using geography type will, in fact, implicitly use meter as units when in use with appropriate functions, and ST_Distance and ST_DWithin will both implicitly cast the second input geometry to the type of the first, if possible – ThingumaBob Sep 25 '18 at 10:14
  • Thanks, it makes senses that for point-on-point proximity searches the index is not beneficial; I simply had a look at the documentation. Also, I was not aware that using ST_Distance and ST_DWithin with the geography use meters as units, but the documentation indeed states it clearly. – Michael Sep 26 '18 at 2:04

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.