矢量地理数据导入PostgreSQL

发表于 分类于 阅读次数： Disqus：

前言

前一段时间，用Flask写了一个Web应用，需要用到一些地理空间的数据，为了更有效地管（zhuang）理（bi）起来，就想用前面提到过的PostgreSQL数据库，这里我们简单介绍一下，如何用Python的相关库，将常见的矢量数据（shp）导入到PostgreSQL中。

前期准备

项目描述

我们将从美国地质调查局（USGS）下载shapefiles，将其导入到PostgreSQL数据库中，能够使用GeoAlchemy 2 ORM 执行地理空间查询，使用SQLAlchemy执行多表关联查询，并在前段显示查询结果。

相关Python库

我们用到的相关库包括：

此外，在使用Flask搭建Web应用的时候还会用到其他的库，例如Jinja2模板、Werkzeug WSGI模块等，此处不在赘述。

下载数据

从USGS的网站上可以下载许多数据，我们这次用到四种数据：

  • 美国各州边界（US States）；
  • 美国各县边界（US County Boundaries）；
  • 美国国会选区（Congressional Districts）；
  • NBA球场位置（NBA Arenas）。

具体的下载地址见 https://www.sciencebase.gov/catalog/item/4f4e4a2ee4b07f02db615738

usgs-data-source

创建数据库

我们使用SQLAlchemy和GeoAlchemy2创建数据库和数据表来存储Web应用用到的数据。

首先导入需要用的模块和方法：

1
2
3
4
5
6
from sqlalchemy import create_engine
from sqlalchemy_utils import database_exists, create_database, drop_database
from sqlalchemy import Column, Integer, String, ForeignKey, Float, MetaData
from sqlalchemy.orm import relationship
from geoalchemy2 import Geometry
from sqlalchemy.ext.declarative import declarative_base

构造连接的字符串，格式如下：
conn_string = '{DBtype}://{user}:{password}@{instancehost}:{port}/{database}'

在创建PostgreSQL数据库之后，我们需要为新创建的数据库添加postgis扩展：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
conn_string = 'postgresql://postgres:postgres@localhost:5432/nba_arenas'
engine = create_engine(conn_string, echo=False)

# Check to ensure that the database doesn't exist
# If it doesn't, create it and generate the PostGIS extention and tables
if not database_exists(engine.url):
    create_database(engine.url)

    # Create a direct connection to the database using the engine.
    # This will allow the new database to use the PostGIS extension.

    conn = engine.connect()
    conn.execute('commit')
    try:
        conn.execute('CREATE EXTENSION postgis')
    except Exception as e:
        print(e)
        print('extension postgis already exists')
    conn.close()

定义数据库表

根据MVC的架构模式，数据库（表）属于模型：模型（Model） 用于封装与应用程序的业务逻辑相关的数据以及对数据的处理方法。“ Model ”有对数据直接访问的权力，例如对数据库的访问。“Model”不依赖“View”和“Controller”。

在Python中，我们通过继承已经预先封装了大部分数据库管理功能的超类编写新的类来定义数据表。通过封装我们不需要关心模型的实现代码就可以实现不同类型RDBMS的管理。SQLAlchemy模型可以用于产生包括SQL Server、Oracle、Postgres SQL和MySQL等不同的数据库表，而GeoAlchemy2则仅适用于PostgreSQL/PostGIS。

1
2
3
4
5
6
# For SQLAlchemy database classes, declarative_base
# allows for inheritance of database methods and properties

# Definde the model Base
Base = declarative_base()
metadata = MetaData()

通过继承基类，我们可以创建不同表的模型类，通过这些模型类分别建立相应的数据库表。其中__tablename__属性定义了数据库表的名称，通过关键字primary_key定义表的主键。

对于空间数据而言，我们还需要用srid定义几何形状（Geometry）的空间参考坐标系，例如srid=4326代表我们使用的是WSG 1984坐标系。

对于具有关系的表，例如County、District和Sate，我们用专门的类和方法来管理表的关系和多表查询。例如，ForeginKey 类为子表提供外键，relationship实现连表查询，backref是反向关联。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
metadata.clear()

# Define the Arena class, which will model the Arena database table
class Arena(Base):
    __tablename__ = 'arena'
    __table_args__ = {'extend_existing': True}
    id =  Column(Integer, primary_key=True)
    name = Column(String)
    longitude = Column(Float)
    latitude = Column(Float)
    geom = Column(Geometry(geometry_type='POINT'))

# Define the County Class
class County(Base):
    __tablename__ = 'county'
    __table_args__ = {'extend_existing': True}
    id = Column(Integer, primary_key=True)
    district = Column(String)
    name = Column(String)
    state_id = Column(Integer, ForeignKey('state.id'))
    geom = Column(Geometry(geometry_type='MULTIPOLYGON', srid=4326))

