Python 等值线生成（TIN三角网）

等值线功能实现

等值点插值

在一条边上按照高程值线性插值，计算出待插值高程的x,y坐标，并将插值的结果实例成一个点类，方便统一管理，如下所示：

#插值函数
def Interpolation(edge,n,i):
    if(edge.BeginPoint.H>edge.EndPoint.H):
            Hmax=edge.BeginPoint
            Hmin=edge.EndPoint
    elif(edge.BeginPoint.H<edge.EndPoint.H):
        Hmin=edge.BeginPoint
        Hmax=edge.EndPoint
    if (n>Hmin.H and n<Hmax.H):
        dy=Hmax.Y-Hmin.Y
        dx=Hmax.X-Hmin.X
        linescale=(n-Hmin.H)/(Hmax.H-Hmin.H)
        x=Hmin.X+linescale*dx
        y=Hmin.Y+linescale*dy
        name="{0}-{1}".format(n,i)
        eqpoint=Point(name,x,y,n)
        return eqpoint
    else:
        return 0

TIN 边界边查找

遍历TIN的每个三角形边，我们认为具有邻接三角形的边在TIN内部，不具有邻接三角形的边是边界边。我们的数据特点是只具有开曲线，所以找到边界边作为某条等值线的生成起点至关重要，代码如下：

#边界边列表
def BorderTri(Triangle_List):
    Border_list=[]
    for i in range(0,len(Triangle_List)):
        if (len(Triangle_List[i].BaseLine.Belonging_Triangle)==1):
            Border_list.append(Triangle_List[i].BaseLine)
        if (len(Triangle_List[i].newLine1.Belonging_Triangle)==1):
            Border_list.append(Triangle_List[i].newLine1)
        if (len(Triangle_List[i].newLine2.Belonging_Triangle)==1):
            Border_list.append(Triangle_List[i].newLine2)
    return  Border_list

单个三角形内生长等值线

在定位的三角形内生长等值线，根据一个三角形有且仅有两个等值点，在三角形的另外两条边上插值第二个等值点，并将所在的边作为新的生长起点，如下所示：

#一个三角形内生长等值线
def Grow_Eq(triangle,eqedge,n):
    point=0
    edges=[triangle.BaseLine,triangle.newLine1,triangle.newLine2]
    for edge in edges:
        if(edge==eqedge):
            continue
        else:
            point=Interpolation(edge,n,edge.Belonging_Triangle[0])
            if(point!=0):
                break
    return edge,point

一条等值线生成

从边界边的一侧按单个三角形生长等值线，直至生长到TIN对侧的另一条边界边一条等值线生长完成。

#查找等值点 输出等值线序列点
def Equivalent_point(Triangle_List,n):
    point=0
    eqlinelist=[]
    Border_list=bor=BorderTri(Triangle_List)
    for borderedge in Border_list:
        edge=borderedge
        point=Interpolation(edge,n,edge.Belonging_Triangle[0])
        if(point!=0):
            eqlinelist.append(point)
            break
    #根据邻接三角形生成等值线
    index=edge.Belonging_Triangle[0]
    [edge,point]=Grow_Eq(Triangle_List[index],edge,n)
    eqlinelist.append(point)
    while(len(edge.Belonging_Triangle)==2):
        index=edge.Belonging_Triangle[1]
        [edge,point]=Grow_Eq(Triangle_List[index],edge,n)
        eqlinelist.append(point)
    return eqlinelist

按等高距生成全部等值线

在TIN的最低点到最高点的高程范围内，按照一定的等高距，由低到高逐步生成每条等值线，直至全局等值线生成，代码如下：

# 全局等高线
def Contour_Line(Point_List,Triangle_List,d):
    eqlinlists=[]
    Hmax=EPSILON
    Hmin=1/EPSILON
    for p in Point_List:
        if(p.H>Hmax):
            Hmax=p.H
        if(p.H<Hmin):
            Hmin=p.H
    n=Hmin+d
    while(n<Hmax):
        eqlinlist=Equivalent_point(Triangle_List,n)
        eqlinlists.append(eqlinlist)
        n=n+d
    return eqlinlists

