How can I use geometry method snapToLine to find new coordinates?

[RebeccaStrauch__GISP]

Using 10.3.1.

From Geometry—Help | ArcGIS for Desktop

I'm probably missing something basic here, but given a polyline, and a point geometry....how do I use snapToLine(in_point) and how do I capture the new point (in a python script)?

I've tested using arcpy.Near_analysis, which will work, but the description sounds more like what I need.  This is a small part of a bigger script.....which I'm now trying using a numpy array.  A general description, given my FC of classified points, I grab the first point...create the contour....snap the point to the contour and split at the point...then go a given distance in each direction along the line, split and create a transect...capture the coordinates at the three points. Previously I had this in Avenue and I'm updating, and adding new features, of course.  There's a bit more to it, but that's a general description for anyone interested.

Fwiw, I'll include my test script (I'll clean it up once I get this figured out...lots of stuff not being used in in this part.)

import os
import numpy as np
from numpy.lib import recfunctions as rfn  
import arcpy
from arcpy import env
from arcpy.sa import *
#from ADFGUtils import *
#from gpdecorators import *

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Set the Geoprocessing environment...
initialSR = arcpy.SpatialReference(3338)
geoSR = arcpy.SpatialReference(4269)

theWorkspace = r'C:\_beartest\Prep.gdb'
arcpy.env.workspace = theWorkspace
arcpy.env.overwriteOutput = True    
inStudy = r'C:\_beartest\Prep.gdb\Unit20A_boundary'
arcpy.env.Extent = inStudy

inDEM = r"C:\_beartest\Unit20Araster.gdb\elevExtent"  # clip later or \elevClip ..
rasterGDB, rasterFile = os.path.split(inDEM)
arcpy.env.snapRaster = inDEM
outContour = arcpy.os.path.join(theWorkspace, "c3")
outContour2 = arcpy.os.path.join(theWorkspace, "c4")

inPts = r'C:\_beartest\Prep.gdb\Pts1_3338'   # arcpy.ListFeatureClasses("Pts*3338*", "Point")
minTrans = '5000 Meters'    
maxTrans = '20000 Meters'
length, unit = maxTrans.split(" ")
halfTrans = int(length) / 2
flatOption = "fishnet"    

nearSearchRadius = "1000 Meters"
location = "LOCATION"    
angle = "ANGLE"
tmpPt = arcpy.os.path.join(theWorkspace, "tmpPt")

tmpRoute = arcpy.os.path.join(theWorkspace, "tmpRoute")    

field_names = "*"
arr = arcpy.da.FeatureClassToNumPyArray(inPts, field_names, "", initialSR)
dtypes = arr.dtype
arr2 = (np.array(arr, dtype=dtypes)).view(np.recarray) 

cnt = 0
while  cnt < 10:
    print("{0} PtID {1} elev: {2}".format(cnt, arr2.PtID[cnt], arr2.elev_ft[cnt]))
    transtype = arr2.TransType[cnt]
    ptID = arr2.PtID[cnt]
    elevM = arr2.elev_m[cnt]
    origX = arr2.Shape[cnt][0]
    origY = arr2.Shape[cnt][0]
    if transtype == "Riparian":
        print("this is a Riparian")
    elif cnt == 9:
        print("this is a {0}".format(transtype))
        ContourList(inDEM, outContour, elevM) 
        ptGeom = arcpy.PointGeometry(arcpy.Point(origX, origY, 0, 0, ptID))
        #nearPt = arcpy.Near_analysis(ptGeom, outContour, nearSearchRadius, "LOCATION", "ANGLE")
        #newPt = outContour.snapToLine(ptGeom)
    cnt += 1

variables of interest:  ptGeom, outContour

Hopefully it is something I'm just overlooking.

[DanPatterson_Retired]

Rebecca... I can't find the numpy stuff right now, but I did do some stuff for lessons on vector projection of a point to a line.  Attached is one of those package things that contains the results.  I have attached the script as a zip file and her are some results for the results... they will make more sense when you open the project.

>>> Points 
 A: 2 2 NaN NaN
 B: 10 6 NaN NaN
 C: 6 7 NaN NaN
 D: 7.2 4.6 NaN NaN
Distances
 A-B:  8.9442719100 
 A-C:  6.4031242374 
 B-C:  4.1231056256 
 C-D:  2.6832815730

Basically A and B form a line.  Point C is the point near the line and D is the intersection point.,