# Define the District class
class District(Base):
    __tablename__ = 'district'
    __table_args__ = {'extend_existing': True}
    id = Column(Integer, primary_key=True)
    district = Column(String)
    name = Column(String)
    state_id = Column(Integer, ForeignKey('state.id'))
    geom = Column(Geometry(geometry_type='MULTIPOLYGON', srid=4326))

# Define the State class
class State(Base):
    __tablename__ = 'state'
    __table_args__ = {'extend_existing': True}
    id = Column(Integer, primary_key=True)
    name = Column(String)
    statefips = Column(String)
    stpostal = Column(String)
    counties = relationship('County', backref='state')
    districts = relationship('District', backref='state')
    geom = Column(Geometry(geometry_type='MULTIPOLYGON', srid=4326))


传统的方法，是在父类中定义一个关系relationship或叫正向引用 Reference，子类只需定义一个外键。单纯定义的relationship(‘子类名’)只是一个正向引用，也就是只能让父类调用子对象。反过来，如果要问children他们的父亲是谁，就不行了。所以，我们还需要一个反向引用 (Back Reference)的声明，让子对象能够知道父对象是谁。

在上面，我们用relationship建立了两个表的关联，并且用backref实现了反向关联，以

state_id = Column(Integer, ForeignKey('state.id')) counties = relationship('County', backref='state')

为例, state为主表，county为从表：

  • 查询state表，返回state_obj，可以通过state_obj.counties查询到county表中外键关联的数据，查主表返回从表里的数据，这个称之为“正向查询”。
  • 查询county表，返回county_obj，可以通过county_obj.state查询到外键关联的state表中数据。查从表返回主表里的数据，这个称之为：“反向查询”。

创建数据库表

我们可以用两种方法来创建数据库表：

  • 一种是调用模型类的 __table__方法的create()函数，分别创建每一个表；
  • 另一种是调用基类Base的metadata方法的create_all()方法，一次性创建所有继承自Base的表。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# Generate the District table from the District class.
# If it already exists, drop it and regenerate it
# try:
#     District.__table__.create(engine)
# except:
#     District.__table__.drop(engine)
#     District.__table__.create(engine)

# # Generate the County table from the County class.
# # If it already exists, drop it and regenerate it
# try:
#     County.__table__.create(engine)
# except:
#     County.__table__.drop(engine)
#     County.__table__.create(engine)


# # Generate the Arena table from the Arena class.
# # If it already exists, drop it and regenerate it
# try:
#     Arena.__table__.create(engine)
# except:
#     Arena.__table__.drop(engine)
#     Arena.__table__.create(engine)

# # Generate the State table from the State class.
# # If it already exists, drop it and regenerate it
# try:
#     State.__table__.create(engine)
# except:
#     State.__table__.drop(engine)
#     State.__table__.create(engine)

# To create all the tables at once, comment out the try/except blocks above
# and uncomment the line below:
Base.metadata.create_all(engine)

pgAdmin

导入数据

利用pgAdmin查看PostgreSQL数据库可以看到，我们已经成功创建了数据库及所需要用的数据表，接下来我们就将下载好的数据写入到数据库中。

导入模块

为了从shapefile中读取数据，我们需要用到pyshp和pygeoif模块，其中pyshp用来读取shapefile中的空间（属性）数据，pygeoif模块则是将地理空间数据转换为Python对象，我们将用它将shapefile中的二进制数据转换为可以被GeoAlchemy2处理的WKT格式。

1
2
3
4
# The pyshapefile module is used to read shapefiles and
# the pygeoif module is used to convert between geometry types
import shapefile
import pygeoif

读取数据以后，我们还要用SQLAlchemy和GeoAlchemy2模块将数据导入到数据库中，除了前面导入的模块以外，我们还需要导入一些新的模块：

1
from sqlalchemy.orm import sessionmaker

使用tkinter提供GUI：

1
2
3
# The built-in Tkinter GUI module allows for file dialogs
from tkinter import filedialog
from tkinter import Tk

我们创建一个session，与engine绑定，用来管理查询和写入。

1
2
Session = sessionmaker(bind=engine)
session = Session()

加载shp数据

打开shapefiles

1
2
3
4
5
6
7
# Read the Arena shapefile using the Reader class of the pyshp module
arena_shp = open('./data/Arenas_NBA/Arenas_NBA.shp', 'rb')
arena_dbf = open('./data/Arenas_NBA/Arenas_NBA.dbf', 'rb')
arena_reader = shapefile.Reader(shp=arena_shp, dbf=arena_dbf)
arena_shapes = arena_reader.shapes()
arena_records = arena_reader.records()
print(arena_reader)
shapefile Reader
    29 shapes (type 'POINT')
    29 records (30 fields)

对county、district和state执行类似的操作，不再重复。

写入数据库

在这里我们通过循环将shapefile中的每条数据顺序读取出来，对于几何属性数据，采用WKT的格式保存：

  • 对于简单的数据（POINT），我们手动构造WKT；
  • 对于复杂的数据（POLYGON)，使用pygeoif先转换为pygeofi MultiPolygon格式，在转换为WKT。

