2

I am trying to find the nearest point on the road network closest to a given point. I imported the road network into PostGIS from QGIS and am using the SQL window in QGIS DB Manager to query the road table. The query I used is:

SELECT * 
FROM public."wa_transformed"
ORDER BY ST_Distance(
ST_StartPoint(geom),
st_SetSrid(ST_MakePoint(-122.2652671, 47.30995661), 4326),
true
) ASC
LIMIT 1;

The idea is to find the distance between all the roads in the road network and find the first item in the sorted list. Image with road, point and result. The red marker is the point (47.309, -122.265) whereas the thick brown line on the left is the query result (which is definitely not the closest source from red point).

I had to transform the roads (original SRID - 4326) to match the OpenStreetMap layer on QGIS using a CRS defined by (+proj=lcc +lat_1=45.83333333333334 +lat_2=47.33333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0000000001 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs). However, I get the same result whether I use the transformed table or original table.

I am using POSTGIS="2.5.2 r17328" [EXTENSION] PGSQL="110" GEOS="3.7.0-CAPI-1.11.0 3.7.1" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.4, released 2018/03/19" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER

The original roads shapefile is here and the transformed shapefile is here

  • 2
    "I had to transform the roads to match the OpenStreetMap layer on QGIS". Something is wrong here, please open another question where you describe the issue you had there and what exactly you did and why. – bugmenot123 Aug 15 at 8:10
  • I am not certain of all the transformations that were applied to the shapefile prior to me getting it. The main purpose of those transformations was to make the network routable for GAMA by joining disconnected roads etc. I am trying to find that and will post it separately. The projection string in the post above is just the QGIS default CRS USER 10000 which seemed to work. – Chintan Pathak Aug 15 at 19:12
3

The actual issue here is the geometry type: your data contains/consists of MultiLineStrings, and those don't have a ST_StartPoint (because it's ambiguous if the geometries were multi part); your original shapefile, however, contains only single part MultiLineStrings.


Your original data is projected in an unknown CRS (you might want to recheck with your data provider, though) extremely close to EPSG:2927, and thus to be considered equal to EPSG:4326 for most use-cases.

However, to use it on a point in EPSG:4326, you need to transform either of them to match the other; here are some options depending on what you did to your data, including the extraction of geometry into plain LineStrings:

  • import your original with target CRS EPSG:2926 and update the geometry column to EPSG:4326:

    ALTER TABLE wa_roads
      ALTER COLUMN geom TYPE GEOMETRY(LINESTRING, 4326)
        USING ST_Transform(ST_GeometryN(geom, 1), 4326)
    ;
    
  • update your wa_roads table (the original) if you have set the wrong CRS (not tried to transform):

    ALTER TABLE wa_roads
      ALTER COLUMN geom TYPE GEOMETRY(LINESTRING, 4326)
        USING ST_Transform(ST_SetSRID(ST_GeometryN(geom, 1), 2927), 4326)
    ;
    
  • reimport if you have applied any transformation (see first bullet point), since it is just easier than to try to backtrace what has happened

  • you can also create a new table using the geometry manipulation sequences above

Make sure you have a spatial index in place; if none has been created, run:

CREATE INDEX sidx_geom__wa_roads  --name it as you like; I use [s]idx_<column>__<table>
  ON wa_roads
  USING GIST (geom)
;

Then run the Nearest Neighbor search using the more optimized <->operator:

SELECT *
FROM   wa_roads
ORDER BY
       geom <-> ST_SetSRID(ST_MakePoint(-122.2652671, 47.30995661), 4326)
LIMIT  1
;

to find the overall closest geometry, or

SELECT *
FROM   wa_roads
ORDER BY
       ST_StartPoint(geom) <-> ST_SetSRID(ST_MakePoint(-122.2652671, 47.30995661), 4326)
LIMIT  1
;

to get the geometry with the closest ST_StartPoint. The answer then is wa_roads."ID" = 1167.

Note that, if the latter is what you need, consider using an index on ST_StartPoint(geom) to power the <-> operator.

  • I had to apply the second ALTER command for it to work. The new WA_roads layer also doesnt need any transformation to match with the OSM layer. The nearest point is quite accurate too. @ThingumaBob : Would you happen to know what the <-> operator is called? I had seen it before but could not find any documentation for it. – Chintan Pathak Aug 15 at 19:08
  • @ChintanPathak PostGIS operators are nameless, strictly speaking, but it's commonly referred to as KNN operator. – ThingumaBob Aug 15 at 19:15

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.