shapefile 数据导入及空间查询

导入数据

范例数据包含纽约市的街道、地铁站等信息，路径: nyc_shapefile

使用 shp2pgsql导入数据

导入 nyc_census_blocks 数据

shp2pgsql \
  -D \
  -I \
  -s 26918 \
  nyc_census_blocks.shp \
  nyc_census_blocks \
  | gsql dbname=yukontutorial

其中：

-D: 使用 dump format 方式，这个要比默认的 insert format 快很多。
-I: 加载数据完成后，为 geometry 类型的列创建一个索引
-s: 数据的 srid 值。
nyc_census_blocks.shp: 要加载的数据集
nyc_census_blocks: 加载后在数据库中创建的表名
gsql dbname=yukontutorial : 将数据导入到 yukontutorial 数据库中

当你看到有如下输出时，则说明数据导入成功

Field popn_total is an FTDouble with width 11 and precision 0
Field popn_white is an FTDouble with width 11 and precision 0
Field popn_black is an FTDouble with width 11 and precision 0
Field popn_nativ is an FTDouble with width 11 and precision 0
Field popn_asian is an FTDouble with width 11 and precision 0
Field popn_other is an FTDouble with width 11 and precision 0
Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]
SET
SET
BEGIN
NOTICE:  CREATE TABLE will create implicit sequence "nyc_census_blocks_gid_seq" for serial column "nyc_census_blocks.gid"
CREATE TABLE
NOTICE:  ALTER TABLE / ADD PRIMARY KEY will create implicit index "nyc_census_blocks_pkey" for table "nyc_census_blocks"
ALTER TABLE
                         addgeometrycolumn
--------------------------------------------------------------------
 public.nyc_census_blocks.geom SRID:26918 TYPE:MULTIPOLYGON DIMS:2
(1 row)

CREATE INDEX
COMMIT
ANALYZE

然后重新连接到 yukontutorial 数据库，输入 \d 查看新创建的表：

yukontutorial=# \d
                                    List of relations
 Schema |           Name            |   Type   | Owner |             Storage
--------+---------------------------+----------+-------+----------------------------------
 public | geography_columns         | view     | omm   |
 public | geometry_columns          | view     | omm   |
 public | geomodel_columns          | view     | omm   |
 public | nyc_census_blocks         | table    | omm   | {orientation=row,compression=no}
 public | nyc_census_blocks_gid_seq | sequence | omm   |
 public | raster_columns            | view     | omm   |
 public | raster_overviews          | view     | omm   |
 public | spatial_ref_sys           | table    | omm   | {orientation=row,compression=no}
(8 rows)

同理，依次导入剩余的数据集：

nyc_neighborhoods.shp
nyc_streets.shp
nyc_subway_stations.shp

你可能会好奇 .shp 文件到底是什么？

通常来说 shp 文件一般指的是 .shp, .shx,.dbf 还有一些其他类型文件且前缀名称一致的文件集合。shp 即 shape file ,通常单独的一个 shp 文件时没有什么用的，必须和其他类型的文件结合使用。

必要文件：

.shp：存储地理要素的几何信息
.shx：存储要素几何图形的索引信息
.dbf: 存储地理要素的属性信息（非几何信息）

可选文件：

.prj：存储空间参考信息，即地理坐标系统信息和投影坐标系统信息。使用 well-known 文本格式进行描述。

那 SRID 又是什么？

大多数导入过程都是不言自明的，但即使是经验丰富的 GIS 专业人员也可能被 SRID 难倒。

SRID 表示 Spatial Reference IDentifier（空间参考标识符)。它定义了我们数据的地理坐标系统和投影的所有参数。

SRID 很方便，因为它将有关地图投影的所有信息（可能非常复杂）打包（更具体的说应该是映射）到一个数字中。

你可以在以下链接中查找我们在上面使用的投影的定义：SRID

数据操作

我们想知道纽约市街道的总长度是多少应该怎么做？

select SUM(ST_LENGTH(GEOM)) from nyc_streets;

你可以看到如下输出：

yukontutorial=# select  SUM(ST_LENGTH(GEOM)) from nyc_streets ;
      sum
---------------
 10418904.7172
(1 row)

这样我们就可以知道纽约市街道总长度为 10418904.7172m

ST_Length(geometry geom)：如果 geom 的类型是 LineString 或者 MultiLineString,那么就返回这个 geometry 的长度。

sum：聚合函数，可以计算 nyc_streets 表中每一行 geometry 对象长度的总和。

我们想知道纽约市最西边的地铁站是哪一个？

select st_x(geom),name from nyc_subway_stations order by st_x(geom) limit 1;

输出结果：

yukontutorial=# select st_x(geom),name from nyc_subway_stations order by st_x(geom) limit 1;
       st_x       |    name
------------------+-------------
 563292.117258056 | Tottenville
(1 row)

ST_X(geometry a_point);: 如果 geometry 是一个点，则返回它的 x 坐标，否则返回 NULL。

现在我们知道位于纽约市最西边的地铁站是 Tottenville。

我们想知道距离 Broad St 地铁站 10m 范围内的街道有那些？

我们可以先查询一下 Broad St 地铁站的具体位置：

select st_astext(geom) from nyc_subway_stations where name='Broad St';

输出结果：

yukontutorial=# select st_astext(geom) from nyc_subway_stations where name='Broad St';
                st_astext
------------------------------------------
 POINT(583571.905921312 4506714.34119218)
(1 row)

现在我们知道了 Borad St 地铁站的具体地点为 POINT(583571.905921312 4506714.34119218)，接下来可以查询距离地铁站 10m 范围内的街道了。

SELECT name
FROM nyc_streets
WHERE ST_DWithin(
        geom,
        ST_GeomFromText('POINT(583571.905921312 4506714.34119218)',26918),
        10
      );

输出结果：

   name
-----------
 Wall St
 Broad St
 Nassau St
(3 rows)

现在我们可以知道距离 Broad St 地铁站 10m 内的街道有这 3 个。

ST_AsText(geometry g1): 返回 WKT 形式的 geometry 数据。

ST_GeomFromText(text WKT, integer srid); : 返回由 WKT 形式表示的 geometry 数据类型。同时指定它的 SRID。

ST_DWithin(geometry g1, geometry g2, double precision distance_of_srid) : 如果 g1 和 g2 相距 distance_of_srid 之内，则返回真，否则返回假。

我们想知道纽约市人口密度最大的社区是哪一个？

SELECT
  n.name,
  Sum(c.popn_total) / (ST_Area(n.geom) / 1000000.0) AS popn_per_sqkm
FROM nyc_census_blocks AS c
JOIN nyc_neighborhoods AS n
ON ST_Intersects(c.geom, n.geom)
GROUP BY n.name, n.geom
ORDER BY 2 DESC limit 1;

输出结果：

       name        |  popn_per_sqkm
-------------------+------------------
 North Sutton Area | 68435.1328377268
(1 row)

从结果中我们可以看到 North Sutton Area 的人口密度是最大的。

ST_Intersects( geometry geomA , geometry geomB ): 返回 geomA 和 geomB 是否相交