How to add spatial reference information in distance calculation?

[RajeshPatnaik]

I have a long polyline with multiple vertices. i want shape length of only a portion of this polyline. For that I have extracted the coordinates of all vertices of this long polyline and created a list. I have removed the coordinates from this list which are not required, so that i have a list of coordinates only of the portion whose shape length is required.

Now i need to calculate the shape length using this coordinates list. For that i have used the following code:

import numpy as np
def dist_from_point_list(point_list):
    line = np.array(point_list, float)
    return np.sqrt(np.sum((line[1:] - line[:-1])**2, -1)).sum()
print(dist_from_point_list(point_list))

When I measure the distance in arcgis pro, the shape length is different from what i get using this function in my code. I have to somehow use the spatial reference/ projection information in my code to get the correct spatial length.

Any help in this regard would be highly appreciated. Thanks.

[Anonymous User]

You could create a temp feature and then measure it by using the shape@length and passing the spatial reference into the cursor.

feature_info = [[[1, 2], [2, 4], [3, 7]],
                [[6, 8], [5, 7], [7, 2], [9, 5]]]

# A list that will hold each of the Polyline objects
features = []

for feature in feature_info:
    # Create a Polyline object based on the array of points
    # Append to the list of Polyline objects
    features.append(
        arcpy.Polyline(
            arcpy.Array([arcpy.Point(*coords) for coords in feature])))

# Persist a copy of the Polyline objects using CopyFeatures
tmp = arcpy.CopyFeatures_management(features, "memory/lineSegment")
length = [f[0] for f in arcpy.da.SearchCursor(tmp, 'SHAPE@LENGTH', spatial_reference=arcpy.SpatialReference(4326))][0]

print(length)

[JohannesLindner]

Convert back to arcpy.Polyline and call getLength():

point_list = [[1000, 1000], [2000, 2000], [3000, 3000], [4000, 4000]]

point_array = arcpy.Array([arcpy.Point(p[0], p[1]) for p in point_list])
polyline = arcpy.Polyline(point_array, spatial_reference=arcpy.SpatialReference(25832))

polyline.getLength("PLANAR", "METERS")  # uses cartesian coordinates
# 4242.64
polyline.getLength("GEODESIC", "METERS")  # considers Earth's curvature
# 4231.36

[Anonymous User]

I like this one better.

[RajeshPatnaik]

 Thanks everyone!
Yes, I tried both the methods. But it gives wrong length. I have attached herewith shape file.
First i used following code to find the  list coodinates.

cable_dict = {}
for row in arcpy.da.SearchCursor(infc, ["OID@", "SHAPE@", "name"]):
    if row[1] != None:
        print("Feature {}:".format(row[0]))
        partnum = 0
        pt_list = []
        # Step through each part of the feature
        for part in row[1]:
            for pnt in part:
                if pnt:
                    pt_list.append((pnt.X, pnt.Y))
                else:
                    print("Interior Ring:")
            partnum += 1
        print(pt_list)
        cable_dict[row[2]] = pt_list

Then i found the coordinate of the feature which i used in the following code.

