Thursday 28 September 2017

postgis - How do I get the closest point on a road to a point when I have imported OSM data using osm2pgsql?


I have been trying to use OSM data to find the closest point to a road for a given point. I imported the OSM data into postgis using osm2pgsql, and I am getting results back, but they are at the absolute edge of my imported data (and the point I'm trying to specify is near the middle) and distance reported is about 8000km. I'm pretty sure it has to be something to do with mixing different coordinate systems but I'm quite new to this so not sure.


The SQL query I'm trying to use is:


SELECT osm_id,

highway,
name,
ref,
ST_X(ST_ClosestPoint(ST_Transform(r.way,4326),
ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326))),
ST_Y(ST_ClosestPoint(ST_Transform(r.way,4326),
ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326))),
ST_Distance_Sphere(ST_ClosestPoint(ST_Transform(r.way,4326),
ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326)),
ST_SetSRID(ST_MakePoint(55.938627,-3.198706),4326))

FROM planet_osm_roads r
ORDER BY 7 ASC
LIMIT 10;

Answer



First of all, ST_Distance_Sphere returns in meters, so you are actually looking at 8km, which might be more reasonable. I suspect, also, that you have you lat/lon the wrong way round -- your point is somewhere off the coast of the Seychelles, not Carlisle, and there are not too many roads in the Indian Ocean.


Also, while the query planner will optimize this away, there is no need to repeatedly enter the same point. Instead, create it in a sub-query and reference it, eg:


SELECT 
osm_id,
highway,
name,

ref,
ST_X(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom),
ST_Y(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom),
ST_Distance_Sphere(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom), point.geom)
FROM
planet_osm_roads r,
(Select ST_SetSRID(ST_MakePoint(-3.198706, 55.938627), 4326) as geom) point
ORDER BY 7 ASC
LIMIT 10;


However, this is still very inefficient. To make better use of a indexes, look at using ST_DWithin (r.way, point, distance) ORDER BY ST_Distance(r.way, point). If you look at the docs for ST_DWithin you will get an idea how it is used.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...