2

I'm trying to write a query that will return all OSM ways which are contained within a bounding box of min/max gps coordinates.

I get these cooridinates from the greater london bounding box on http://wiki.openstreetmap.org/wiki/Bounding_Box

From googling I have this query, but it doesn't return any results:

SELECT 
    planet_osm_roads.* 
FROM planet_osm_roads 
WHERE ST_Contains( 
    ST_SetSRID(
        ST_MakeBox2D(
            ST_Point(0.3340155, 51.2867602), ST_Point(-0.5103751, 51.6918741)
        ), 
        4326
    ), 
    planet_osm_roads.way
)

Any idea what I'm doing wrong?

Update

So I've swapped the points around, tested that the polygon created is what I expect by getting the geometry as geojson and visualising it on geojson.io - it looks fine, however I still get 0 results and it happens basically instantly (which seems weird). The table is definitely fully populated! I also changed to use ST_Intersects but that hasn't made it return anything either

SELECT 
    planet_osm_roads.* 
FROM planet_osm_roads 
WHERE ST_Intersects(
    planet_osm_roads.way, 
    ST_SetSRID(
        ST_MakeBox2D(
    ST_Point(0.3340155, 51.2867602),
            ST_Point(-0.5103751, 51.6918741)
        ), 
        4326
    )
)

It gives me this GeoJSON when I SELECT ST_AsGeoJSON

{
  "type": "Polygon",
  "coordinates": [
    [
      [-0.5103751, 51.2867602],
      [-0.5103751, 51.6918741],
      [0.3340155, 51.6918741],
      [0.3340155, 51.2867602],
      [-0.5103751, 51.2867602]
    ]
  ]
}

Which as I said, looks correct on geojson.io :)

Update #2

I've copied a query from another SO question that gets the closest point on a road or something, with a point in London, it returns some results roughly as I expect

SELECT 
 osm_id,
 highway,
 name,
 ref, 
 ST_X(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom)),      
 ST_Y(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom)),
 ST_Distance_Sphere(ST_ClosestPoint(ST_Transform(r.way, 4326), point.geom), point.geom)
FROM 
  planet_osm_roads r,
  (Select ST_SetSRID(ST_MakePoint(-0.1367229, 51.5377733), 4326) as geom) point
ORDER BY 7 ASC
LIMIT 10;

returns

41383228;"trunk";"Camden Street";"A400";-0.137183850111055;51.5373791549638;54.195918203
253353304;"primary";"Royal College Street";"A5202";-0.13568782069246;51.538616921758;118.002625807
17513372;"primary";"Royal College Street";"A5202";-0.134746843460286;51.5374616146768;140.996741591
166369355;"primary";"Royal College Street";"A5202";-0.136110935219222;51.5391336013284;157.068830304
253352352;"primary";"Royal College Street";"A5202";-0.136578148998492;51.5397000983131;214.484260552
6247747;"primary";"Pratt Street";"A503";-0.140021211819466;51.5376420920654;228.588127843
211950233;"primary";"Bayham Street";"A503";-0.140021211819466;51.5376420920654;228.588127843
253350830;"primary";"St. Pancras Way";"A5202";-0.135177855133607;51.5396274637255;232.22039167
11266300;"primary";"Royal College Street";"A5202";-0.136734905015571;51.5398841985372;234.722999647
253350834;"primary";"St Pancras Way";"A5202";-0.135399101773244;51.5397738652997;240.557396864

The point is actually closest to Mandela Street which runs next to Camden Street, but I'll worry about that after I have my initial query working.

My data came from download england osm data and running this command:

osm2pgsql --create --database england_osm -C 6144 ~/Downloads/england-latest.osm

googlemap

7
  • 1
    Swap the two points in ST_MakeBox2D? Doc says the first point should be the lower left one.
    – mkennedy
    Commented Jan 30, 2015 at 22:28
  • 1
    @mkennedy. The docs do say that, but it makes no difference, in practice. Check out, select st_astext(st_makebox2d(st_makepoint(50,50), st_makepoint(5,5)) ); Commented Jan 30, 2015 at 23:36
  • What Christopher says is true, but I am surprised that no roads are contained in a polygon covering all of London -- I hope so, I'm going there in a few hours by bus :D. The query looks fine, so perhaps you could post a couple of rows of the actual geometry values? Commented Jan 30, 2015 at 23:41
  • Thanks guys, I've updated with observations from your comments :) Commented Jan 31, 2015 at 10:58
  • Would you be able to provide the link the you got england-latest.osm from, so we are dealing with the same data. I fancy loading osm for the UK and might have some insights after doing so. Commented Feb 2, 2015 at 17:02

3 Answers 3

1

The query looks good but if you want all ways which are partly in a bounding box you have to use the st_intersects function.

If you want to show the error of your defined bounding box and you don't have any software you can put you bounding box in an ST_AsGeoJson function to get an geojson string as result. The easiest way to check your geometry is to enter the geojson result in an online viewer like http://geojson.io/

1
  • Thanks, useful tool! I've changed to ST_Intersects but still no luck! Commented Jan 31, 2015 at 10:58
1

So, it was because I hadn't specified that points should be stored in latitude and longitude when running osm2pgsql.

The command I used successfully was

osm2pgsql -c -l -s -d england_osm ~/Downloads/england-latest.osm.pbf
1
  • Ha ha. Good lesson. I have made similar mistakes so many times Commented Feb 6, 2015 at 14:01
0
-- roads inside bounding box
SELECT l.osm_id
FROM planet_osm_line l,
  (
    SELECT ST_Transform(ST_MakeEnvelope(min_x, min_y, max_x, max_y, 4326),3857) AS box
    FROM (
           SELECT
             least(x1, x2)    AS min_x,
             greatest(x1, x2) AS max_x,
             least(y1, y2)    AS min_y,
             greatest(y1, y2) AS max_y
           FROM (SELECT
                   :lon1 AS x1, 
                   :lon2 AS x2,
                   :lat1  AS y1,
                   :lat2  AS y2) cr
         ) r) bb
WHERE l.way && bb.box -- intersection (use @ for contain)
  and l.highway IS NOT NULL;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.