feature_info = [[(469741.90079999994, 233745.38800000027), (469741.4121000003, 233745.19480000064), (469741.1469999999, 233745.08999999985), (469742.0535000004, 233743.93400000036), (469742.95999999996, 233742.77800000086), (469743.1457000002, 233741.93960000016), (469744.3021, 233736.71729999967), (469744.6370000001, 233735.20500000007), (469744.9084000001, 233734.10789999925), (469744.95710000023, 233733.9110000003), (469745.82519999985, 233730.40210000053), (469745.9199999999, 233730.0189999994), (469747.1799999997, 233724.37299999967), (469746.95799999963, 233723.80000000075), (469747.4139999999, 233722.3440000005), (469747.56900000013, 233721.86690000072), (469747.6469999999, 233721.73399999924), (469750.0284000002, 233713.38040000014), (469750.4979999997, 233711.73299999908), (469751.28199999966, 233710.97399999946), (469751.4919999996, 233710.39699999988), (469825.6730000004, 233504.0979999993), (469825.26400000043, 233504.92899999954), (469813.48900000006, 233531.3560000006), (469813.3729999997, 233531.6129999999), (469801.8559999997, 233561.35799999908), (469800.2570000002, 233563.26400000043), (469797.0800000001, 233570.77600000054), (469793.90299999993, 233578.28800000064), (469792.6009999998, 233580.41200000048), (469797.7457999997, 233582.2851999998), (469799.33100000024, 233582.86240000091), (469801.6150000002, 233583.69400000013), (469794.9391999999, 233599.6598000005), (469792.40390000027, 233605.72309999913), (469789.87600000016, 233611.7689999994), (469789.3559999997, 233613.37900000066), (469786.5020000003, 233620.32799999975), (469783.6370000001, 233627.3230000008), (469780.784, 233634.3440000005), (469777.8969999999, 233641.32200000063), (469775.2419999996, 233648.27800000086), (469772.51400000043, 233655.12189999968), (469769.9369999999, 233661.77999999933), (469765.5580000002, 233672.3900000006), (469764.8053000001, 233674.33830000088), (469761.61199999973, 233682.60400000028), (469757.9497999996, 233691.7310000006), (469757.62590000033, 233692.5382000003), (469756.5080000004, 233695.3239999991), (469754.17619999964, 233702.33110000007), (469753.43620000035, 233704.55440000072), (469752.5036000004, 233707.35669999942), (469751.57100000046, 233710.15909999982), (469751.4919999996, 233710.39699999988), (469750.87399999984, 233710.21829999983), (469914.5789999999, 233428.50999999978), (469914.88609999977, 233428.49760000035), (469916.05510000046, 233428.45069999993), (469918.4615000002, 233426.21199999936), (469921.8501000004, 233423.05959999934), (469930.1710000001, 233473.81299999915), (469929.324, 233474.09200000018), (469928.5959999999, 233470.64100000076), (469927.5410000002, 233466.39790000021), (469923.9411000004, 233456.26070000045), (469922.7379999999, 233452.87299999967), (469920.59530000016, 233446.5427000001), (469920.07799999975, 233445.01449999958), (469918.75739999954, 233441.11319999956), (469918.43099999987, 233440.1490000002), (469916.59669999965, 233433.80859999917), (469916.2450000001, 233432.59300000034), (469914.65809999965, 233428.70390000008), (469914.5789999999, 233428.50999999978), (469914.44600000046, 233428.54900000058), (469914.0614, 233428.11099999957), (469908.3674999997, 233421.62629999965), (469905.82799999975, 233418.73399999924), (469906.2189999996, 233418.75), (469906.31900000013, 233418.66200000048), (469912.2860000003, 233413.43989999965), (469912.1730000004, 233413.30599999987), (469911.7759999996, 233413.22199999914), (469904.25320000015, 233419.67679999955), (469901.9397, 233421.66169999912), (469899.4305999996, 233423.81460000016), (469889.17399999965, 233432.61500000022), (469880.9966000002, 233439.99540000036), (469875.89800000004, 233444.59699999914), (469861.96800000034, 233458.4539999999), (469852.60720000044, 233468.59249999933), (469845.44099999964, 233476.35400000028), (469838.574, 233483.48179999925), (469835.5120000001, 233486.66000000015), (469834.8194000004, 233488.11759999953), (469832.3590000002, 233493.29499999993), (469830.5712000001, 233496.18359999917), (469825.6730000004, 233504.0979999993), (469824.4938000003, 233500.06020000018)]]
features = []

for feature in feature_info:
    # Create a Polyline object based on the array of points
    # Append to the list of Polyline objects
    features.append(arcpy.Polyline(arcpy.Array([arcpy.Point(*coords) for coords in feature])))

# Persist a copy of the Polyline objects using CopyFeatures
tmp = arcpy.CopyFeatures_management(features, "memory/lineSegment")
length = [f[0] for f in arcpy.da.SearchCursor(tmp, 'SHAPE@LENGTH', spatial_reference=arcpy.SpatialReference(27700))][0]

print(length)

The length found is : 1074.598630135699
whereas the shape length from arcgis pro is : 471.116806
I have attched herewith the shape file which i used my testing.

Why this difference in length??

[DanPatterson]

# -- z is your array of points
diff = z[1:] - z[:-1]
np.sum(np.sqrt(np.einsum('ij,ij->i', diff, diff)))
1074.598814739643

You have duplicate points in the geometry, so I suspect that there shouldn't be and you are connecting segments that shouldn't be.  You need to check your geometry (plot it or something)

uniq, idx, cnts = np.unique(z, return_index=True, return_counts=True, axis=0)
uniq.shape
(97, 2)

np.where(cnts>1)[0]
 array([21, 55, 79], dtype=int64)

# duplicates for 21, 55, 79

 I think separate parts are being connected that shouldn't be 

[JohannesLindner]

Like Dan said, you are connecting things in your code that aren't connected in your feature.

This is the faulty part of your code:

        pt_list = []
        # Step through each part of the feature
        for part in row[1]:
            for pnt in part:
                if pnt:
                    pt_list.append((pnt.X, pnt.Y))

You go though all parts of a multiplart polyline and append all points into a single array. But that introduces connections ( => line segments => additional length) that aren't there in your feature.

Take the northern part of your polyline (I exploded it into single parts):

This is part 5:

This is part 6:

This is part 7:

The red vertices are the ends of the line parts.

In your code, you loop through the line parts and append the points from start to end of the line part. That means that you suddenly have a virtual line from the end of part 6 to the start of part 7:

This accounts for 36.5 meters of your error. Add in the other virtual lines you introduced and you got your difference.

[RajeshPatnaik]

Thanks @DanPatterson @JohannesLindner  for pointing out the reason of difference in length.
The reason of the difference in length is clear. 🙂
I have to calculate the length of each part separately and finally I have to sum them up. This will give correct length.
Thanks all.