rpostgresqlpostgisopenstreetmaposm2pgsql

Query all primary roads from PostGis Database filled with OSM Data within a boundary

I imported OpenStreetMap data through osm2pgsql into PgSQL (PostGIS) and I would like to get an SF object from the data containing all primary roads (geometry) within a an area (bbox) into R.

I got lost since I would like to have also relations and nodes and im not sure if only a query on planet_osm_roads will be sufficient and how the imported data structure is different to osm xml data im normaly working with. I understand it is probably a bit broader question but

I would appreciate a start so to say to understand the query language better.

This is my approach but sadly i get an error

conn <- RPostgreSQL::dbConnect("PostgreSQL", host = MYHOST,
                               dbname = "osm_data", user = "postgres", password = MYPASSWORD)
pgPostGIS(conn)


a<-pgGetGeom(conn, c("public", "planet_osm_roads"), geom = "way", gid = "osm_id",
          other.cols = FALSE, clauses  = "WHERE highway = 'primary' && ST_MakeEnvelope(11.2364353533134, 47.8050651144447,  11.8882527806375, 48.2423300001326)")
a<-st_as_sf(a)

This is an error i get:

Error in postgresqlExecStatement(conn, statement, ...) : 
  RS-DBI driver: (could not Retrieve the result : ERROR:  syntax error at or near "ST_MakeEnvelope"
LINE 2: ...lic"."planet_osm_roads" WHERE "way" IS NOT NULL   ST_MakeEnv...
                                                             ^
)
Error in pgGetGeom(conn, c("public", "planet_osm_roads"), geom = "way",  : 
  No geometries found.
In addition: Warning message:
In postgresqlQuickSQL(conn, statement, ...) :
  Could not create execute: SELECT DISTINCT a.geo AS type 
                        FROM (SELECT ST_GeometryType("way") as geo FROM "public"."planet_osm_roads" WHERE "way" IS NOT NULL   ST_MakeEnvelope(11.2364353533134, 47.8050651144447,  11.8882527806375, 48.2423300001326)) a;

This is the db:

osm_data=# \d
                  List of relations
  Schema  |        Name        |   Type   |  Owner   
----------+--------------------+----------+----------
 public   | geography_columns  | view     | postgres
 public   | geometry_columns   | view     | postgres
 public   | planet_osm_line    | table    | postgres
 public   | planet_osm_nodes   | table    | postgres
 public   | planet_osm_point   | table    | postgres
 public   | planet_osm_polygon | table    | postgres
 public   | planet_osm_rels    | table    | postgres
 public   | planet_osm_roads   | table    | postgres
 public   | planet_osm_ways    | table    | postgres
 public   | spatial_ref_sys    | table    | postgres
 topology | layer              | table    | postgres
 topology | topology           | table    | postgres
 topology | topology_id_seq    | sequence | postgres




schema_name         table_name geom_column geometry_type     type
1      public    planet_osm_line         way    LINESTRING GEOMETRY
2      public   planet_osm_point         way         POINT GEOMETRY
3      public planet_osm_polygon         way      GEOMETRY GEOMETRY
4      public   planet_osm_roads         way    LINESTRING GEOMETRY




                      Table "public.planet_osm_roads"
       Column       |           Type            | Collation | Nullable | Default 
--------------------+---------------------------+-----------+----------+---------
 osm_id             | bigint                    |           |          | 
 access             | text                      |           |          | 
 addr:housename     | text                      |           |          | 
 addr:housenumber   | text                      |           |          | 
 addr:interpolation | text                      |           |          | 
 admin_level        | text                      |           |          | 
 aerialway          | text                      |           |          | 
 aeroway            | text                      |           |          | 
 amenity            | text                      |           |          | 
 area               | text                      |           |          | 
 barrier            | text                      |           |          | 
 bicycle            | text                      |           |          | 
 brand              | text                      |           |          | 
 bridge             | text                      |           |          | 
 boundary           | text                      |           |          | 
 building           | text                      |           |          | 
 construction       | text                      |           |          |
Solution

  • Your query looks just fine. Check the following example:

    WITH planet_osm_roads (highway,way) AS (
  VALUES 
    ('primary','SRID=3857;POINT(1283861.57 6113504.88)'::geometry), --inside your bbox
    ('secondary','SRID=3857;POINT(1286919.06 6067184.04)'::geometry)  --somewhere else ..
)
SELECT highway,ST_AsText(way)
FROM planet_osm_roads
WHERE 
   highway IN ('primary','secondary','tertiary') AND
   planet_osm_roads.way && ST_Transform(
  ST_MakeEnvelope(11.2364353533134,47.8050651144447,11.8882527806375,48.2423300001326, 4326),3857
);

 highway |          st_astext           
---------+------------------------------
 primary | POINT(1283861.57 6113504.88)

    This image illustrates the BBOX and the points used in the example above

    enter image description here

    However, there are a few things to consider.

    A partial index is an index built over a subset of a table; the subset is defined by a conditional expression (called the predicate of the partial index). The index contains entries only for those table rows that satisfy the predicate. Partial indexes are a specialized feature, but there are several situations in which they are useful.

    So try something like this:

    CREATE INDEX idx_planet_osm_roads_way ON planet_osm_roads USING gist(way) 
WHERE highway IN ('primary','secondary','tertiary');

    And also highway needs to be indexed. So try this ..

    CREATE INDEX idx_planet_osm_roads_highway ON planet_osm_roads (highway);

    .. or even another partial index, in case you can't delete the other data but you don't need it for anything:

    CREATE INDEX idx_planet_osm_roads_highway ON planet_osm_roads (highway) 
WHERE highway IN ('primary','secondary','tertiary');

    You can always identify bottlenecks and check if the query planer is using your index with EXPLAIN.

    Further reading