利用Gdal的Postgis数据地图导出思路

环境:

 windows vista

 python2.6

 psycopg2

 PostgreSQL8.3

 POstGIS1.4.0

 GDAL-1.6.1-py2.6-win32.egg

 GDAL-1.6.1-py2.6-win32.exe

 

 首先利用Postgis的空间函数切割地图数据：参考：《postgis的地图切割方案》

 数据库部分表结构：

 CREATE TABLE zhandian<br> (<br> city_id character varying(64) NOT NULL,<br> objectid character varying(64) NOT NULL,<br> user_label character varying(255) NOT NULL,<br> the_geom geometry,<br> CONSTRAINT enforce_dims_the_geom CHECK (ndims(the_geom) = 2),<br> CONSTRAINT enforce_geotype_the_geom CHECK (geometrytype(the_geom) = ‘POINT’::text OR the_geom IS NULL),<br> CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = (-1))<br> )<br> WITH (OIDS=FALSE);<br> ALTER TABLE zhandian OWNER TO postgres;

 利用psycopg2取出相应的数据，然手写入shp文件：

 以下是单行记录导入shp文件的演示：(transfer.py）

from osgeo import ogr<br> import psycopg2<br> import sys

ogr.RegisterAll()

#加载驱动，定义生成的文件,如ESRI Shapefile格式的zhandian.shp
driver=ogr.GetDriverByName(‘ESRI Shapefile’)

ds=driver.CreateDataSource(‘E:/pyWorkspace/zhandian.shp’)

layer=ds.CreateLayer(‘zhandian’,geom_type=ogr.wkbPoint)

#定义字段
field1=ogr.FieldDefn(‘User_label’,ogr.OFTString)
field1.SetWidth(100)
field2=ogr.FieldDefn(‘objectid’,ogr.OFTString)
field2.SetWidth(64)

layer.CreateField(field1)
layer.CreateField(field2)

featureDefn=layer.GetLayerDefn()

feature=ogr.Feature(featureDefn)

#数据库查询
conn=psycopg2.connect(‘host=localhost port=5432 password=postgres dbname=gisdb_jl user=postgres’)
cur=conn.cursor()
sql=”SELECT user_label,objectid,ST_ASGML(the_geom) FROM zhandian;”
cur.execute(sql)
result=cur.fetchone()
conn.commit()
cur.close()
conn.close()

#组织数据
feature.SetField(0,result[0])
feature.SetField(1,result[1])
point=ogr.CreateGeometryFromGML(result[2]);
feature.SetGeometry(point)

layer.CreateFeature(feature)

#释放内存
feature.Destroy()
ds.Destroy()

  

为解决中文问题，在dos下 先运行：

set PGCLIENTENCODING=GBK

然后python transfer.py

取出的数据：

 最后的shp效果图：

 

转载自：https://blog.csdn.net/no_cross_no_crown/article/details/6268970