可视化结果

可视化代码

'''*******************窗体********************'''
# 生成一个窗体界面 命名为window
window=tk.Tk()
window.title('TIN')
window.geometry('1000x750')

# 在窗体上生成一块画布 用于绘制导线图
canvas = tk.Canvas(window,width=800,height=600,bg = 'white')
canvas.place(x=190,y=10)

#画等值线
def Draw_Equivalent_line(eqpointlists,Point_List):
    xyminmax=XYMinMax(Point_List)
    for eqpointlist in eqpointlists:
        for i in range(0,len(eqpointlist)-1):
            gxgy1=GaussToScreenCor(xyminmax,eqpointlist[i].X,eqpointlist[i].Y)
            gxgy2=GaussToScreenCor(xyminmax,eqpointlist[i+1].X,eqpointlist[i+1].Y)
            oval = canvas.create_oval(gxgy1[0]-1, gxgy1[1]-1, gxgy1[0]+1, gxgy1[1]+1)
            ova2 = canvas.create_oval(gxgy1[0]-1, gxgy1[1]-1, gxgy1[0]+1, gxgy1[1]+1)
            Connect_Point(gxgy1,gxgy2,'red')

#********调试********
path="D:\等高线点数据.txt"
point_list=ReadDataTXT(path)
Net=CreatTIN(point_list)
Draw_TIN(point_list,Net[0])
eqs=Contour_Line(point_list,Net[1],30)
Draw_Equivalent_line(eqs,point_list)
window.mainloop()
'''*******************所有代码在此之上********************'''

等值线（未平滑）

全部代码

#插值函数
def Interpolation(edge,n,i):
    if(edge.BeginPoint.H>edge.EndPoint.H):
            Hmax=edge.BeginPoint
            Hmin=edge.EndPoint
    elif(edge.BeginPoint.H<edge.EndPoint.H):
        Hmin=edge.BeginPoint
        Hmax=edge.EndPoint
    if (n>Hmin.H and n<Hmax.H):
        dy=Hmax.Y-Hmin.Y
        dx=Hmax.X-Hmin.X
        linescale=(n-Hmin.H)/(Hmax.H-Hmin.H)
        x=Hmin.X+linescale*dx
        y=Hmin.Y+linescale*dy
        name="{0}-{1}".format(n,i)
        eqpoint=Point(name,x,y,n)
        return eqpoint
    else:
        return 0
#边界边列表
def BorderTri(Triangle_List):
    Border_list=[]
    for i in range(0,len(Triangle_List)):
        if (len(Triangle_List[i].BaseLine.Belonging_Triangle)==1):
            Border_list.append(Triangle_List[i].BaseLine)
        if (len(Triangle_List[i].newLine1.Belonging_Triangle)==1):
            Border_list.append(Triangle_List[i].newLine1)
        if (len(Triangle_List[i].newLine2.Belonging_Triangle)==1):
            Border_list.append(Triangle_List[i].newLine2)
    return  Border_list
#一个三角形内生长等值线
def Grow_Eq(triangle,eqedge,n):
    point=0
    edges=[triangle.BaseLine,triangle.newLine1,triangle.newLine2]
    for edge in edges:
        if(edge==eqedge):
            continue
        else:
            point=Interpolation(edge,n,edge.Belonging_Triangle[0])
            if(point!=0):
                break
    return edge,point
