PostGIS - Shapefiles

创建一个新数据库来存储地理空间数据：

psql pglab -c "create database geodb"

创建扩展：

psql geodb -c "create extension postgis;"
psql geodb -c "create extension pgrouting;"

确保扩展程序安装正确：

psql geodb -c "select * from pg_available_extensions where name in ('postgis','pgrouting');"

image-20231111023736603

下载Shapefile

shapefile 是指名称以 .shp、.shx、.dbf 结尾且具有公共前缀名称（例如 nyc_hotspots）的文件集合：

  • .shp—形状格式；特征几何本身
  • .shx——形状索引格式；特征几何的位置索引
  • .dbf——属性格式；

可选文件包括：

  • .prj——投影格式；坐标系和投影信息。

本节我们将：

  • 从公共 s3 下载纽约市出租车行程数据。我们将使用 2015 年 12 月创建的数据。
  • 从 NYC OpenData 平台下载 NYC 数据集。使用shp2pgsql程序将 shapefile 转换为几何格式插入到 PostgreSQL/PostGIS 数据库的 SQL

执行以下命令：

sudo yum -y install postgis curl
cd $HOME/environment
mkdir postgis && cd postgis

aws s3 cp s3://nyc-tlc/csv_backup/yellow_tripdata_2015-12.csv .

curl -L -o boros.zip https://data.cityofnewyork.us/api/geospatial/tqmj-j8zm?method=export\&format=Shapefile
unzip -o boros.zip -d boros

# `shp2pgsql`程序通过将形状数据从二进制数据转换为一系列 SQL 命令，然后在数据库中运行以加载数据，使形状数据可在 PostGIS 中使用。 `-s`标志告诉程序`spatial reference identifier (SRID)`数据是什么
shp2pgsql -c -D -s 4326-i -I boros/geo_export*.shp public.nyc_boros > nyc_boros.sql

curl -L -o subwaylines.zip https://data.cityofnewyork.us/api/geospatial/3qz8-muuu?method=export\&format=Shapefile
unzip -o subwaylines.zip -d subwaylines
shp2pgsql -c -D -s 4326 -i -I subwaylines/geo_export*.shp public.nyc_subwaylines > nyc_subwaylines.sql

curl -L -o nyc_subways.zip https://data.cityofnewyork.us/api/geospatial/arq3-7z49?method=export\&format=Shapefile
unzip -o nyc_subways.zip -d nyc_subways
shp2pgsql -c -D -s 4326 -i -I nyc_subways/geo_export*.shp public.nyc_subways > nyc_subways.sql


curl -L -o nyc_streets.zip  https://data.cityofnewyork.us/api/geospatial/exjm-f27b?method=export\&format=Shapefile
unzip -o nyc_streets.zip -d nyc_streets
shp2pgsql -c -D -s 4326 -i -I nyc_streets/geo_export*shp public.nyc_streetlines > nyc_streetlines.sql

curl -L -o neighbor.zip https://data.cityofnewyork.us/api/geospatial/cpf4-rkhq?method=export\&format=Shapefile
unzip -o neighbor.zip -d neighbor
shp2pgsql -c -D -s 4326 -i -I neighbor/geo_export*.shp public.nyc_neighborhood > nyc_neighborhood.sql

curl -L -o nyc_hotspots.zip https://data.cityofnewyork.us/api/geospatial/a9we-mtpn?method=export\&format=Shapefile
unzip -o nyc_hotspots.zip -d nyc_hotspots
shp2pgsql -c -D -s 4326 -i -I nyc_hotspots/geo_export*.shp public.nyc_hotspots > nyc_hotspots.sql

使用 SRID 查询 shapefile

我们可以使用Spatial Reference IDentifier (SRID)shapefile 来查询在线数据库以获取地图的定义。SRID 很方便，因为它将有关地图投影（可能非常复杂）的所有信息打包到一个数字中。

可以通过在在线数据库查找来查看地图投影的定义，打开https://epsg.io/4326 或者直接在 PostGIS 数据库中对表spatial_ref_sys进行查询：

