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.

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

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

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

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.