下面是部分代码作为示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
# Iterate through the Arena data read from the shapefile
for count, record in enumerate(arena_records):

    # Instantiate an arena from the Arena class and populate the fields
    arena = Arena()
    arena.name = record[6]
#     print(arena.name)

    # Get the district geometry associated with the current record and use
    # indexing to get the longitude and latitude for the point.
    # Insert the coordinates into a Well Known Text string template using string formatting
    point = arena_shapes[count].points[0]
    arena.longitude = float(point[0])
    arena.latitude = float(point[1])
    arena.geom = 'SRID=4326;POINT({0} {1})'.format(arena.longitude, arena.latitude)
    session.add(arena)
session.commit()


# Iterate through the District data read from the shapefile
# This uses the STFIPS data to query the State table and find the related state
for count, record in enumerate(district_records):

    # Instantiate a district from the District class and populate the fields
    district = District()
    district.district = record[0]
    district.name = record[1]
    state = session.query(State).filter_by(statefips=record[4]).first()
    district.state_id = state.id
#     print(district.name, district.district)

    # Get the district geometry associated with the current record and convert it
    # into a pygeoif MultiPolygon, then into Well Known Text
    # for insertion into the data table
    district_geo = district_shapes[count]
    gshape = pygeoif.MultiPolygon(pygeoif.geometry.as_shape(district_geo))
    district.geom = 'SRID=4326;{0}'.format(gshape.wkt)

    session.add(district)
    if count % 50 == 0:
        session.commit()
session.commit()

# Close the session and dispose of the engine connection to the database
session.close()
engine.dispose()

结果测试

1
2
3
4
engine = create_engine(conn_string)
Session = sessionmaker(bind=engine)
session = Session()
Base = declarative_base()

简单查询

我们可以从数据库中查询所有的体育馆信息，然后将名字和经纬度存储到pandas的DataFrame中。

1
import pandas as pd
1
2
3
4
arenas = session.query(Arena).all()
arenas_list = pd.DataFrame([(arena.id, arena.name, arena.longitude, arena.latitude)
                            for arena in arenas], columns=['id', 'name', 'lon', 'lat'])
arenas_list.head()
id name lon lat
0 1 American Airlines Arena -80.187580 25.781170
1 2 ARCO Arena -121.518690 38.649320
2 3 EnergySolutions Arena -111.901025 40.768206
3 4 Fed Ex Forum -90.051300 35.138100
4 5 Quicken Loans Arena -81.688700 41.496480

空间查询

既然是空间数据，而且我们也为PostgreSQL添加了PostGIS扩展，所以我们测试一下空间查询的结果，在这里我们分别测试空间包含（ST_Contains）和空间相交（ST_Intersects），为了实现这一查询，我们用到了SQLAlchemy的filter方法和GeoAlchemy2 ORM定义的Geometry列。

1
2
3
4
5
6
7
8
arena = arenas[0]
county=session.query(County).filter(
    County.geom.ST_Contains(arena.geom)).first()
print(county.name)
if county != None:
    district=session.query(District).filter(
        District.geom.ST_Intersects(county.geom)).first()
print(district.name)
Miami-Dade
Alcee L. Hastings

关联查询

前面，我们以state为主表、county/district为从表，建立了关联，下面我们可以看一下如何进行关联查询。我们使用“state”作为backref的参数，所以可以通过{county_obj}.state实现逆向查询。而正向查询则更为简单，直接调用{state_obj}.counties即可。

1
2
3
state = county.state
print(state.name)
print([sc.name for sc in state.counties])
Florida
['Hardee', 'Okeechobee', 'Hillsborough', 'Indian River', 'St. Lucie', 'Manatee', 'DeSoto', 'Glades', 'Martin', 'Miami-Dade', 'Lake', 'Sumter', 'Seminole', 'Volusia', 'Orange', 'Citrus', 'Hernando', 'Polk', 'Osceola', 'Pasco', 'Brevard', 'Pinellas', 'Highlands', 'Hendry', 'Sarasota', 'Charlotte', 'Lee', 'Palm Beach', 'Broward', 'Collier', 'Monroe', 'Jackson', 'Holmes', 'Washington', 'Gadsden', 'Madison', 'Hamilton', 'Calhoun', 'Liberty', 'Columbia', 'Baker', 'Nassau', 'Suwannee', 'Okaloosa', 'Walton', 'Santa Rosa', 'Escambia', 'Leon', 'Duval', 'Lafayette', 'Bradford', 'Union', 'Jefferson', 'Clay', 'Wakulla', 'Bay', 'Alachua', 'Gilchrist', 'Gulf', 'St. Johns', 'Taylor', 'Franklin', 'Marion', 'Flagler', 'Putnam', 'Dixie', 'Levy']

小结

我们简单介绍了从shapefile中读取数据，然后存储到PostgreSQL中，并且可以实现空间查询和多表关联查询，可以看到只要了解空间数据的结构，同样可以将其他格式的数据例如GeoJSON，存储到PosttgreSQL中。

luoxiaohei