# -*- coding: utf-8 -*-

"""
@File    : 17_面矢量数据叠加分析.py
@Author  : fungis@163.com
@notice  :
"""

from osgeo import ogr, osr
from osgeo import gdal
import os


def calculated_proportion(in_path1, in_path2, out_path):
    """
    面图层相交，属性字段保留，两个shp文件必须是同一个坐标系
    :param in_path1: 主shp文件
    :param in_path2: 用来相交的shp文件
    :param out_path: 输出结果路径
    :return:
    """
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.Open(in_path1, 1)
    layer_1 = data_source.GetLayer()
    crs = layer_1.GetSpatialRef()
    epsg1 = crs.GetAttrValue('AUTHORITY', 1)

    r_data_source = driver.Open(in_path2, 1)
    layer_2 = r_data_source.GetLayer()
    crs = layer_2.GetSpatialRef()
    epsg2 = crs.GetAttrValue('AUTHORITY', 1)
    if epsg1 != epsg2:
        print("坐标系不一致")
    else:
        print("坐标系一致")

    # 新建DataSource，Layer
    out_ds = driver.CreateDataSource(out_path)
    out_lyr = out_ds.CreateLayer("out", layer_1.GetSpatialRef(), layer_1.GetGeomType())
    # 遍历原始的Shapefile文件给每个Geometry做Buffer操作
    # current_union = layer[0].Clone()
    print('the length of layer:', len(layer_1))
    if len(layer_1) == 0:
        return

    p = ["SKIP_FAILURES=YES", "PROMOTE_TO_MULTI=YES", "INPUT_PREFIX=1", "METHOD_PREFIX=2"]
    # gdal.TermProgress_nocb 控制台输出进度
    layer_1.Intersection(layer_2, out_lyr, p, gdal.TermProgress_nocb)
    out_lyr.SyncToDisk()
    del data_source, r_data_source, out_ds


if __name__ == '__main__':
    shp_path = r"E:\gisData\test2\绍兴市.shp"
    shp_path2 = r"E:\gisData\test2\浙江省.shp"
    out_path = r"E:\gisData\test2\out.shp"
    calculated_proportion(shp_path, shp_path2, out_path)