WKT创建shapefile、shapefile输出WKT

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.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'
f.write(WGS84_WKT)
f.close()

写在最后

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