CopyData.py 5.6 KB
# coding=utf-8
#author:        4N
#createtime:    2021/6/11
#email:         nheweijun@sina.com
import copy
from osgeo.ogr import *
from osgeo import gdal,ogr
import uuid
import time
import os
import json


def get_info_from_sqlachemy_uri(uri):
    parts = uri.split(":")
    user = parts[1][2:]

    password_list = parts[2].split("@")
    if password_list.__len__() > 2:
        password = "@".join(password_list[:-1])
    else:
        password = parts[2].split("@")[0]
    host = parts[2].split("@")[-1]
    port = parts[3].split("/")[0]
    database = parts[3].split("/")[1]

    return user, password, host, port, database

def open_pg_data_source(iswrite, uri):
    """
    # 获取PostGIS数据源
    :return:
    """
    db_conn_tuple = get_info_from_sqlachemy_uri(uri)
    fn = "PG: user=%s password=%s host=%s port=%s dbname=%s " % db_conn_tuple
    driver = ogr.GetDriverByName("PostgreSQL")
    if driver is None:
        raise Exception("打开PostgreSQL驱动失败,可能是当前GDAL未支持PostgreSQL驱动!")
    ds = driver.Open(fn, iswrite)
    if ds is None:
        raise Exception("打开数据源失败!")
    return ds


def move_digui(coords,offx,offy):
    if isinstance(coords[0],list):
        for coord in coords:
            move_digui(coord, offx, offy)
    else :
        coords[0] += offx
        coords[1] += offy

def move(geo,offx,offy):


    geojson_str = geo.ExportToJson()
    geojson = json.loads(geojson_str)

    coords = geojson["coordinates"]
    move_digui(coords, offx, offy)
    og = ogr.CreateGeometryFromJson(json.dumps(geojson))

    return og


def copydata(bound,input,output):
    work_dir = os.path.dirname(os.path.abspath(__file__))

    shp_driver: Driver = ogr.GetDriverByName("ESRI Shapefile")
    bound_ds : DataSource = shp_driver.Open(bound, 0)
    bound_layer :Layer = bound_ds.GetLayer(0)

    bf = bound_layer.GetNextFeature()
    bound_geom = bf.GetGeometryRef()



    gdb_driver: Driver = ogr.GetDriverByName("FileGDB")
    # gdb_ds: DataSource = gdb_driver.Open(gdb, 0)
    # layer: Layer = gdb_ds.GetLayerByName(layer_name)

    input_ds: DataSource = shp_driver.Open(input, 0)
    layer: Layer = input_ds.GetLayer(0)

    if not layer:
        raise Exception("打开数据失败!")


    output_ds: DataSource = gdb_driver.CreateDataSource(output)


    gdb_layer:Layer = output_ds.CreateLayer(layer.GetName(),layer.GetSpatialRef(),layer.GetGeomType())

    gdb_layer.CreateFields(layer.schema)




    #
    bound_exent = bound_layer.GetExtent()

    layer_extent = layer.GetExtent()

    xxrange=layer_extent[1]-layer_extent[0]
    yrange = layer_extent[3]-layer_extent[2]


    left = (int((layer_extent[0] - bound_exent[0]) / xxrange) + 1) * (-1)
    right = int((bound_exent[1] - layer_extent[1]) / xxrange) + 1
    up = int((bound_exent[3] - layer_extent[3]) / yrange) + 1
    down =  (int((layer_extent[2] - bound_exent[2]) / yrange) + 1) * (-1)


    begin = time.time()
    count=1


    for f in layer:

        new_f = copy.copy(f)
        for xt in range(left,right+1,1):
            for yt in range(down,up+1,1):
                out_g = move(f.GetGeometryRef(),xxrange*xt,yrange*yt)

                if out_g.Intersect(bound_geom):
                    out_g = out_g.Intersection(bound_geom)
                    new_f.SetGeometry(out_g)
                    new_f.SetFID(count)
                    gdb_layer.CreateFeature(new_f)
                    count += 1
                    if count % 10000 == 0:
                        # print(count)
                        with open(os.path.join(work_dir, "copy.txt"), "w") as fi:
                            fi.write("已完成{}".format(count))

    print(time.time()-begin)

    input_ds.Destroy()

def envelop_2_polygon(env):
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(env[0], env[2])
    ring.AddPoint(env[0], env[3])
    ring.AddPoint(env[1], env[3])
    ring.AddPoint(env[1], env[2])
    ring.AddPoint(env[0], env[2])
    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    return poly

def clip():

    input = r"E:\Data\copy\x.shp"
    output = r"E:\Data\copy\xianorigin.shp"
    gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
    shp_driver: Driver = ogr.GetDriverByName("ESRI Shapefile")

    input_ds: DataSource = shp_driver.Open(input, 0)
    layer: Layer = input_ds.GetLayer(0)

    bound_geo = envelop_2_polygon([112.84,112.952,22.914,22.985])

    if not layer:
        raise Exception("打开数据失败!")

    output_ds: DataSource = shp_driver.CreateDataSource(output)
    #
    # out_layer: Layer = output_ds.CreateLayer(layer.GetName(), layer.GetSpatialRef(), layer.GetGeomType())
    # out_layer.CreateFields(layer.schema)
    for f in layer:
        new_f:Feature = copy.copy(f)
        out_g = f.GetGeometryRef()
        if out_g.Intersect(bound_geo):
            out_g = out_g.Intersection(bound_geo)
            new_f.SetGeometry(out_g)
            new_f.SetField("ttt",1)
            # new_f.UnsetField("OBJECTID")
            # out_layer.CreateFeature(new_f)

    # output_ds.Destroy()



if __name__ == '__main__':
    # bound_shp = "E:\Data\广东边界\gdsample.shp"
    # input = r"E:\Data\copy\origin.shp"
    # output = r"E:\Data\copy\re.shp"
    # copydata(bound_shp,input,output)
    #
    # bound_shp = "/root/CopyData/gdsample.shp"
    # input = r"/root/CopyData/origin.shp"
    # output = r"/root/CopyData/re.gdb"
    # copydata(bound_shp,input,output)
    clip()