# Creating best fit line from point data using ArcGIS for Desktop?

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. 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")
plines = arcpy.mapping.ListLayers(mxd,"lines")

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.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)
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()); 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.
– Rex
Oct 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.

1. Use "Add XY Coordinates (Data Management)" to add X and Y coordinates fields.

3. Use View Menu > Graphs > Create Graphs.

4. 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

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)
y.append(item)
regression = numpy.polyfit(x, y, 1)
m = float(regression)
b = float(regression)
x1 = (YMin-b)/m
x2 = (YMax-b)/m
feature_info = [[x1, YMin], [x2,YMax]]
lineArray = arcpy.Array()
for x,y in feature_info:
features = arcpy.Polyline(lineArray)
arcpy.CopyFeatures_management(features, regressionline)
arcpy.Clip_analysis(regressionline, outFCmbg, outFC)
print "Finished Without Errors!"
except:
tb = sys.exc_info()
tbinfo = traceback.format_tb(tb)
pymsg = "PYTHON ERRORS:\nTraceback info:\n" + tbinfo + "\nError Info:\n" + str(sys.exc_info())
msgs = "ArcPy ERRORS:\n" + arcpy.GetMessages(2) + "\n"