「GIS教程」Python利用WKT创建shapefile、shapefile输出WKT

最近在做一个网格化管理系统时，本准备使用PostgreSQL数据库来存储网格WKT数据（方便直接导出shp文件，定期发布地图服务），奈何它不能很好地运行后端框架的表结构，遂放弃选择继续使用MySQL。那用MySQL导出shp便成了难题，经过我的一番操作成功解决该问题：利用Java定时任务，调用Python脚本，遍历MySQL网格化数据表，将WKT转为shapefile要素类，废话不多说给大家贴代码瞅
下。

WKT创建shapefile

由于单位数据不能公开，这里就用个表格简单模拟下。首先是读取表格数据，创建要素类和属性表结构，其次写入投影信息，最后创建要素，关闭数据流即可。

效果预览

实现代码

import osr
import datetime
import pandas as pd
from osgeo import ogr

local_file = r'./results/wkt20211201.xlsx'; # 输入数据源
# 输出要素类
out_shapefile = r"./results/shp/result_{}.shp".format(datetime.datetime.now().strftime('%Y%m%d') )
# 创建要素类
outDriver = ogr.GetDriverByName("ESRI Shapefile")
outDataSource = outDriver.CreateDataSource(out_shapefile)
# 设置投影
srs = osr.SpatialReference()
# 4326-GCS_WGS_1984; 4490-GCS_China_Geodetic_Coordinate_System_2000 srs.ImportFromEPSG(4326)
outLayer = outDataSource.CreateLayer("states_extent", srs, geom_type=ogr.wkbPolygon)

# 创建3个字段并设置类型和长度
aField = ogr.FieldDefn("name", ogr.OFTString) # 设置字段名次与类型 aField.SetWidth(20) # 设置字段长度
outLayer.CreateField(aField) # 创建字段
bField = ogr.FieldDefn("adcode", ogr.OFTString) bField.SetWidth(100)
outLayer.CreateField(bField)
cField = ogr.FieldDefn("wkt", ogr.OFTString)
# 字符串最大字段254
cField.SetWidth(254)
outLayer.CreateField(cField)
# 写记录函数
def addPoly(param_a, param_b, param_c, wktpoly):
# 创建几何
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn) feature.SetField("name",str(param_a).encode("gbk").decode('ISO-8859-1')) #
编码转换
    feature.SetField("adcode", str(param_b))
    feature.SetField("wkt", param_c)
    feature.SetGeometry(wktpoly)
    outLayer.CreateFeature(feature)
    feature = None
# 读取WKT数据表
columns_list = pd.read_excel(local_file) # 写入记录至要素类
for i in range(len(columns_list)):
field_name = columns_list['name'][i] field_adcode = columns_list['adcode'][i] field_wkt = columns_list['WKT'][i]
# WKT转为图形
    wktpoly = ogr.CreateGeometryFromWkt(field_wkt)
    addPoly(field_name, field_adcode, field_wkt, wktpoly)
outDataSource = None

shapefile导出WKT表格

遍历要素类的中属性字段，选择输出特定的字段及WKT文本，存储到excel表格中。

输入数据

实现代码

import datetime
import pandas as pd
from osgeo import ogr
# 读取数据
shapefile = r"./data-use/shp/广东省.shp"
out_file = r'./results/wkt' + str(datetime.datetime.now().strftime('%Y%m%d')) + '.xlsx' # 输出表格
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
layerDefinition = layer.GetLayerDefn()
# 获取字段名称列表
fields = []
for i in range(layerDefinition.GetFieldCount()):
    fields.append(layerDefinition.GetFieldDefn(i).GetName())
print(fields)
# 获取字段值和相应的几何 record = []
for feature in layer:
    lst = []
    name = feature.GetField("NAME")
    lst.append(name)
    adcode = feature.GetField("adcode")
    lst.append(adcode)
        geom = feature.GetGeometryRef()
    geomwkt = geom.ExportToWkt()
    lst.append(geomwkt)
    record.append(lst)
# 写入excel中
df = pd.DataFrame(record, columns=['name', 'adcode', 'WKT']) df.to_excel(out_file, encoding='gbk')

表格数据转shp

类似ArcMap中的x、y转点，不过利用Python可以进行批处理操作，我一般将爬取的POI转为矢量 数据，发布地图服务，在网页中调用;也可以通过geopandas库调用在线底图，叠加POI矢量shp形成 下图效果。

实现代码

import shapefile # 使用pyshp import pandas as pd
# 表格文件存放位置
local_file = r'./data-use/table/data_poi_hang-zhou-shi-ATM.csv' # 新建数据存放位
result_address = './results/shp/test_01.shp'
df = pd.read_csv(local_file)
# 开始操作数据
file = shapefile.Writer(result_address)
# 创建字段
columns_list = list(df.head(0)) # 获取字段列名 print(columns_list)
# 添加字段(较长的字段放宽长度)
for i in range(len(columns_list)):
        if columns_list[i] in ['business_area', 'cityname', 'id', 'name', 'pname']: file.field(columns_list[i], 'C', '100') # C代表数据类型为字符串，长度为100
    elif columns_list[i] in ['lat', 'lon']:
            file.field(columns_list[i], 'N', '31', decimal=10) # 数值类型double # 添加属性表
for num in range(len(df['lon'])):
    file.point(df['lon'][num], df['lat'][num])
    # address,business_area,cityname,id,lat,lon,name,pname,tel,type
    file.record(df['business_area'][num], df['cityname'][num], df['id'][num],df['lat'][num],df['lon'][num], df['name'][num], df['pname'][num])
file.close()

# 定义投影方法1(通过epsg来定义)
# proj = osr.SpatialReference()
# proj.ImportFromEPSG(4326) # 4326-GCS_WGS_1984; 4490- GCS_China_Geodetic_Coordinate_System_2000
# proj.ExportToWkt()
# 写入投影
f = open(result_address.replace(".shp", ".prj"), 'w')
# f.write(proj)
# 定义投影方法2
# 如果是别的修改WKT代码即可
WGS84_WKT = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.2572 23563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]' f.write(WGS84_WKT)
f.close()

写在最后

以上就是全部内容 ，推荐大家使用miniconda来管理Python环境，添加国内镜像源，下载速度很 快，简直不要太好用。全部代码和示例数据已上传至码云仓库，欢迎大家下载查看:https://gitee.co m/fungiser/python-shapefile-operate