从 PostGis 数据库中查询边界内充满 OSM 数据的所有主要道路
Query all primary roads from PostGis Database filled with OSM Data within a boundary
我通过 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, 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;
这是数据库：
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 | | |
我明白了。
这就是如何从 11.2364353533134,47.8050651144447,11.8882527806375,48.2423300001326 BBOX:
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,47.8050651144447,11.8882527806375,48.2423300001326,4326)
,3857)
,planet_osm_roads.way);")
toc()
RPostgreSQL::dbDisconnect(con)
我必须验证是否真的考虑了 BBOX 值。我不确定。
您的查询看起来很好。检查以下示例：
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, consider simply usingST_DWithin提供一个区域。它基本上会做同样的事情，但你只需要提供一个参考点和距离。 - 取决于
highway类型在 tableplanet_osm_roads中的分布，并且考虑到您的查询 总是 寻找primary、secondary或tertiary高速公路，创建partial index可以显着提高性能。正如文档所说：
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.
所以尝试这样的事情：
CREATE INDEX idx_planet_osm_roads_way ON planet_osm_roads USING gist(way)
WHERE highway IN ('primary','secondary','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','secondary','tertiary');
您始终可以识别瓶颈并检查查询规划器是否正在使用您的索引 EXPLAIN。
进一步阅读