#查找等值点 输出等值线序列点
def Equivalent_point(Triangle_List,n):
    point=0
    eqlinelist=[]
    Border_list=bor=BorderTri(Triangle_List)
    for borderedge in Border_list:
        edge=borderedge
        point=Interpolation(edge,n,edge.Belonging_Triangle[0])
        if(point!=0):
            eqlinelist.append(point)
            break
    #根据邻接三角形生成等值线
    index=edge.Belonging_Triangle[0]
    [edge,point]=Grow_Eq(Triangle_List[index],edge,n)
    eqlinelist.append(point)
    while(len(edge.Belonging_Triangle)==2):
        index=edge.Belonging_Triangle[1]
        [edge,point]=Grow_Eq(Triangle_List[index],edge,n)
        eqlinelist.append(point)
    return eqlinelist
# 全局等高线
def Contour_Line(Point_List,Triangle_List,d):
    eqlinlists=[]
    Hmax=EPSILON
    Hmin=1/EPSILON
    for p in Point_List:
        if(p.H>Hmax):
            Hmax=p.H
        if(p.H<Hmin):
            Hmin=p.H
    n=Hmin+d
    while(n<Hmax):
        eqlinlist=Equivalent_point(Triangle_List,n)
        eqlinlists.append(eqlinlist)
        n=n+d
    return eqlinlists
'''*******************窗体********************'''
# 生成一个窗体界面 命名为window
window=tk.Tk()
window.title('TIN')
window.geometry('1000x750')

# 在窗体上生成一块画布 用于绘制导线图
canvas = tk.Canvas(window,width=800,height=600,bg = 'white')
canvas.place(x=190,y=10)

# 画出数据点
def Draw_Point(Point_List):
    xyminmax=XYMinMax(Point_List)
    for point in Point_List:
        gxgy=GaussToScreenCor(xyminmax,point.X,point.Y)
        oval = canvas.create_oval(gxgy[0]-1, gxgy[1]-1, gxgy[0]+1, gxgy[1]+1)
        #canvas.create_text(gxgy[0]-13,gxgy[1],text=point.PointName)

# 画出TIN网
def Connect_Point(list1,list2,color):
    line = canvas.create_line(list1[0], list1[1],list2[0],list2[1],fill = color)
def Draw_TIN(Point_List,Line_List):
    xyminmax=XYMinMax(Point_List)
    for line in Line_List:
        gxgy1=GaussToScreenCor(xyminmax,line.BeginPoint.X,line.BeginPoint.Y)
        gxgy2=GaussToScreenCor(xyminmax,line.EndPoint.X,line.EndPoint.Y)
        oval = canvas.create_oval(gxgy1[0]-1, gxgy1[1]-1, gxgy1[0]+1, gxgy1[1]+1)
        ova2 = canvas.create_oval(gxgy1[0]-1, gxgy1[1]-1, gxgy1[0]+1, gxgy1[1]+1)
        Connect_Point(gxgy1,gxgy2,'black')
#画等值线
def Draw_Equivalent_line(eqpointlists,Point_List):
    xyminmax=XYMinMax(Point_List)
    for eqpointlist in eqpointlists:
        for i in range(0,len(eqpointlist)-1):
            gxgy1=GaussToScreenCor(xyminmax,eqpointlist[i].X,eqpointlist[i].Y)
            gxgy2=GaussToScreenCor(xyminmax,eqpointlist[i+1].X,eqpointlist[i+1].Y)
            oval = canvas.create_oval(gxgy1[0]-1, gxgy1[1]-1, gxgy1[0]+1, gxgy1[1]+1)
            ova2 = canvas.create_oval(gxgy1[0]-1, gxgy1[1]-1, gxgy1[0]+1, gxgy1[1]+1)
            Connect_Point(gxgy1,gxgy2,'red')
#********调试********
path="D:\等高线点数据.txt"
point_list=ReadDataTXT(path)
Draw_Point(point_list)
Net=CreatTIN(point_list)
Draw_TIN(point_list,Net[0])
eqs=Contour_Line(point_list,Net[1],10)
Draw_Equivalent_line(eqs,point_list)
window.mainloop()
'''*******************所有代码在此之上********************'''