2

I have a postgis point table with a large list of cities.

I need to find duplicates in it, but only if they are 10km close each other, which means that is the same city.

For example, I have 5 points called 'Eneby', and suppose that 4 of this are 10Km close each others, and the 5th is 300Km away, what will mean that is another Eneby town. The query in this case should only have to show 4 points.

So, I've made this query, but is not showing the expected result:

Select p1.id, p2.id, * from points as p1, points as p2 where p1.name='Eneby'
and st_dwithin(p1.geom,p2.geom,0.1, true)
and p1.id=p2.id

The result is the same 5 points, no matter what distance (in degrees) I set. I'm trying with different distance values, but is always the same result.

Does anybody know if is there any mistake in this query?

I've been looking through this forum, and found some similar questions, but can't find exactly this problem.

3

Assuming that id is the primary key for your points table, and that you want to group point ids which are 10km apart in order to select a best point from each collection. Note that your distance parameter in this case should be the assumed diameter of your city, not the radius. This code will first, for every point, create an ordered array of points within the diameter sharing a same city name.

WITH point_arrays AS(
SELECT p1.name, ARRAY_AGG(p2.id) ORDER BY p2.id
FROM points p1, points p2
WHERE p1.name = p2.name 
AND st_dwithin(p1.geom,p2.geom,0.1, true)
AND p1.name = 'NAME' --omit if you want to do this for all city points
GROUP BY p1.name)

Then, because of redundancy in the first block, find the distinct array of point ids for every city name.

SELECT DISTINCT name, point_array
FROM point_arrays
2

It might be a solution :

1 create buffers radius 5 for each point, keep attributes ; 2 select buffers intersecting each other with same name ; 3 keep the "best" cities ;

What do you think ?

5
  • Good idea, but still don't know how to intersect two polygons from the same layer – Pablo Pardo Sep 16 '14 at 10:50
  • You can add the same layer 2 times : then you can intersect one with the other. I can't try at this moment but if I remember well, you can do that. Tell me if it works. – Leehan Sep 16 '14 at 18:20
  • If I intersect the same layer, I'll always get the same layer. – Pablo Pardo Sep 16 '14 at 19:44
  • Sorry, I will try tomorrow morning at my office. – Leehan Sep 16 '14 at 20:02
  • Can you spare sample of your datas ? – Leehan Sep 17 '14 at 12:01
2

I think your p1.id=p2.id is the problem. If you change it to p1.id<>p2.id and add a p1.name='Eneby, I think you'll get what you are after.

with cities as 
(
select 1 ID, 'Eneby' as CityName, ST_MakePoint(10, 10) geom
union
select 2 ID, 'Eneby' as CityName, ST_MakePoint(11, 11) geom
union
select 3 ID, 'Eneby' as CityName, ST_MakePoint(11, 12) geom
union 
select 4 ID, 'City B' as CityName, ST_MakePoint(50, 50) geom
union 
select 5 ID, 'City B' as CityName, ST_MakePoint(51, 52) geom
union
select 6 ID, 'Eneby' as CityName, ST_MakePoint(500, 500) geom
)
select ST_DWithin(p1.geom,p2.geom,300) iswithin, ST_Distance(p1.geom,p2.geom) dist,  p1.ID,p2.ID ID2,   p1.CityName
    from cities p1, cities p2
    where p1.ID<>p2.ID
    and p1.CityName='Eneby' and p2.CityName='Eneby'
    and ST_DWithin(p1.geom,p2.geom,300)

And if you wanted to do all the cities, I think you could get all the distances between the points and then group the city names by the name and whether they are with x distance of one another:

with cities as 
(
select 1 ID, 'City A' as CityName, ST_MakePoint(10, 10) geom
union
select 2 ID, 'City A' as CityName, ST_MakePoint(11, 11) geom
union
select 3 ID, 'City A' as CityName, ST_MakePoint(11, 12) geom
union 
select 4 ID, 'City B' as CityName, ST_MakePoint(50, 50) geom
union 
select 5 ID, 'City B' as CityName, ST_MakePoint(51, 52) geom
union
select 6 ID, 'City A' as CityName, ST_MakePoint(500, 500) geom
)
select a.CityName, within300,  ST_AsText(ST_Collect(geom2)) as geomcoll
from 
(
    select c1.geom, c2.geom geom2 ,c1.ID, c2.ID ID2, c1.CityName, c2.CityName CityName2, ST_Distance( c1.geom,c2.geom), CASE WHEN ST_Distance( c1.geom,c2.geom) <  300 THEN 1 ELSE 0 END  within300
    from cities c1, cities c2
    where c1.ID <> c2.ID 
) a
where CityName = CityName2
group by a.CityName, within300
1
  • 1
    To avoid doing the distance both ways you could use c1.ID < c2.ID – MickyT Sep 15 '14 at 20:58

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.