[GIS] Searching planet_osm_point by longitude and latitude

latitude longitudeosm2pgsqlpostgispostgresqlquery

I used osm2pgsql to import data into my PostgreSQL database and it put gave me the following 4 tables, planet_osm_point, planet_osm_line, planet_osm_polygon and planet_osm_roads.

I have looked at the schema for osm2pgsql and inside the planet_osm_point table there is a column 'way' of type point. I tried to query it as an array like it says in the documentation here (http://www.postgresql.org/docs/current/static/functions-geometry.html):

It is possible to access the two component numbers of a point as
though the point were an array with indexes 0 and 1. For example, if
t.p is a point column then SELECT p[0] FROM t retrieves the X
coordinate and UPDATE t SET p[1] = … changes the Y coordinate. In
the same way, a value of type box or lseg can be treated as an array
of two point values.

Anyone know the correct query the planet_osm tables from osm2pgsql based on long/lat?
For example if I wanted to retrieve the information for Big Ben (London) I would query along the lines of:

    SELECT * FROM planet_osm_point 
    WHERE ST_Y(ST_Transform(way, 4326)) = '-0.12' 
      AND ST_X(ST_Transform(way, 4326)) = '51.5';

There is no error message, it comes out with the table names, but no data is output either. (Its an empty table). What am I doing wrong here?

Perhaps it is too precise, i.e. there are no coordinates that match it exactly? Perhaps I need to try being more general, say 51 < long < 52 and -0.5 < lat < 0. Although this would throw back multiple results so would not be quite as accurate.

Best Answer

To search for points around that location in London, you'd want a query like this

SELECT * FROM planet_osm_polygon WHERE ST_Intersects(way, ST_Buffer(ST_Transform(ST_SetSRID(ST_Point(-0.12,51.5),4326), 3857), 500); This will take the point, transform it into the web mercator projection (EPSG 3857) used by default in osm2pgsql, buffer it by 500 mercator meters, and find points that lie within that polygon.

A more accurate way to do this would be

SELECT * FROM planet_osm_polygon WHERE ST_Intersects(way, ST_Transform(ST_Buffer( ST_SetSRID(ST_Point(-0.12,51.5),4326)::geography,500)::geometry, 3857) );

This converts to geography so you can expand the point by 500 true meters.

A third, and probably quickest way, is

SELECT * FROM planet_osm_polygon WHERE ST_DWithin(way, ST_Transform(ST_SetSRID(ST_Point(-0.12,51.5),4326), 3857), 500);

Of course in practice you'll have some other WHERE conditions to apply, and be using a more complicated real-world geometry to search against.

Some common mistakes are

  • Trying to compare geometries of different projections
  • Reprojecting the geometry in the planet_osm_point table, preventing index usage
  • Mixing up latitude and longitude.

It's worth noting that PostGIS geometries are completely distinct from PostgreSQL geometric types. The appropriate page to see available functions is Chapter 8 of the PostGIS reference.

Related Question