(三十五)arcpy开发&计算polyline的起点和结束点

今天我们来学习一下,使用arcpy来寻找polyline要素类的起始端点和结束端点。首先,程序开始会对要素类polyline进行遍历,然后利用半正矢公式对两个端点进行计算。这里主要利用到半正矢公式,以及在arcpy中创建要素类字段,并添加相应的值。我们来看一下具体的实现代码。

import arcpy
import sys
import math

def haversine(point1, point2):
    lat1, long1 = point1
    lat2, long2 = point2
    diffLat = math.radians(lat2 - lat1)
    diffLong = math.radians(long2 - long1)
    radius = 6378137
    a = (math.sin(diffLat / 2) * math.sin(diffLat / 2) +
         math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *
         math.sin(diffLong / 2) * math.sin(diffLong / 2))
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    d = c * radius
    return d

def find_distant_pts(allpts, beginpts, endpts, second_check=False):
    beg = beginpts[0]
    end = endpts[0]
    biggest_distance = 0
    if len(allpts) > 2:
        for pt in allpts:
            if allpts.count(pt) == 1:
                for pt2 in allpts:
                    if allpts.count(pt2) == 1 or second_check is True:
                        test_distance = haversine(pt, pt2)
                        if test_distance > biggest_distance:
                            biggest_distance = test_distance
                            if pt in beginpts:
                                beg = pt
                                end = pt2
                            else:
                                beg = pt2
                                end = pt
    return beg, end

def calc_endpoint_coord(in_polyline):
    beginlat = "BegLat"
    beginlong = "BegLong"
    endlat = "EndLat"
    endlong = "EndLong"


    field_list = [beginlat, beginlong, endlat, endlong]
    existing = [field.name for field in arcpy.ListFields(in_polyline)]
    field_error = False
    for f in field_list:
        if f not in existing:
            arcpy.AddField_management(in_polyline, f, "DOUBLE")
        else:
            arcpy.AddError("ERROR - field named {} already exists".format(f))
            field_error = True
    if field_error:
        sys.exit()


    with arcpy.da.UpdateCursor(in_polyline, ['SHAPE@', beginlong, beginlat, endlong, endlat]) as cur:
        for row in cur:

            beginpts = [(part.getObject(0).Y, part.getObject(0).X) for part in row[0]]
            endpts = [(part.getObject(part.count - 1).Y, part.getObject(part.count - 1).X) for part in row[0]]

            allpts = beginpts + endpts

            beg, end = find_distant_pts(allpts, beginpts, endpts)

            if beg == end:
                beg, end = find_distant_pts(allpts, beginpts, endpts, True)

            row[1] = beg[1]
            row[2] = beg[0]
            row[3] = end[1]
            row[4] = end[0]
            cur.updateRow(row)

in_polyline=r"C:\Users\qin\Desktop\zxy\shp_lb.shp"
calc_endpoint_coord(in_polyline)

在代码的中,注意使用for循环的便捷形式。

beginpts = [(part.getObject(0).Y, part.getObject(0).X) for part in row[0]]
            endpts = [(part.getObject(part.count - 1).Y, part.getObject(part.count - 1).X) for part in row[0]]

最后,我们来看一下具体的实现结果。



                               更多内容,请微信扫二维码关注公众号,或者加入arcpy开发qq学习群:487352121

                                                                     

转载自:https://blog.csdn.net/u010608964/article/details/89412976

You may also like...