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

Question

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.

enter image description here

Answer 1

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()

enter image description here

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

• Rex: 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.

Answer 2

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.

2. Select your features (as you already did)

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

Answer 3

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

• FelixIP: +1 nice one, although I'd use fit (Y,X) because it goes S=>N
• FelixIP: If I knew that polyfit exists, the code would be half the size.

Answer 4

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.