2

I am using ST_Distance in postgis to get the closest point in a table to each polygon in another table. This query is taking a long time to run (2 hours with 18,000 records in the polygon table and 7,500 records in the point table). I am calculating the distance in miles and also limiting the search to points that are within 100 miles (160934 meters) of each polygon. What else could I do to speed this up? The geometries are all in SRID 4326, and I have created spatial indexes on both tables. I'm running postgresql 9.5 with postgis 2.2.2. Here is the query:

SELECT DISTINCT ON(g1.gid) g1.*, g2.plant_name As closest_utility_name,
g2.primsource As closest_utility_type,
round(ST_Distance(g1.geom::geography,g2.geom::geography) * 0.000621371) As closest_utility_miles
INTO parcels_utility
FROM parcels As g1, powerplants_us_201603 As g2
WHERE ST_DWithin(g1.geom, g2.geom, 160934)
ORDER BY g1.gid, closest_utility_miles;

Explain provides this query plan:

"Unique  (cost=11006.38..11006.60 rows=45 width=713)"
"  ->  Sort  (cost=11006.38..11006.49 rows=45 width=713)"
"        Sort Key: g1.gid, (round((_st_distance((g1.geom)::geography, (g2.geom)::geography, '0'::double precision, true) * '0.000621371'::double precision)))"
"        ->  Nested Loop  (cost=0.15..11005.14 rows=45 width=713)"
"              ->  Seq Scan on parcels g1  (cost=0.00..1692.01 rows=18001 width=653)"
"              ->  Index Scan using powerplants_us_201603_gix on powerplants_us_201603 g2  (cost=0.15..0.51 rows=1 width=60)"
"                    Index Cond: (geom && st_expand(g1.geom, '160934'::double precision))"
"                    Filter: ((g1.geom && st_expand(geom, '160934'::double precision)) AND _st_dwithin(g1.geom, geom, '160934'::double precision))"

Based on the comment below I've updated the ST_DWithin function to use the geography type:

SELECT DISTINCT ON(g1.gid) g1.*, g2.plant_name As closest_utility_name,
g2.primsource As closest_utility_type,
round(ST_Distance(g1.geom::geography,g2.geom::geography) * 0.000621371) As closest_utility_miles
INTO parcels_utility
FROM parcels As g1, powerplants_us_201603 As g2
WHERE ST_DWithin(g1.geom::geography, g2.geom::geography, 160934)
ORDER BY g1.gid, closest_utility_miles;

It has been running now for about an hour, the query plan is below:

"Unique  (cost=612038648.76..612038650.68 rows=200 width=455)"
"  ->  Sort  (cost=612038648.76..612038649.72 rows=383 width=455)"
"        Sort Key: g1.gid, (round((_st_distance((g1.geom)::geography, (g2.geom)::geography, '0'::double precision, true) * '0.000621371'::double precision)))"
"        ->  Nested Loop  (cost=0.00..612038632.33 rows=383 width=455)"
"              Join Filter: (((g1.geom)::geography && _st_expand((g2.geom)::geography, '160934'::double precision)) AND ((g2.geom)::geography && _st_expand((g1.geom)::geography, '160934'::double precision)) AND _st_dwithin((g1.geom)::geography, (g2.geom)::g (...)"
"              ->  Seq Scan on parcels g1  (cost=0.00..20630.55 rows=152755 width=395)"
"              ->  Materialize  (cost=0.00..430.86 rows=7524 width=60)"
"                    ->  Seq Scan on powerplants_us_201603 g2  (cost=0.00..393.24 rows=7524 width=60)"
  • Please edit the question to include necessary details: What indexes have you created on the tables? What SRID(s) are used in the geometries? What query plan does EXPLAIN return? What exact versions of PosgreSQL and PostGIS are in use? – Vince Oct 7 '16 at 0:48
  • Thanks for your comment, that is important information I left out - I've edited my post. – kflaw Oct 7 '16 at 1:05
  • 1
    Are you sure the DWithin isn't looking for shapes within 161k degrees? – Vince Oct 7 '16 at 1:14
  • 1
    The cost went from 11k to 612m, which means the cast to geography made it impossible to use the spatial index (bad juju) – Vince Oct 7 '16 at 3:03
  • 1
    Is your data really global? If not, it would pay off (very much) to store your data in a local projection like UTM. – tilt Oct 7 '16 at 10:39
3

You're mixing units in your original attempt, and improperly using geography in your second. I'll provide a geography solution, though if your data are local the very highest performance will be using planar (not 4326) geometry.

-- Change the data type in the tables
ALTER TABLE powerplants_us_201603 ALTER COLUMN geom TYPE geography USING geography(geom);
ALTER TABLE parcels ALTER COLUMN geom TYPE geography USING geography(geom);
ALTER TABLE powerplants_us_201603 RENAME COLUMN geom TO geog;
ALTER TABLE parcels RENAME COLUMN geom TO geog;

-- Build real geography indexes
CREATE INDEX parcels_geog_x ON parcels USING GIST (geog);
CREATE INDEX powerplants_geog_x ON powerplants_us_201603 USING GIST (geog);

-- Run the (simpler now) query
SELECT DISTINCT ON (g1.gid)
  g1.*, 
  g2.plant_name AS closest_utility_name,
  g2.primsource AS closest_utility_type,
  round(ST_Distance(g1.geog, g2.geog) / 1609.34) AS closest_utility_miles
FROM parcels AS g1, 
JOIN powerplants_us_201603 AS g2
ON ST_DWithin(g1.geog, g2.geog, 100 * 1609.34)
ORDER BY g1.gid, closest_utility_miles;
  • Thanks for your helpful comments and answer. I'm learning something new all the time! – kflaw Oct 18 '16 at 2:48
0

We don't have the dataset available, but I'll give it a shot. Generally, ST_Distance is slow, because it needs to do the math each time without using any spatial index. And it needs to do it in 100 miles radius, which is quite a lot (I guess).

How about this LATERAL JOIN based query? You could experiment with geometry/ST_Transform to see which works better, I actually have very little experience with geography type myself.

SELECT
g1.*,
g2.closest_utility_name,
g2.closest_utility_type
FROM parcels g1
JOIN LATERAL (
SELECT
    g2.plant_name closest_utility_name,
    g2.primsource closest_utility_type
FROM powerplants_us_201603 g2
ORDER BY
    g1.geom::geography <-> g2.geom::geography
LIMIT 1
) g2 ON true;
  • 1
    This promiscuous casting to geography is only going to harm you, since if the original data is geometry, it's very likely there is not geography index, and the <-> operator is only going to give you speed if it has a native index on it. – Paul Ramsey Oct 11 '16 at 22:24

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.