问题描述
我通过 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
的索引。
进一步阅读
- Getting all Buildings in range of 5 miles from specified coordinates
- Buffers (Circle) in PostGIS
- How can I set data from on table to another according spatial relation geometries in these tables
我想通了。
这就是你如何从 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 值……我不确定。