How to get distance between a point and a list of highways?

Question

I'm new to GIS. I'm trying to determine whether a list of GPS points are within a highway/road (in terms of meters). I have downloaded a OSM file and populated my database using osm2pgsql.

I have the following tables:

                     List of relations
 Schema |            Name            |   Type   |  Owner   
--------+----------------------------+----------+----------
 public | geography_columns          | view     | map_user
 public | geometry_columns           | view     | map_user
 public | planet_osm_line            | table    | map_user
 public | planet_osm_point           | table    | map_user
 public | planet_osm_polygon         | table    | map_user
 public | planet_osm_roads           | table    | map_user
 public | raster_columns             | view     | map_user
 public | raster_overviews           | view     | map_user
 public | spatial_ref_sys            | table    | map_user
 public | gps_coordinates            | table    | map_user
 public | gps_coordinates_id_seq     | sequence | map_user

My coordinates table looks like so:

id | longitude  | latitude  |                        geom                        

I populated my table with a few longitude/latitude points, then ran the following function to generate a geometry value:

UPDATE vehicle_coordinates SET geom = ST_SetSRID(ST_MakePoint(longitude, latitude), 4326);

I initially tried using the PostGIS function ST_DWITHIN like so:

SELECT gps_coordinates.id, highway.name
FROM gps_coordinates, planet_osm_line highway
WHERE highway.highway IS NOT NULL
AND ST_DWITHIN(gps_coordinates.geom, highway.way, 5);

But no results were returned. I don't think I can just use meters in this function as I'm using latitude/longitude as my base unit. I'm stuck and not sure where to go from here. Can I use St_Distance in my WHERE clause? It states that it returns the distance in meters. But I'm not getting any rows returned:

SELECT gps_coordinates.id, highway.name
FROM gps_coordinates, planet_osm_line highway
WHERE highway.highway IS NOT NULL
ST_DistanceSphere(vehicle_coordinates.geom, highway.way) <= 5;

It could be that my data is bad (I got it from GEOFABRIK, so should be good) or that I'm confusing my units here.

Answer 1

By default, osm2pgsql transforms the raw geographic coordinates (EPSG:4326 - WGS84) to EPSG:3857; to populate the geom column in your gps_coordinates table accordingly (i.e. defined as GEOMETRY(POINT, 3857)), run

UPDATE gps_coordinates
  SET geom = ST_Transform(ST_SetSRID(ST_MakePoint(longitude, latitude), 4326), 3857)
;

After that, running

SELECT a.id,
       b.name
FROM   gps_coordinates AS a
JOIN   planet_osm_line AS b
  ON   ST_DWithin(a.geom, b.way, 5)
WHERE  b.highway IS NOT NULL
;

will return all b.name of lines per a.id that are within 5 meters of each a.ids point (meaning that, if there are 3 lines near a point, you'll get three rows with the same a.id and different b.name).


Alternatively, you can transform on-the-fly; with the geom column of gps_coordinates defined as GEOMETRY(POINT, 4326), and populated with

UPDATE gps_coordinates
  SET geom = ST_SetSRID(ST_MakePoint(longitude, latitude), 4326)
;

running

SELECT a.id,
       b.name
FROM   gps_coordinates AS a
JOIN   planet_osm_line AS b
  ON   ST_DWithin(ST_Transform(a.geom, 3857), b.way, 5)
WHERE  b.highway IS NOT NULL
;

will yield the same result.

---

However:

EPSG:3857 is pretty useless for anything else than serving map (image) tiles, especially doing anything related to measurements, and you should avoid it if possible; in fact, if you don't have a specific reason to store geometries in any (projected) CRS, store them in a geographic coordinate system (i.e. EPSG:4326 - WGS84 Lat/Lon); PostGIS is really good with those (and you can transform on-the-fly), but as you suspected correctly, the unit degree that comes with them is equally useless for measurements.

While Projections have their limits outside of their defined regions (distortions), PostGIS has the GEOGRAPHY type along the standard GEOMETRY type you are currently using that is specifically made to work with spheroidal geographic coordinates, and e.g. measurements are as precise as they can get (at a slight cost in processing time). This is a wonderful alternative to finding a suitable projection for your (possibly multiple) ROI, and you can cast between types on-the-fly (e.g. way::GEOGRAPHY) to access their functionality.

I suggest to get a good read on PostgreSQL/PostGIS (spatial) index behaviour and the GEOGRAPHY type.

Answer 2

The documentation for ST_Distance states: For geometry type returns the minimum 2D Cartesian distance between two geometries in projected units (spatial ref units).

Therefore, it is not returning meters like you assumed, it is returning degrees. Depending on your latitude, one degree is equal to around 111 KM!

https://en.wikipedia.org/wiki/Decimal_degrees

You can use a rough estimate of 5m = 0.00005 degrees.

Use ST_Transform to transform between crs's. ST_SetSRID only defines the crs.

The crs for both layers must match or you will not be able to do any spatial analysis between the two.

If you transform into a projected coordinate system it will improve your distance calculations.

If you use a projected coordinate system with units of meters you can do your distance calcs in meters rather than degrees.