2

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.

2
  • you need to make sure your tables are in the same SRID -- and do not use 4326 for distance calcualtions. try converting to more local SRID
    – ziggy
    Oct 7 '19 at 20:47
  • Yeah I just noticed that. The table from the OSM file is 3857, so I converted my coordinates table to 3857 as well. But I'm getting 0 rows returned.
    – foadster
    Oct 7 '19 at 21:32
3

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.

12
  • Thank you, I'll give it a try. I don't see a way to import my OSM file with osm2pgsql and have it use the raw geographic coordinates (EPSG:4326). Anyway to do this?
    – foadster
    Oct 8 '19 at 14:51
  • I tried your query, no luck. I get the error: "ERROR: column b.geom does not exist". I know I need to use the highway type within planet_osm_line. See the FROM clause in my original query "planet_osm_line highway"
    – foadster
    Oct 8 '19 at 16:30
  • osm2pgsql has the -l switch that preserves the Lat/Lon. alternatively, you can alter the column type and transform, or transform on-the-fly like in the example with the other geom column. you need to replace geom with the actual name of the geometry column...geom is very common, but the OSM tables have way I think?
    – geozelot
    Oct 8 '19 at 16:45
  • Yeah you're right, my bad, didn't catch the geom/way in your query. Fixed it and it works, but still no results. I may reload the OSM file using the -l switch to see if that helps. Unsure why I'm not getting any results.
    – foadster
    Oct 8 '19 at 16:48
  • 1) are you sure those points are within 5m of any line? Note: planet_osm_lines are mainly admin borders and other things, not roads; just in case roads is what you want to check against, use planet_osm_roads 3) could you post the output from SELECT ST_AsEWKT(geom) FROM gps_coordinates LIMIT 1; and the same with way and the OSM table
    – geozelot
    Oct 8 '19 at 16:54
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.

6
  • If I used a value of 5 for example, wouldn't it return everything within a 500 km radius? I'm getting 0 rows returned, and in theory I would get some sort of result.
    – foadster
    Oct 7 '19 at 21:33
  • Correct. The distance is just one of your problems. How are you setting your crs? You need to 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 and if you use a projection with units of meters you can do your distance calcs in meters rather than degrees. Does that help?
    – jbalk
    Oct 7 '19 at 22:46
  • Yeah that helps. So after I imported my OSM file, I checked the SRID type by doing a SELECT ST_SRID on the geom column (highway.way). It returned 3857. So I had to drop my custom coordinates table and redefine the geom column as 3857. Even after doing this, I still receive no results. To answer your statements, my CRS for both layers match now correct? So my calculations should work
    – foadster
    Oct 8 '19 at 7:02
  • This may be a silly question, but have you plotted your points and ways on a map? Do they visually line up how you would expect? It may be that they don't line up because of an error somewhere in the data?
    – jbalk
    Oct 8 '19 at 17:18
  • Also from your comment I am not sure that you've defined your geom column correctly. What is the native projection for the coordinates table? You need to define the geom as the native projection first. So if it is 4326, define the column as 4326, then when you use that geometry, use it with st_transform like: st_transform(coordinates.geom,3857). You cannot just arbitrarily assign a projection to the data and expect it to transform itself automagically.
    – jbalk
    Oct 8 '19 at 17:22

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.