If you are interested and in no rush, I have Near replicated in a couple of variants, one of the obvious is the arcpy python route and the other entailed numpy (can't remember the details yet.)

This should give you some ideas.  The only trick is in limiting the search radius of the from point to  the candidate to features.

[RebeccaStrauch__GISP]

Thanks Dan,

That should be helpful.  I think it doesn't like the variable/objects going in to the command.

>>> newPt = (outContour.snapToLine(ptGeom))
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
AttributeError: 'str' object has no attribute 'snapToLine'
>>>

I'm thinking it doesn't like my outContour feature class.  But I need to digest the vector_projection_arcpy script a little more and see what geometry worked for you.  I think I just need to step away from it for a while. 

[DanPatterson_Retired]

are you using python 2.7? I am using 3.4.x  I will have to look at it in the morning and I get off the iThingy

PS   The demo works with points I created in the script... it is not the version for using with other files!  It just demonstrates the principle, you will have to change the output file names in the script.  This package thing might be an issue, the shapefiles where in a folder called shapefiles in the folder that the mxd was in.  Also, I don't use the arcmap python IDE, I have used the script in pyscripter and pythonwin and the script must be in the same folder as the mxd...if that package thing produces one... it is the first time I used them

[RebeccaStrauch__GISP]

yes, using python 2.7 with Desktop 10.3.1.

And yes, I realize you are hardcoding the points/lines in the demo.  Just need to translate my generated contour line as input.  I do think that is the issue....that is, the polyline geometry, not the point.  looked at it too long this afternoon.....brain is shot.

[NeilAyres]

Rebecca,

in those geometry methods, outContour would have to be a geometry object, not an entire feature class.

I would load those 2 sets of geometries into a dict or something then get the logic right of cycling through them.

[RebeccaStrauch__GISP]

Thanks Neil.  The outContour is actually created on the fly....the ptGeom object I am creating as I cycle thru the points.  The Near_analysis gets me the coord I want in the Near_X, Near_Y fields, and I could use this to split the line and go from there, but then I need to copy all the attributes over.....do my other steps...then the final transect out to a new file.

The snapToLine looked promising since it looks to 1) actually move the point...which would be helpful, but not mandatory, and 2) keep the attributes.

In my sleep (where many of my programming problems are solved) I think I know how I should approach it (getting the outContour to a geom object, as you mentioned).  I'll give that a go.  Thanks for your input.

[DanPatterson_Retired]

PS

You might be interested in this... which is part of the bigger Near picture.  It entails using Numpy to determine euclidean distances between point objects

16 ... NumPy vs SciPy... making a point

I tended to focus more on the speed issue between pure numpy and the implementation in SciPy... which wasn't available unless you did a special install... partial conclusion for 50,000,000 million origin-destination points... and my rant

......

1 e7 0.196187973    0.137131929

5 e8 0.960922956    0.664637804

Ok ... are you missing something like me?  Look at the last line... 50 million (50,000,000) Euclidean origin - destination calculations ... one origin and 50,000,000 destinations... yes with a square root calculation too... and in real world (MTM zone 9 coordinates) covering random points in a 50 km x 50 km range. Hmmm.... Ok!  I get it Numpy using einsum takes 44.5% longer to calculate!! That must be it!  Now I am onboard with all those time testers extolling the virtues of one algorithm over an other or one piece of software as being the latest in time saving efforts.

Einsum isn't the only implementation for distance, and now that pro comes with scipy, determining distances is trivial

[RebeccaStrauch__GISP]

ya...the speed difference in this case isn't that critical....it's getting the data out.....and making sure I have it documented enough so the next person can understand it.  Sometimes sticking to the standard arcpy commands are the best way to make sure it is readable to others....especially if the next person is an =ologist type, not a GIS person.

I'm thinking my issue is one of two things..

ptGeom = arcpy.PointGeometry(arcpy.Point(origX, origY, 0, 0, ptID))

so I tried a temp point with just the arcpy.Point part of it.  Or more likely it doesn't like line the contour.  You example shows

baseLine = arcpy.Polyline(arcpy.Array([A,B]),SR)

So, I'm thinking I need to isolate the contour line segments and grab the geometry of that...instead of having it be a featureclass.  I may give that a try.

[DanPatterson_Retired · Accepted solution]

Attached project and a folder with 2 shapefiles and a script to demo how to get the below.

All it demonstrates is how to get the lines connecting a bunch of random points to a polyline...

This is not a full implementation obviously but it is part of the thought process.