AnsweredAssumed Answered

Arcpy to determine distance

Question asked by smkruzik_pa on Nov 30, 2016

Eons ago, our group would identify points for registration forms using a USGS quad broken into 9 segments.  From the top right of one of the segments, it would basically say 9,645 feet left, 7,123 feet down (or something like that).  I've been tasked with converting these to true lat/lon in decimal degrees.  I can usually figure out a python way of doing something like this, but working with the geometry and unit difference is really proving to be difficult for me.  Could anyone please help or lead me in the right direction?

 

I need to write a tool that the user can enter in the original x,y distance in feet.  The tool would then; from the top right corner of the quad; plot the point based on the difference (basically, just slope from the origin) and then return what the true lat/lon would be in decimal degrees.  So in essence, I'm taking an origin of perhaps 41.111, -75.111.  User then puts in 2000 feet across, 1500 feet down.  Python takes care of the slope and plots the point based on the values.  It then returns the new lat/lon; say 41.002, -75.189.

 

here is what I had just from using the ArcGIS help snippets.  I was using the arcpy window to observe the process, but once it's working, I'd move it to a stand-alone script.  The part in the code after the triple hash is where I'm having trouble.  I try and grab the top right vertex point of the selected quad/sector polygon and get it's x,y to use as the origin.  I sort of have it, but can't seem to figure out how to isolate it to use in the math.  Even then, I'm not sure the python functions that I need to use to then get a new location based on the slope in feet from the origin.

 

I apologize for not having the below formatted. I don't use this often, so I don't know how to get it wrapped into the python code block.  I'm sure glad ESRI makes it easy to use emoticons in here though because that certainly makes things easier!(not sure which one is the sarcasm one)

 

import arcpy
from arcpy import mapping
import os

 

mxd = mapping.MapDocument("CURRENT")
FC = mapping.ListLayers(mxd)[2] #this is the quad cut into 9 sectors and is identified on the registration forms

 

# the following 2 variables would be populated by the user, but for now are dummy values for testing

Quad = "'" + "TAMAQUA" + "'"

Sector = 3

 

# find the record in the feature class based on above

arcpy.SelectLayerByAttribute_management(FC,"NEW_SELECTION","[QUAD_NAME] = " + Quad + "AND [SECT_NUM] =" + str(Sector))

 

### this is where it gets tricky...

for row in arcpy.da.SearchCursor(FC, ["OID@", "SHAPE@"]):
     print "Feature {0}:".format(row[0])
     partnum = 0
     for part in row[1]:
         print "Part {0}:".format(partnum)
         for pnt in part:
             if pnt:
                 print "{0}, {1}".format(pnt.X, pnt.Y)
             else:
                 print "Interior Ring"
         partnum += 1

 

 

 

The results I get are:

Feature 3930:
Part 0:
-75.916347055, 40.875087773
-75.874642422, 40.875088341
-75.874642491, 40.833417869
-75.916347289, 40.833417228
-75.916347055, 40.875087773

Outcomes