I have GPS points that show the locations of surface bubbles from a research dive along a transect. Wind, waves, and currents have conspired to add error to these data. I want to create a single best fit line from the points highlighted in blue in the attached picture using ArcGIS 10.3 but I do not know how to do that.
4 Answers
Unfortunately solution by Farid Cher uses regression analysis. It minimises either (X-distance)^2 to line, or (Y-distance)^2, depending on what values were picked for Y axis. It seems that you’d like to minimise distance to line from points.
Complete solution can be found here: https://math.stackexchange.com/questions/839464/how-to-find-a-line-that-minimizes-the-average-squared-perpendicular-distance-fro but it’s to much effort.
Approximate solution can be achieved by using average of XY regression and YX regression lines.
Try this script:
import arcpy, traceback, os, sys
import numpy as np
try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
mxd = arcpy.mapping.MapDocument("CURRENT")
points = arcpy.mapping.ListLayers(mxd,"points")[0]
plines = arcpy.mapping.ListLayers(mxd,"lines")[0]
g=arcpy.Geometry()
geometryList=arcpy.CopyFeatures_management(points,g)
geometryList=[p.firstPoint for p in geometryList]
SX,SY,SX2,SXY,SY2=0,0,0,0,0
minX=geometryList[0].X
maX=minX
N=len(geometryList)
for p in geometryList:
SX+=p.X;SX2+=p.X*p.X;SY+=p.Y;SXY+=p.X*p.Y;SY2+=p.Y*p.Y
if p.X<minX:minX=p.X
if p.X>maX:maX=p.X
# y regression
A=np.array([[SX,N],[SX2,SX]])
B=np.array([SY,SXY])
(a,c)=np.linalg.solve(A,B)
# X regression
A=np.array([[SY,N],[SY2,SY]])
B=np.array([SX,SXY])
(A,C)=np.linalg.solve(A,B)
a=(a+1/A)/2
c=(c-C/A)/2
p1=arcpy.Point(minX,a*minX+c)
arr=arcpy.Array(p1)
p2=arcpy.Point(maX,a*maX+c)
arr.add(p2)
pLine=arcpy.Polyline(arr)
curT = arcpy.da.InsertCursor(plines,"SHAPE@")
curT.insertRow((pLine,))
del mxd
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
Note, script will work on selection.
On the example shown average distance to Y regression line was 444 m, distance to 'Min line' was 421 m
-
One note. In your code, I believe you need to closeout your cursor with 'del curT' or use it within a 'with' statement or else it won't commit the row changes.– RexOct 28, 2022 at 13:48
This is not possible with ArcGis built-in tools to draw a regression line fit to your point features geographically.
Instead you should use Graph tool to create a regression line.
Use "Add XY Coordinates (Data Management)" to add X and Y coordinates fields.
Select your features (as you already did)
Use View Menu > Graphs > Create Graphs.
Fill the parameters (select X, Y coordinate fields) and then add a new function of type slope
FelixIP, thanks! That was great! It worked like a charm. I did code a polynomial regression for ArcGIS (not nearly as elegant as yours...my code below). There was a slight difference in the line.
try:
'''This Tool will take a feature class of points and it will create a
best fit line based on a polynomial regression of the x, and y values.
This tool only works on projected datasets.
The output regression line will be clipped to a minimum bounding circle
around the points.
This tool was developed and tested with an ArcGIS 10.3 with Python 2.7
This tool will work with an ESRI Basic level license.'''
import sys, traceback
print "go"
import arcpy, numpy
arcpy.env.overwriteOutput = True
#Set your data paths here....
inFC = r"C:\gTemp\divetract.shp"
outFC = r"C:\gTemp\aaaregressionline.shp"
outFCmbg = r"in_memory\outFCmbg"
regressionline = r"in_memory\regressionlinetemp"
arcpy.MinimumBoundingGeometry_management(inFC, outFCmbg, "CIRCLE", "ALL", "", "NO_MBG_FIELDS")
desc = arcpy.Describe(outFCmbg)
extent = desc.extent
XMin = extent.XMin
YMin = extent.YMin
XMax = extent.XMax
YMax = extent.YMax
array = arcpy.da.FeatureClassToNumPyArray(inFC, ["SHAPE@XY"])
x = []
y = []
for item in array:
x.append(item[0][0])
y.append(item[0][1])
regression = numpy.polyfit(x, y, 1)
m = float(regression[0])
b = float(regression[1])
x1 = (YMin-b)/m
x2 = (YMax-b)/m
feature_info = [[x1, YMin], [x2,YMax]]
lineArray = arcpy.Array()
for x,y in feature_info:
lineArray.add(arcpy.Point(x,y))
lineArray.add(lineArray.getObject(0))
features = arcpy.Polyline(lineArray)
arcpy.CopyFeatures_management(features, regressionline)
arcpy.Clip_analysis(regressionline, outFCmbg, outFC)
print "Finished Without Errors!"
except:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info()[1])
msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"
arcpy.AddError(pymsg)
arcpy.AddError(msgs)
print pymsg + "\n"
print msgs
If you want a quick approximation, you can create your own line by adding a segment feature by placing two vertices on the map that correspond to the start and end of your transects. Next, examine the coordinates of the vertices, either with the Information button (i in a circle) or by looking deeper at the table of vertices. With the start and end coordinates, you can calculate the slope of the line and its intercept.
This works best when you're creating lines for display purposes and can be faster than figuring it out with precision if you have <50 lines, but if you are interested in doing this for hundreds of lines, a more automated approach would be better.