从边界内填充 OSM 数据的 PostGis 数据库查询所有主要道路

如何解决从边界内填充 OSM 数据的 PostGis 数据库查询所有主要道路

我通过 osm2pgsql 将 OpenStreetMap 数据导入到 Pgsql (PostGIS) 我想从包含的数据中获取一个 SF 对象一个区域（bBox）内的所有主要道路（几何）都变成了R。

我迷路了，因为我也想拥有关系和节点并且我不确定仅对planet_osm_roads 进行查询是否就足够了，以及导入的数据结构与我通常使用的osm xml 数据有何不同。我知道这可能是一个更广泛的问题，但

如果能更好地理解查询语言，我将不胜感激。

这是我的方法，但遗憾的是我遇到了错误

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)

这是我得到的错误：

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,: 
  No geometries found.
In addition: Warning message:
In postgresqlQuicksql(conn,...) :
  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,48.2423300001326)) a;

这是数据库：

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                      |           |          |

解决方法

您的查询看起来不错。检查以下示例：

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)

此图说明了 BBOX 和上面示例中使用的点

查看 documentation 以了解有关 bbox 交集运算符 && 的更多信息。

但是，有一些事情需要考虑。

如果您自己创建 BBOX 以便为 ST_Contains 提供一个区域，请考虑简单地使用 ST_DWithin。它基本上会做同样的事情，但你只需要提供一个参考点和距离。
根据表 highway 中 planet_osm_roads 类型的分布，并考虑到您的查询将总是寻找 primary、secondary或 tertiary 高速公路，创建 partial index 可以显着提高性能。正如文档所说：

部分索引是建立在表子集上的索引；子集由条件表达式定义（称为部分索引）。索引只包含那些表行的条目满足谓词。部分索引是一个特殊的功能，但在几种情况下它们很有用。

所以尝试这样的事情：

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

而且 highway 也需要被索引。所以试试这个..

CREATE INDEX idx_planet_osm_roads_highway ON planet_osm_roads (highway);

.. 甚至另一个部分索引，以防您无法删除其他数据但您不需要它：

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

您始终可以识别瓶颈并检查查询规划器是否使用带有 EXPLAIN 的索引。

进一步阅读

我想通了。这就是你如何从 PostGIS 数据库中获取一个 SF 对象，该数据库在 11.2364353533134,48.2423300001326 BBOX 中填充了 OSM 数据：

library(sf)
library(RPostgreSQL)
library(tictoc)
pw <- MYPASSWORD
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv,dbname = "osm_data",host = "localhost",port = 5432,user = "postgres",password = pw)


tic()
sf_data = st_read(con,query = "SELECT osm_id,name,highway,way
  FROM planet_osm_roads
WHERE highway = 'primary' OR highway = 'secondary' OR highway = 'tertiary'AND ST_Contains(
  ST_Transform(
    ST_MakeEnvelope(11.2364353533134,3857),planet_osm_roads.way);")
toc()

RPostgreSQL::dbDisconnect(con)

我必须验证是否真的考虑了 BBOX 值……我不确定。