topic Re: How to add spatial reference information in distance calculation? in Python Questions

How to add spatial reference information in distance calculation?

RajeshPatnaik — Thu, 17 Mar 2022 11:55:55 GMT

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.

Re: How to add spatial reference information in distance calculation?

Anonymous User — Thu, 17 Mar 2022 13:06:23 GMT

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)

Re: How to add spatial reference information in distance calculation?

JohannesLindner — Thu, 17 Mar 2022 13:16:56 GMT

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

Re: How to add spatial reference information in distance calculation?

Anonymous User — Thu, 17 Mar 2022 13:31:01 GMT

I like this one better.

Re: How to add spatial reference information in distance calculation?

RajeshPatnaik — Thu, 17 Mar 2022 14:58:14 GMT

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??

Re: How to add spatial reference information in distance calculation?

DanPatterson — Thu, 17 Mar 2022 15:35:24 GMT

# -- 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

Re: How to add spatial reference information in distance calculation?

JohannesLindner — Fri, 18 Mar 2022 12:45:02 GMT

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.

Re: How to add spatial reference information in distance calculation?

RajeshPatnaik — Fri, 18 Mar 2022 18:53:05 GMT

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.