1

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.

  • Welcome to GIS.SE. Could you show us -- via the edit button, above -- what query you used to get that error? (Also, your long-lat for London are reversed.) – Martin F Apr 23 '14 at 17:28
  • I used this query: select * from planet_osm_point where way[0] = '51.5' and way[1] = '-0.12' – Quanqai Apr 23 '14 at 17:32
  • Still looks like you're mixing up coordinate associations: X <=> long. Y <=> lat. – Martin F Apr 23 '14 at 19:48
  • Ah, yeah, I tried it both ways to no avail. 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. – Quanqai Apr 23 '14 at 19:53
1

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.

  • Thank you for this explanation, really cleared up some things that were confusing me. Plus I like the method of expanding the polygon from a central point, much easier than me having to write methods to get upper and lower bounds from a central point. – Quanqai Apr 27 '14 at 12:01
  • As a general rule, if you ever use ST_X or ST_Y in a WHERE clause, you're doing something wrong and will have lost a lot of performance. – Paul Norman Apr 29 '14 at 5:47
  • 1
    Was wondering why I couldn't get an output as I just substituted these statements into my Java code. So I tried putting them straight into pgAdmin III and I got the error, in all cases, ERROR: Operation on mixed SRID geometries SQL state: XX000. Do you know what causes this? – Quanqai Apr 30 '14 at 18:33
  • Okay, this seems to be a unique issue. I have played around with the statements some more. I can change the longitude to whatever and it works, but for some reason if I change the latitude so that it is larger than 52.4 it throws the mixed SRID error. – Quanqai Apr 30 '14 at 18:51
0

Note that while we usually talk of latitude and longitude, in that order, there are different conventions for this ordering, even within GIS, cartography and surveying.

The convention in PostGIS, and others, is longitude then latitude relating to X and Y, respectively.

Thus, for London try:

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

But that assumes you will have an exact match of coordinate values. If you wish to do a range query you need to specify a search box or circle. For example:

What's the fastest way to do a bounding box query in PostGIS

Problem with using a circle as polygon for querying PostGIS

How to find nearby objects using Easting and Northing range

  • These links are useful but I am confused by the ST_MakeEnvelope method - I tried the following statement: SELECT *from planet_osm_point where way && ST_MakeEnvelope(51.4, 0.1, 51.6, 0.2); and still there were no points output (also tried with @ as well as &&). – Quanqai Apr 23 '14 at 20:37
  • Again, check your X-Y order, and sign. – Martin F Apr 23 '14 at 20:42
  • Ah, I got confused between what yours said and what the answerer beneath yours said (with ST_MakeEnvelope(l,b,r,y,SRID)). Yours was clearer. I think this is right: select *from planet_osm_point where way && ST_MakeEnvelope(-0.2, 51, -0.1, 52); However still no output but if I do it manually so that it is SELECT * FROM planet_osm_point WHERE (ST_Y(ST_Transform(way,4326)) > 51 AND ST_Y(ST_Transform(way,4326)) < 52) AND (ST_X(ST_Transform(way,4326)) > -0.2 AND ST_X(ST_Transform(way,4326)) < -0.1) I get the right output. – Quanqai Apr 23 '14 at 20:54
  • Doing ST_X(...)=... and ST_Y(...)=... checks will prevent any index usage, and kill any performance gains you would have gotten from indexes. – Paul Norman Apr 26 '14 at 8:08

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.