psql geodb -c "select srtext from spatial_ref_sys where srid = 4326;"

image-20231111024309218

将空间数据文件导入 PostgreSQL

在数据库中创建表trips：

psql geodb
CREATE TABLE trips (
  id serial primary key,
  vendor_id text,
  pickup_datetime timestamp without time zone,
  dropoff_datetime timestamp without time zone,
  passenger_count numeric,
  trip_distance numeric,
  pickup_longitude numeric,
  pickup_latitude numeric,
  rate_code text,
  store_and_fwd_flag text,
  dropoff_longitude numeric,
  dropoff_latitude numeric,
  payment_type text,
  fare_amount numeric,
  surcharge numeric,
  mta_tax numeric,
  tip_amount numeric,
  tolls_amount numeric,
  improvement_surcharge numeric,
  total_amount numeric,
  pickup_location geometry(Point, 4326),
  dropoff_location geometry(Point, 4326)
);
\q

使用 shp2pgsql 命令生成的 sql 文件创建其他表:

psql geodb -f nyc_boros.sql
psql geodb -f nyc_subwaylines.sql
psql geodb -f nyc_subways.sql
psql geodb -f nyc_streetlines.sql
psql geodb -f nyc_neighborhood.sql
psql geodb -f nyc_hotspots.sql

image-20231111024412679

使用COPY命令将出租车行程数据加载到表中, 可能需要几分钟才能完成:

cd ~/environment/postgis
cat yellow_tripdata_2015-12.csv | psql geodb -c "COPY trips(vendor_id, pickup_datetime, dropoff_datetime, passenger_count, trip_distance, pickup_longitude, pickup_latitude, rate_code, store_and_fwd_flag, dropoff_longitude, dropoff_latitude, payment_type, fare_amount, surcharge, mta_tax, tip_amount, tolls_amount, improvement_surcharge, total_amount) FROM stdin CSV HEADER;"

image-20231111024801216

目前所有经纬度数据均为NUMBER格式，需要转换为PostGIS Geometry以供进一步分析。以下命令将更新整个表, 可能需要几分钟才能完成：

psql geodb
update trips set pickup_location = ST_SetSRID(ST_MakePoint(pickup_longitude, pickup_latitude), 4326),
dropoff_location = ST_SetSRID(ST_MakePoint(dropoff_longitude, dropoff_latitude), 4326);
\q

image-20231111102308332

在表trips的列上创建索引dropoff_datetime以加快搜索速度：

psql geodb
create index trips_dropoff_datetime on trips (date_trunc('day',dropoff_datetime));
\q

image-20231111102359195

该nyc_streetlines表加载了来自 geom 列的 MULTILINESTRING 数据。对于额外的路由计算，我们将添加额外的列, 以将 MULTILINESTRING 数据转换为 LINESTRING 数据

psql geodb
alter table nyc_streetlines add geomline geometry;
update nyc_streetlines set geomline = ST_LineMerge(geom);
\q

image-20231111102415342

在表中创建拓扑nyc_streetlines以绘制街道地图。pgRouting扩展了 PostGIS/PostgreSQL 地理空间数据库以提供地理空间路由和其他分析功能。pgRouting 的 pgr_createTopology功能提供了一种创建路由网络拓扑的方法：

  • 为每个道路链接分配源标识符和目标标识符。
  • 创建一个相关的顶点表。
  • 创建基本索引。
psql geodb
ALTER TABLE nyc_streetlines ADD COLUMN source integer;
ALTER TABLE nyc_streetlines ADD COLUMN target integer;
CREATE INDEX nyc_streetlines_source_idx ON nyc_streetlines (source);
CREATE INDEX nyc_streetlines_target_idx ON nyc_streetlines (target);
SELECT pgr_createTopology('nyc_streetlines', 0.0001, 'geomline', 'gid');
\q

image-20231111105058768

列出创建的所有表：

psql geodb -c "\dt"

image-20231111105116454

小结

在本节中我们下载了 shapefile，并使用 shp2pgsql 实用程序加载了 shapefile， 创建了多个表。