GeoDjango是一款支持地理数据存储和空间查询的Django衍生项目。 这份教程假定你熟知Django,如果你是刚刚接触Django,那么请您先阅读 Django 教程 以获得Django的一些基础知识。同时,GeoDjango运行除了Django之外,仍 需要一些其他条件。– 请参考 GeoDjango 安装文档 以获详情。
这份教程将利用 世界地图 这个GIS应用来指导你的学习。 [1] 这份教程中部分代码受到了或者直接来自 GeoDjango 基础实例应用 这个项目的启发。 [2]
首先,项目需要创建一个空间数据库。如果利用PostgreSQL和PostGIS,则可利用如下的 命令从 空间数据库模板 (spatial database template) 来建立数据库:
$ sudo su - postgres
$ createdb -T template_postgis -O geo geodjango
$ exit
利用 postgres 超级用户来使用 createdb 命令,因为建立PostGIS 数据库需要更高的权利。选项 -O geo 指定 createdb 利用 geo 用户作为数据库的拥有者。当然您也可替换使用其他用户名。
MySQL和Oracle不必建立空间数据库模板,所以您不必在意以上步骤。
利用 django-admin.py 脚本来创建一个 geodjango 项目:
$ django-admin.py startproject geodjango
当 geodjango 项目初始化后,在里面建立 world Django 实例:
$ cd geodjango
$ python manage.py startapp world
geodjango 的项目设置存储在文件 settings.py 中,编辑数据库连接:
DATABASE_ENGINE = 'postgresql_psycopg2'
DATABASE_NAME = 'geodjango'
DATABASE_USER = 'geo'
另外,编辑 INSTALLED_APPS 设置,使之包含 django.contrib.admin 、 django.contrib.gis 、以及 geodjango.world (刚刚新建的项目):
INSTALLED_APPS = (
'django.contrib.auth',
'django.contrib.contenttypes',
'django.contrib.sessions',
'django.contrib.sites',
'django.contrib.admin',
'django.contrib.gis',
'geodjango.world'
)
世界地图数据保存于 zip 压缩文件 . 在 world 实例中创建一个数据目录,下载并解压世界地图数据:
$ mkdir world/data
$ cd world/data
$ wget http://thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip
$ unzip TM_WORLD_BORDERS-0.3.zip
$ cd ../..
这份世界地图zip压缩文件包含了一整套文件,诸如 ESRI Shapefile ,它是众所周知的地理数据格式。解压过后,您会发现还包括:
GDAL的 ogrinfo 是用来检验shapefile或者其他矢量数据的元数据 (metadata)的工具:
$ ogrinfo world/data/TM_WORLD_BORDERS-0.3.shp
INFO: Open of `world/data/TM_WORLD_BORDERS-0.3.shp'
using driver `ESRI Shapefile' successful.
1: TM_WORLD_BORDERS-0.3 (Polygon)
``ogrinfo`` 提示说shapefile只有一个图层,并且该图层含有多边形数据。为了获得更多信息,制定图层名称,并利用选项 ``-so`` 来获取概述信息::
$ ogrinfo -so world/data/TM_WORLD_BORDERS-0.3.shp TM_WORLD_BORDERS-0.3
INFO: Open of `world/data/TM_WORLD_BORDERS-0.3.shp'
using driver `ESRI Shapefile' successful.
Layer name: TM_WORLD_BORDERS-0.3
Geometry: Polygon
Feature Count: 246
Extent: (-180.000000, -90.000000) - (180.000000, 83.623596)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]]
FIPS: String (2.0)
ISO2: String (2.0)
ISO3: String (3.0)
UN: Integer (3.0)
NAME: String (50.0)
AREA: Integer (7.0)
POP2005: Integer (10.0)
REGION: Integer (3.0)
SUBREGION: Integer (3.0)
LON: Real (8.3)
LAT: Real (7.3)
这份详细的概述信息表明图层中地理要素的数量为246,以及其地理几何边界,空间参考系统(“SRS WKT”) 同每一个属性字段的详细信息。例如, FIPS: String (2.0) 表明有一个名为 FIPS 的字符串字段,最大字符长度为2,类似的 LON: Real (8.3) 是一个浮点型字段,最大可容纳8个字节以及最多3个浮点位。尽管这些信息可以很容易的在 世界地图 的网站上找到,但是该命令可以在没有提供元数据时获取这些信息。
现在利用 ogrinfo 来检验世界地图数据,创建一个GeoDjango模型来表达数据:
from django.contrib.gis.db import models
class WorldBorders(models.Model):
# 基于世界地图shapefile文件中的属性字段而构建的Django字段
name = models.CharField(max_length=50)
area = models.IntegerField()
pop2005 = models.IntegerField('Population 2005')
fips = models.CharField('FIPS Code', max_length=2)
iso2 = models.CharField('2 Digit ISO', max_length=2)
iso3 = models.CharField('3 Digit ISO', max_length=3)
un = models.IntegerField('United Nations Code')
region = models.IntegerField('Region Code')
subregion = models.IntegerField('Sub-Region Code')
lon = models.FloatField()
lat = models.FloatField()
# GeoDjango制定的一个地理几何字段 (MultiPolygonField),
# 并用GeoManager实例覆盖默认的管理器。
mpoly = models.MultiPolygonField()
objects = models.GeoManager()
class Meta:
# 制定复数名.
verbose_name_plural = "World Borders"
# 返回表达模型的字符串。
def __unicode__(self):
return self.name
两个需要注意的问题:
用户在模型中声明的地理几何字段,其默认的空间参考系统是WGS84 (其 SRID 是 4326) – 换句话说,经度和纬度是用”度”做单位。如果你想用其他的坐标系统,你可以制定地理几何字段的 srid 。
定义好模型之后,需要同步空间数据库。首先,让我们看看 WorldBorders 模型产生的SQL语句:
$ python manage.py sqlall world
这个命令将输出如下的SQL语句:
BEGIN;
CREATE TABLE "world_worldborders" (
"id" serial NOT NULL PRIMARY KEY,
"name" varchar(50) NOT NULL,
"area" integer NOT NULL,
"pop2005" integer NOT NULL,
"fips" varchar(2) NOT NULL,
"iso2" varchar(2) NOT NULL,
"iso3" varchar(3) NOT NULL,
"un" integer NOT NULL,
"region" integer NOT NULL,
"subregion" integer NOT NULL,
"lon" double precision NOT NULL,
"lat" double precision NOT NULL
)
;
SELECT AddGeometryColumn('world_worldborders', 'mpoly', 4326, 'MULTIPOLYGON', 2);
ALTER TABLE "world_worldborders" ALTER "mpoly" SET NOT NULL;
CREATE INDEX "world_worldborders_mpoly_id" ON "world_worldborders" USING GIST ( "mpoly" GIST_GEOMETRY_OPS );
COMMIT;
如果这段SQL语句符合要求,你可以通过执行 syncdb 命令来建立数据表:
$ python manage.py syncdb
Creating table world_worldborders
Installing custom SQL for world.WorldBorders model
利用 syncdb 命令可以方便的创建一个管理员用户;如果感兴趣,可以参考 createsuperuser 命令。
这一节将向您展示如何利用 LayerMapping 工具将世界地图的shapefile文件导入GeoDajngo模型。导入数据的方式有很多种,除去使用刚才提到的方法,您还可以利用:
前文中谈到利用 ogrinfo 来获取shapefile文件的内容。GeoDjango内置了GDAL强大的OGR库 – 换句话说,您可以利用支持OGR的Python API来处理矢量数据。
首先,开启Django shell:
$ python manage.py shell
如果 世界地图 数据像前文所述的那样下载后,便可以用python的内置模块 os 来确定其路径:
>>> import os
>>> from geodjango import world
>>> world_shp = os.path.abspath(os.path.join(os.path.dirname(world.__file__),
... 'data/TM_WORLD_BORDERS-0.3.shp'))
现在,世界地图shapefile文件会被GeoDjango的 DataSource 接口打开:
>>> from django.contrib.gis.gdal import *
>>> ds = DataSource(world_shp)
>>> print ds
/ ... /geodjango/world/data/TM_WORLD_BORDERS-0.3.shp (ESRI Shapefile)
数据源对象 (Data source) 可以包含多个空间数据图层,然后,shapefile文件只允许有一个图层:
>>> print len(ds)
1
>>> lyr = ds[0]
>>> print lyr
TM_WORLD_BORDERS-0.3
您同样可看到该图层的地理类型以及所包含的地理要数个数:
>>> print lyr.geom_type
Polygon
>>> print len(lyr)
246
Note
shapefile并不允许对地理几何类型的定义。这份shapefile文件和大多数shapefile文件一样实际上包含的是 MultiPolygon 作为其地理几何要素类型。您需要注意的是 GeoDajngo的 PolygonField 字段并不会接受 MultiPolygon 类型的数据 – 这也是在模型中使用 MultiPolygonField 的原因。
Layer 也可以指定空间参考系对象。 – 如果有的话, 则 srs 属性将返回 SpatialReference 对象:
>>> srs = lyr.srs
>>> print srs
GEOGCS["GCS_WGS_1984",
DATUM["WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]]
>>> srs.proj4 # PROJ.4 representation
'+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '
在此,我们注意到shapfile文件使用的是流行的WGS84空间参考系统 – 也就是说利用度作为单位表达经纬度。
另外,shapefile文件中属性数据也可能含有数据。在世界地图图层中,将包含如下属性数据
>>> print lyr.fields
['FIPS', 'ISO2', 'ISO3', 'UN', 'NAME', 'AREA', 'POP2005', 'REGION', 'SUBREGION', 'LON', 'LAT']
检测每个字段的OGR字段类型 (例如,一个字段是整形还是字符串):
>>> [fld.__name__ for fld in lyr.field_types]
['OFTString', 'OFTString', 'OFTString', 'OFTInteger', 'OFTString', 'OFTInteger', 'OFTInteger', 'OFTInteger', 'OFTInteger', 'OFTReal', 'OFTReal']
您可以读取图层中的每一个要素并提取其地理几何字段 (通过 geom 字段获得) 以及其属性字段 (该字段 值 可通过 get() 方法) 读取:
>>> for feat in lyr:
... print feat.get('NAME'), feat.geom.num_points
...
Guernsey 18
Jersey 26
South Georgia South Sandwich Islands 338
Taiwan 363
Layer objects may be sliced:
>>> lyr[0:2]
[<django.contrib.gis.gdal.feature.Feature object at 0x2f47690>, <django.contrib.gis.gdal.feature.Feature object at 0x2f47650>]
每个不同的要素可以通过其要素ID获取:
>>> feat = lyr[234]
>>> print feat.get('NAME')
San Marino
例如,San Marino的地理几何信息被提取并解压为WKT格式或者GeoJSON格式:
>>> geom = feat.geom
>>> print geom.wkt
POLYGON ((12.415798 43.957954,12.450554 ...
>>> print geom.json
{ "type": "Polygon", "coordinates": [ [ [ 12.415798, 43.957954 ], [ 12.450554, 43.979721 ], ...
我们将进一步在 world 实例中创建一个python脚本文件 load.py ,加入如下语句:
import os
from django.contrib.gis.utils import LayerMapping
from models import WorldBorders
world_mapping = {
'fips' : 'FIPS',
'iso2' : 'ISO2',
'iso3' : 'ISO3',
'un' : 'UN',
'name' : 'NAME',
'area' : 'AREA',
'pop2005' : 'POP2005',
'region' : 'REGION',
'subregion' : 'SUBREGION',
'lon' : 'LON',
'lat' : 'LAT',
'mpoly' : 'MULTIPOLYGON',
}
world_shp = os.path.abspath(os.path.join(os.path.dirname(__file__), 'data/TM_WORLD_BORDERS-0.3.shp'))
def run(verbose=True):
lm = LayerMapping(WorldBorders, world_shp, world_mapping,
transform=False, encoding='iso-8859-1')
lm.save(strict=True, verbose=verbose)
解释:
紧接着,从 geodjango 项目目录运行Django shell:
$ python manage.py shell
然后,导入 load 模块,运行 run 让 LayerMapping 做剩下的工作:
>>> from world import load
>>> load.run()
至此,您已经看到如何定义地理模型和如何利用 LayerMapping 工具导入shapefile数据, 其实您还可以进一步利用 ogrinspect 管理命令来自动化您的工作. ogrinspect 命令将支持自省某个GDAL支持的矢量数据源(例如shapefile文件)并自动产生一个模型定义和 LayerMapping 字典数据。
基本使用方法如下:
$ python manage.py ogrinspect <data_source> <model_name> [options]
data_source 指示了GDAL支持的矢量数据的路径,``model_name`` 是模型的名称。options将用于进一步定义模型。
举例而言,如下的命令将自动产生 WorldBorders 模型和先前的映射字典数据:
$ python manage.py ogrinspect world/data/TM_WORLD_BORDERS-0.3.shp WorldBorders --srid=4326 --mapping --multi
解释:
该命令将输出如下的语句,该语句将被复制到GeoDjango应用实例的 models.py 脚本文件中:
# 这是一份由 ogrinspect 自动产生的 Django 模型文件。
from django.contrib.gis.db import models
class WorldBorders(models.Model):
fips = models.CharField(max_length=2)
iso2 = models.CharField(max_length=2)
iso3 = models.CharField(max_length=3)
un = models.IntegerField()
name = models.CharField(max_length=50)
area = models.IntegerField()
pop2005 = models.IntegerField()
region = models.IntegerField()
subregion = models.IntegerField()
lon = models.FloatField()
lat = models.FloatField()
geom = models.MultiPolygonField(srid=4326)
objects = models.GeoManager()
# 为 WorldBorders 模型自动产生的 `LayerMapping` 字典数据
worldborders_mapping = {
'fips' : 'FIPS',
'iso2' : 'ISO2',
'iso3' : 'ISO3',
'un' : 'UN',
'name' : 'NAME',
'area' : 'AREA',
'pop2005' : 'POP2005',
'region' : 'REGION',
'subregion' : 'SUBREGION',
'lon' : 'LON',
'lat' : 'LAT',
'geom' : 'MULTIPOLYGON',
}
GeoDjango 扩展了Django 的ORM模型,并允许使用空间查找。 让我们做一个实验,看看 WorldBorder 模型是否包含某个兴趣点。首先,激活管理shell:
$ python manage.py shell
定义一个兴趣点 [3]:
>>> pnt_wkt = 'POINT(-95.3385 29.7245)'
pnt_wkt 字符串代表一个点在经度 -95.3385 纬度 29.7245. 地理数据类型定义为 Well Known Text (WKT), 它是 Open Geospatial Consortium (OGC)所推荐的地理数据类型。 [4] 导入 WorldBorders 模型, 并利用 pnt_wkt 作为参数执行 contains 查找:
>>> from world.models import WorldBorders
>>> qs = WorldBorders.objects.filter(mpoly__contains=pnt_wkt)
>>> qs
[<WorldBorders: United States>]
这里我们将获得一个 GeoQuerySet 对象,该对象只包含一个要素,这也正是我们所需要寻找的。 同样的,一个 GEOS 地理对象 也可以能被使用 – 比如用 intersects 空间查找来获取 WorldBorders 对象的实例 San Marino ,该空间查找是和 get() 方法一同使用:
>>> from django.contrib.gis.geos import Point
>>> pnt = Point(12.4604, 43.9420)
>>> sm = WorldBorders.objects.get(mpoly__intersects=pnt)
>>> sm
<WorldBorders: San Marino>
空间查找 contains 和 intersects 属于 GeoDjango 数据库接口 ,查看文档获得更多详细资料。
当查询空间数据时,如果地理数据拥有不同的坐标投影,GeoDjango将自动转换。在接下来的的例子中, 坐标系统将表达为 EPSG SRID 32140 , 地图投影坐标系统 仅 在南德州使用,其单位为米而非度。:
>>> from django.contrib.gis.geos import *
>>> pnt = Point(954158.1, 4215137.1, srid=32140)
注意 pnt 对象也可以使用EWKT进行构建,”extended” 代表WKT可以包含SRID:
>>> pnt = GEOSGeometry('SRID=32140;POINT(954158.1 4215137.1)')
当使用 GeoDjango 的 ORM 模型时,将自动把地理要素值转化为SQL,这样可以使开发者在更高层次上进行工作:
>>> qs = WorldBorders.objects.filter(mpoly__intersects=pnt)
>>> qs.query.as_sql() # Generating the SQL
('SELECT "world_worldborders"."id", "world_worldborders"."name", "world_worldborders"."area",
"world_worldborders"."pop2005", "world_worldborders"."fips", "world_worldborders"."iso2",
"world_worldborders"."iso3", "world_worldborders"."un", "world_worldborders"."region",
"world_worldborders"."subregion", "world_worldborders"."lon", "world_worldborders"."lat",
"world_worldborders"."mpoly" FROM "world_worldborders"
WHERE ST_Intersects("world_worldborders"."mpoly", ST_Transform(%s, 4326))',
(<django.contrib.gis.db.backend.postgis.adaptor.PostGISAdaptor object at 0x25641b0>,))
>>> qs # printing evaluates the queryset
[<WorldBorders: United States>]
GeoDjango可以将对地理数据进行标准化的表达。除了访问地理字段,GeoDjango可以创建 GEOS 地理对象 , 并给地理对象强大的功能,例如序列化地理数据:
>>> sm = WorldBorders.objects.get(name='San Marino')
>>> sm.mpoly
<MultiPolygon object at 0x24c6798>
>>> sm.mpoly.wkt # WKT
MULTIPOLYGON (((12.4157980000000006 43.9579540000000009, 12.4505540000000003 43.9797209999999978, ...
>>> sm.mpoly.wkb # WKB (as Python binary buffer)
<read-only buffer for 0x1fe2c70, size -1, offset 0 at 0x2564c40>
>>> sm.mpoly.geojson # GeoJSON (requires GDAL)
'{ "type": "MultiPolygon", "coordinates": [ [ [ [ 12.415798, 43.957954 ], [ 12.450554, 43.979721 ], ...
包括所有GEOS库提供的高级地理几何操作:
>>> pnt = Point(12.4604, 43.9420)
>>> sm.mpoly.contains(pnt)
True
>>> pnt.contains(sm.mpoly)
False
GeoDjango 提供 Django Admin 以使用户通过 OpenLayers 构建并修改地理要素。
让我们再进一步,在 world 应用实例中创建 admin.py 文件,加入如下代码:
from django.contrib.gis import admin
from models import WorldBorders
admin.site.register(WorldBorders, admin.GeoModelAdmin)
然后,编辑 geodjango 项目目录下的 urls.py 文件:
from django.conf.urls.defaults import *
from django.contrib.gis import admin
admin.autodiscover()
urlpatterns = patterns('',
(r'^admin/(.*)', admin.site.root),
)
启动 Django 开发服务器:
$ python manage.py runserver
最后,您能通过 http://localhost:8000/admin/ 进行浏览,通过用 sync 所生成的管理用户名进行登录。您可以浏览任何一个 WorldBorders 实例 – 您可以点选任何国家的边界甚至拖动节点。
通过 OSMGeoAdmin ,GeoDjango 将在管理模块中使用 Open Street Map 图层。这将提供更多的内容 (包括道路和街道详情) 而不是 OpenLayers 默认提供的 Vector Map Level 0 WMS服务 (由 Metacarta 支持)。
请注意:
万事俱备则在 admin.py 文件中用如下代码替换 OSMGeoAdmin:
admin.site.register(WorldBorders, admin.OSMGeoAdmin)
脚注
| [1] | Special thanks to Bjørn Sandvik of thematicmapping.org for providing and maintaining this data set. |
| [2] | GeoDjango basic apps was written by Dane Springmeyer, Josh Livni, and Christopher Schmidt. |
| [3] | Here the point is for the University of Houston Law Center . |
| [4] | Open Geospatial Consortium, Inc., OpenGIS Simple Feature Specification For SQL, Document 99-049. |