python ogr读写shp、dbf、转换为geojson工具
#!/usr/bin/python# coding=utf-8import gdalimport osimport sysfrom osgeo import ogrimport time# 当前工作目录# path = os.getcwd()# print("\n工作地址: ", path)# outpath = str(path) + "/LD_SRC_DATA_MESH_JSON"# prin
·
#!/usr/bin/python
# coding=utf-8
import gdal
import os
import sys
from osgeo import ogr
import time
# 当前工作目录
# path = os.getcwd()
# print("\n工作地址: ", path)
# outpath = str(path) + "/LD_SRC_DATA_MESH_JSON"
# print("\n输出地址: ", outpath)
# 父目录
# parent = os.path.join(path, os.pardir)
# 源数据处理排除文件名称
source_shp_name_filter = ["Node"]
# 源数据处理类型增加文件
source_shp_type_filter = ["poi.dbf"]
# 源数据属性需转换字典
source_shp_properties_filer = ["XMIN", "YMIN", "XMAX", "YMAX", "X_COORD", "Y_COORD", "X_ENTR", "Y_ENTR"]
# 源数据位置
file_path = "D:/vector-tile-program/16Q38M_LD_Src_Data_MESH"
out_path = "D:/vector-tile-program/16Q38M_LD_Src_Data_MESH_OUT"
# 父目录
# print("\n父目录:", os.path.abspath(parent))
def do_layerPoint(layer):
ftr = layer.ResetReading()
ftr = layer.GetNextFeature()
# print('point num:', layer.GetFeatureCount())
# print('extent:', layer.GetExtent())
cc = 1
while ftr:
# print cc
cc += 1
pt = ftr.GetGeometryRef().GetPoint(0)
g = ftr.GetGeometryRef()
# print g#,g.ExportKML()
if pt[0] > 1000 or pt[1] > 1000:
g.SetPoint(0, round(pt[0] / 3600.0, 6), round(pt[1] / 3600.0, 6))
# print(pt)
# print(g)
'''
ng = ogr.Geometry(ogr.wkbPoint)
print pt
ng.SetPoint(0,pt[0]+40,pt[1])
ftr.SetGeometry(ng)
'''
layer.SetFeature(ftr)
ftr = layer.GetNextFeature()
def do_layerLine(layer):
ftr = layer.ResetReading()
ftr = layer.GetNextFeature()
while ftr:
g = ftr.GetGeometryRef()
cnt = g.GetPointCount()
cc = 0
while cc < cnt:
# print(g.GetPoint(cc))
cc += 1
for n in range(cnt):
pt = g.GetPoint(n)
if pt[0] > 1000 or pt[1] > 1000:
g.SetPoint(n, round(pt[0] / 3600.0, 6), round(pt[1] / 3600.0, 6))
layer.SetFeature(ftr)
ftr = layer.GetNextFeature()
def merge_dbf(layer, file):
# 清除历史
merge_dbf_buff = {}
# 获取poiname
path = file.split("/")
path.remove("poi.dbf")
path = "/".join(path) + "/PoiName.dbf"
# print("\n" + path)
driver = ogr.GetDriverByName('ESRI Shapefile')
inds = driver.Open(path, 0) # 0 - read , 1 - write
poi_name_layer = inds.GetLayer()
if poi_name_layer.GetFeatureCount() == 0:
return
# print("name:" + poi_name_layer.GetName())
layer_defn = poi_name_layer.GetLayerDefn()
feature_count = layer_defn.GetFieldCount()
# 获得poiname字典
for i in range(0, layer.GetFeatureCount()):
feature = poi_name_layer.GetFeature(i)
# print(feature)
fields = {}
for feat in range(feature_count):
field = layer_defn.GetFieldDefn(feat)
field_index = layer_defn.GetFieldIndex(field.GetNameRef())
# print(field.GetNameRef(), feature.GetField(field_index))
fields[field.GetNameRef()] = feature.GetField(field_index)
merge_dbf_buff[i] = fields
# print(merge_dbf_buff)
# 添加到poi
if layer.GetFeatureCount() == 0:
return
# print("name:" + layer.GetName())
layer_defn_l = layer.GetLayerDefn()
feature_count_l = layer_defn_l.GetFieldCount()
# 添加poi name
for i in range(0, layer.GetFeatureCount()):
feature = layer.GetFeature(i)
# print(feature)
for bu in merge_dbf_buff:
if feature.GetField("GUID") == merge_dbf_buff[bu]["GUID"]:
# print(feature.GetField("GUID"), merge_dbf_buff[bu]["GUID"])
for gu in merge_dbf_buff[bu]:
if gu != "GUID":
field_index = layer_defn_l.GetFieldIndex(gu)
if field_index < 0:
# print(bu, gu)
field_defn = ogr.FieldDefn(gu, ogr.OFTString)
layer.CreateField(field_defn, 1)
# print(merge_dbf_buff)
# 添加poi name 值
for i in range(0, layer.GetFeatureCount()):
feature = layer.GetFeature(i)
# print(feature)
for bu in merge_dbf_buff:
if feature.GetField("GUID") == merge_dbf_buff[bu]["GUID"]:
for gu in merge_dbf_buff[bu]:
if gu != "GUID":
fieldIndex = layer_defn_l.GetFieldIndex(gu)
# print(fieldIndex)
if fieldIndex > 0:
sss = merge_dbf_buff[bu][gu]
feature.SetField(fieldIndex, merge_dbf_buff[bu][gu])
layer.SetFeature(feature)
# print(gu, sss)
def do_layerDbf_and_setProperties(layer):
# 属性坐标秒度转换
layer_defn = layer.GetLayerDefn()
feature_count = layer_defn.GetFieldCount()
# fields = []
for feat in range(feature_count):
field = layer_defn.GetFieldDefn(feat)
# fields.append(field.GetNameRef())
# log(str(fields))
# for a in fields:
name = field.GetNameRef()
# source_shp_properties_filer
for pro in source_shp_properties_filer:
if name != pro:
continue
# 获取所需字段的序号index
# print(pro)
field_index = layer_defn.GetFieldIndex(pro)
# print(fieldIndex)
for i in range(0, layer.GetFeatureCount()):
feature = layer.GetFeature(i)
# 新的属性内容
code = feature.GetField(field_index)
if code > 1000:
code = code / 3600.0
feature.SetField(field_index, round(code, 6))
layer.SetFeature(feature)
# 合并poi+poiname
# print("\n" + layer.GetName())
def do_layerPolygon(layer):
ftr = layer.ResetReading()
ftr = layer.GetNextFeature()
while ftr:
g = ftr.GetGeometryRef()
cnt = g.GetGeometryCount()
for n in range(cnt):
gg = g.GetGeometryRef(n)
for m in range(gg.GetPointCount()):
pt = gg.GetPoint(m)
# print(pt)
if pt[0] > 1000 or pt[1] > 1000:
gg.SetPoint(m, round(pt[0] / 3600.0, 6), round(pt[1] / 3600.0), 6)
layer.SetFeature(ftr)
ftr = layer.GetNextFeature()
def do_shpfile(file):
# print file
# print('ready file:', file)
layer_name = file.split("/")[-1].split(".")[0]
if shp_filter(layer_name) is False:
return
driver = ogr.GetDriverByName('ESRI Shapefile')
# shp = driver.Open('e:/shp_data/points.shp',1) # 0 - read , 1 - write
inds = driver.Open(file, 1) # 0 - read , 1 - write
layer = inds.GetLayer()
if layer.GetFeatureCount() == 0:
return
gtype = layer.GetLayerDefn().GetGeomType()
if file.lower().find('province') == -1:
pass # return
# 统一属性秒度转换
do_layerDbf_and_setProperties(layer)
if gtype == ogr.wkbPoint:
# print("=== type:wkbPoint ===")
do_layerPoint(layer)
elif gtype == ogr.wkbLineString:
# print("=== type:wkbLineString ===")
do_layerLine(layer)
elif gtype == ogr.wkbPolygon:
# print("=== type:wkbPolygon ===")
do_layerPolygon(layer)
elif gtype == 100:
merge_dbf(layer, file)
else:
log('unknown type:' + str(file))
inds.SyncToDisk()
inds.Destroy()
# 转换到geojson
to_geojson(file, file.split(layer_name)[0] + "/" + layer_name + ".geojson")
def convert(shpdir):
files = os.listdir(shpdir)
for file in files:
# 普通shp文件处理
if file.lower().find('.shp') > -1:
# print(file.lower().find('.shp'), shpdir)
# if file == 'points.shp':
do_shpfile(shpdir + "/" + file)
# 加入dbf文件处理
for t in source_shp_type_filter:
if file.lower().find(t) > -1:
# print(file.lower().find(t), shpdir + "/" + t)
do_shpfile(shpdir + "/" + t)
def getFatherDir(path):
return os.path.join(path, os.pardir)
def create_dir():
isExists = os.path.exists(out_path)
if not isExists:
os.mkdir(out_path)
# 检测排除不需要转换的shp
def shp_filter(name):
boo = True
for filename in source_shp_name_filter:
if filename == name:
boo = False
return boo
def to_geojson(vector, output):
# 打开矢量图层
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
shp_ds = ogr.Open(vector)
shp_lyr = shp_ds.GetLayer(0)
# 创建结果Geojson
base_name = os.path.basename(output)
out_driver = ogr.GetDriverByName('GeoJSON')
out_ds = out_driver.CreateDataSource(output)
if out_ds.GetLayer(base_name):
out_ds.DeleteLayer(base_name)
out_lyr = out_ds.CreateLayer(base_name, shp_lyr.GetSpatialRef())
out_lyr.CreateFields(shp_lyr.schema)
out_feat = ogr.Feature(out_lyr.GetLayerDefn())
# 生成结果文件
for feature in shp_lyr:
out_feat.SetGeometry(feature.geometry())
for j in range(feature.GetFieldCount()):
out_feat.SetField(j, feature.GetField(j))
out_lyr.CreateFeature(out_feat)
del out_ds
del shp_ds
# print("to geojson success")
# 进度
def process_bar(percent, start_str='', end_str='', total_length=0):
bar = ''.join(["▋"] * int(percent * total_length)) + ''
bar = '\r' + start_str + bar.ljust(total_length) + ' {:0>4.1f}%|'.format(percent * 100) + end_str
print(bar, end='', flush=True)
def log(text_log):
text_log = "\n" + time.strftime("%Y-%m-%d %H:%M:%S", time.localtime()) + " " + text_log + "."
print(text_log, "")
def get_dir_list(file_dir):
# 通过 listdir 得到的是仅当前路径下的文件名,不包括子目录中的文件,如果需要得到所有文件需要递归
return os.listdir(file_dir)
def do_vector_tile():
pass
def do_shp_to_json(shp_file_path):
index = 0
for p in shp_file_path:
index += 1
# 处理开始
convert(p)
# 进度显示
# log("progress " + str(round(100 / len(shp_file_path), 2) * index) + "%,shptojson")
process_bar(round(100 / len(shp_file_path), 2) * index / 100, start_str='转换 ', end_str=' 100% ',
total_length=100)
if __name__ == '__main__':
# convert( 'e:/shp_data' )
print("------------------------------------")
print("----------- 3 秒后作业开始-----------")
print("------------------------------------")
time.sleep(3)
# 创建导出文件夹
# create_dir()
log("工作目录: " + file_path)
log("输出地址: " + out_path)
dir_list = get_dir_list(file_path)
# 字典记录当前目录
dir_dictionary = {}
for d in dir_list:
dir_arr = get_dir_list(file_path + "/" + d)
dir_dictionary[d] = dir_arr
log("作业目录:" + str(dir_dictionary))
# 作业开始
shp_file_path = []
for d in dir_dictionary:
for dd in dir_dictionary[d]:
shp_file_path.append(file_path + "/" + d + "/" + dd)
log("作业开始:")
if sys.argv[1:]:
convert(sys.argv[1])
else:
# convert('D:/vector-tile-program/16Q38M_LD_Src_Data_MESH/N/N52F048004')
# print(get_dir_list("D:/vector-tile-program/16Q38M_LD_Src_Data_MESH"))
# 执行批量转换
do_shp_to_json(shp_file_path)
# 执行sh切片开始
# do_vector_tile(shp_file_path)
更多推荐
已为社区贡献1条内容
所有评论(0)