導航

import shapely.geometry as geometry

from shapely.geometry import Polygon
import numpy as np
# 數據格式轉換
from osgeo import gdal_array, ogr
from numpy import *

# 提取邊界點
def eage_polygon(wkt):
    strwkt = str(wkt)
    pixels = []
    if strwkt[0] == "M":
        point = strwkt[14:-2].split(",")
        for k in range(0, len(point)):
            zb = point[k].split(" ")

            if zb[0][0] == "(" and zb[0][1] == "(":
                pixels.append([ float(zb[0][2:]), float(zb[1])])
            elif zb[0][0] == "(" and zb[0][1] != "(":
                pixels.append( [float(zb[0][1:]), float(zb[1])])
            elif zb[1][-1] == ")" and zb[1][-2] != ")":
                pixels.append( [float(zb[0]), float(zb[1][:-1])])
            elif zb[1][-1] == ")" and zb[1][-2] == ")":
                pixels.append( [float(zb[0]), float(zb[1][:-2])])
            else:
                pixels.append( [float(zb[0]), float(zb[1])])
    if strwkt[0] == "P":
        point = strwkt[9:-1].split(",")

        for k in range(0, len(point)):
            zb = point[k].split(" ")

            if zb[0][0] == "(":
                pixels.append( (float(zb[0][1:]), float(zb[1])))

            elif zb[1][-1] == ")":
                pixels.append(( float(zb[0]), float(zb[1][:-1])))
            else:
                pixels.append( (float(zb[0]), float(zb[1])))
    return pixels

# 判斷點是否在多邊形內,并剔除粗差點
def isinpolygon(flat_points,geometry_polygon):
    in_shape_points=[]
    geo=eage_polygon(geometry_polygon)
    for pt in flat_points:
        # make a point and see if it's in the polygon
        if geometry.Point(pt[:2]).within(Polygon(geo)):
            in_shape_points.append(pt[2])

    '''
    四分法粗差點剔除
    '''
    in_shape_points=sort(in_shape_points)
    n=len(in_shape_points)
    q1=in_shape_points[int((n-1)/4+1)]
    q3 = in_shape_points[int(3*(n - 1) / 4 + 1)]
    vmin=q1-1.5*(q3-q1)
    vmax=q3+1.5*(q3-q1)
    in_shape_points=np.array(in_shape_points)
    aaa=in_shape_points[in_shape_points>vmin]
    bbb=aaa[aaa<vmax]
    return bbb


# 按半徑裁剪散點
def radius_cut(point_zone,R,xyzData):
    dist = xyzData - np.array([point_zone[0], point_zone[1], 0])
    x = (dist.T[0] * dist.T[0] + dist.T[1] * dist.T[1]).T
    return xyzData[np.where(x < R)]

# 按shp圖版統計散點信息
def part_statistic_message(shp_path, point_txt):
    #讀取點云
    xyzData = []
    with open(point_txt, "r") as f:
        for line in f.readlines():
            line = line.strip('\n')  # 去掉列表中每一個元素的換行符
            a = line.split(",")
            xyzData.append([float(a[0]), float(a[1]), float(a[2])])
    xyzData = np.array(xyzData)

    driver = ogr.GetDriverByName("ESRI Shapefile")
    dataSource = driver.Open(shp_path, 1)
    layer = dataSource.GetLayer()
    new_field = ogr.FieldDefn("build_H1", ogr.OFTReal)
    layer.CreateField(new_field)
    t = int(layer.GetFeatureCount())
    for i in range(0, t):
        try:
            feat = layer.GetFeature(i)
            geom = feat.GetGeometryRef()
            minX = geom.GetEnvelope()[0]
            minY = geom.GetEnvelope()[2]
            maxX = geom.GetEnvelope()[1]
            maxY = geom.GetEnvelope()[3]
            geometry_polygon=feat.geometry()
            point_zone=[(minX+maxX)/2,(minY+maxY)/2]#X=L,Y=B
            R=((maxX-minX)/2)*((maxX-minX)/2)+((maxY-minY)/2)*((maxY-minY)/2)
            flat_points=radius_cut(point_zone,R,xyzData)
            in_points=isinpolygon(flat_points,geometry_polygon)
            print("均值",np.mean(in_points))
            feat.SetField("build_H1", mean(in_points))
            layer.SetFeature(feat)
        except:
            pass



if __name__ == '__main__':
    shp_path=r"E:\area.shp"
    point_txt=r"E:\point.txt"
    part_statistic_message(shp_path,point_txt)