5

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

4

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

| improve this 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

| improve this answer | |
1

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
| improve this answer | |
  • +1 nice one, although I'd use fit (Y,X) because it goes S=>N – FelixIP Sep 30 '15 at 5:01
  • If I knew that polyfit exists, the code would be half the size. – FelixIP Sep 30 '15 at 8:31
0

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.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.