topic How can I use geometry method snapToLine to find new coordinates? in Python Questions

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

RebeccaStrauch__GISP — Sun, 12 Dec 2021 00:31:53 GMT

Using 10.3.1.

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.

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

DanPatterson_Retired — Sun, 12 Dec 2021 00:31:56 GMT

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.

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

RebeccaStrauch__GISP — Sun, 12 Dec 2021 00:31:59 GMT

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.

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

DanPatterson_Retired — Tue, 09 Feb 2016 03:26:32 GMT

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

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

DanPatterson_Retired — Tue, 09 Feb 2016 03:28:45 GMT

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

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

RebeccaStrauch__GISP — Tue, 09 Feb 2016 03:35:38 GMT

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.

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

RebeccaStrauch__GISP — Tue, 09 Feb 2016 03:49:15 GMT

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.

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

NeilAyres — Tue, 09 Feb 2016 07:39:32 GMT

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.

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

DanPatterson_Retired — Tue, 09 Feb 2016 11:08:37 GMT

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.

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

RebeccaStrauch__GISP — Tue, 09 Feb 2016 15:56:32 GMT

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.

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

RebeccaStrauch__GISP — Tue, 09 Feb 2016 16:04:22 GMT

Thanks for the script, I'll take a look at it a bit more. I did have a couple SearchCursor in one of my versions, but was rethinking it trying to use the numpy arrays.

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

DanPatterson_Retired — Tue, 09 Feb 2016 16:47:24 GMT

I have the numpy version.... just can't find it yet, it is in my notes. I worked on the searchcursor logic first, then speed things up using numpy. I hope this is not mission critical but an incremental for interest need, since I am on a variety of other tangents.

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

RebeccaStrauch__GISP — Tue, 09 Feb 2016 16:54:31 GMT

Only mission critical for me, and I have enough to go on now. Fresh day.....fresher mind. Many ways to get any task done. If you eventually write something up on your blog...that would be a good read.

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

WesMiller — Tue, 09 Feb 2016 20:09:23 GMT

Have you seen this thread arcpy - Standalone Python script to split a polyline with a point layer - Geographic Information Systems Stack Exchange

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

RebeccaStrauch__GISP — Tue, 09 Feb 2016 20:18:17 GMT

I think I did see that one yesterday in my many searches. What it will come down to is a combo of a few....and I'm starting to swing back the direction of using Near_analysis and then maybe Join Field. But, keep having other things pop up, so a start and stop testing session. The line I snap to will be created on the fly, (contour) which is why most of the semi- OTB (out of the box) solutions won't work.

I'll post what I end up with. But thanks for the link. It doesn't hurt for me to look at it again.

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

RebeccaStrauch__GISP — Wed, 10 Feb 2016 16:19:59 GMT

I'm marking Dan's response as the closest to the answer, since he had a script that helped me figure it out. I also found some of my older python scripts that I was doing something similar.

I find the description and examples in the Help lacking. For me, sample code goes a long way.

Thanks for all the suggestions.

edit: btw, no code cleaned up enough to make if worth posting yet.