#!/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)

类库参考:https://gdal.org/python/osgeo.ogr-module.html

Logo

CSDN联合极客时间,共同打造面向开发者的精品内容学习社区,助力成长!

更多推荐