ST_Distance does not give nearest node

Question

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. . 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

Answer 1

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.