nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Calculating distance between a table of specific points xy and specific polylines

NinaKruse1

Hello,

I have a question which I can only solve with a code I guess, so I need your help.

Background is, that in a context of population survey, people where asked to tell their zipcode, later one they set a point on googlemaps which result their adress coordinates and so after geocoding their point of adress (including street, zip code, ...). Unfortunetly the two zip codes (first variable and that from the geocoding) do not coincide in 25% of the data. Now I want to calculate the distance between the geocode point (xy) and the zip code border of the zip code the particular person named at the beginning of the survey.

As I have a huge amout of data my question is now how I can calculate the distance between the geocode point and the zip code border which the person told.

So ich have the attribute table + shapefile with the xy coordinates and one with the zip codes.

I thank you for any advice!!

Finde mehr Beiträge mit dem Tag

Accepted answers

All comments

ChrisSnyder

You could write a script to do this, but how about just using the Near tool: http://resources.arcgis.com/en/help/main/10.2/index.html#//00080000001q000000

ChrisSnyder

"and the zip code border of the zip code the particular person named at the beginning of the survey"

Okay, I guess you would have to write a script.

Probably a better/faster way to do this (how many people pnts do you have)? I wonder if the geom objectects still get corrupted if they get put into a Python dictionary?? Perfomace trick would be to access the geometry of the zip code polys as quickly as possible...

Anyway, some code (UNTESTED!) might look like this:

#assumes the pnts and zipcode polys are in the same PRJ, distances are in PRJ's map units...
peoplePntsFC = r"C:\temp\people.shp" #a pnt FC, with a field that has the claimed zip code
zipCodeFC = r"C:\temp\zipcodes.shp" #polygon FC of zipcodes
arcpy.AddField(peoplePntsFC, "DIST_TO_ZIP", "DOUBLE")
updateRows = arcpy.da.UpdateCursor(peoplePntsFC, ["SHAPE@", "ZIP_CODE" "DIST_TO_ZIP"])
for updateRow in updateRows:
   pntObj = updateRow[0]
   zipCode = updateRow[1]
   zipPolyObj = arcpy.da.SearchCursor(zipCodeFC, ["SHAPE@"], "ZIP_CODE = " + str(zipCode)).next()[0]
   updateRow[2] = pntObj.distanceTo(zipPolyObj)
   updateRows.updateRow(updateRow)
del updateRow, updateRows

NinaKruse1

Hey, thanks for your answer, but the point is that I am not looking for the shortest distance to the nearest zip code, but for the distance between a SPECIFIC xy point and a SPECIFIC zipcode border one person named...

You could write a script to do this, but how about just using the Near tool: http://resources.arcgis.com/en/help/main/10.2/index.html#//00080000001q000000

NinaKruse1

Hello, thanks very much for your time to answer and the code!

I have to work with 3100 datasets.

As I have no experience with pyton yet 😞 , I don t know if it can work with the pyton dict.
Could you give me more explanations in the code, as you did at the beginning lines of the code? With that I can learn what the commandlines mean exactly and I know what I do

. Than I will test you code with pleasure.

Thanks!!, Nina

Okay, I guess you would have to write a script.

Probably a better/faster way to do this (how many people pnts do you have)? I wonder if the geom objectects still get corrupted if they get put into a Python dictionary?? Perfomace trick would be to access the geometry of the zip code polys as quickly as possible...

Anyway, some code (UNTESTED!) might look like this:
#assumes the pnts and zipcode polys are in the same PRJ, distances are in PRJ's map units...
peoplePntsFC = r"C:\temp\people.shp" #a pnt FC, with a field that has the claimed zip code
zipCodeFC = r"C:\temp\zipcodes.shp" #polygon FC of zipcodes
arcpy.AddField(peoplePntsFC, "DIST_TO_ZIP", "DOUBLE")
updateRows = arcpy.da.UpdateCursor(peoplePntsFC, ["SHAPE@", "ZIP_CODE" "DIST_TO_ZIP"])
for updateRow in updateRows:
   pntObj = updateRow[0]
   zipCode = updateRow[1]
   zipPolyObj = arcpy.da.SearchCursor(zipCodeFC, ["SHAPE@"], "ZIP_CODE = " + str(zipCode)).next()[0]
   updateRow[2] = pntObj.distanceTo(zipPolyObj)
   updateRows.updateRow(updateRow)
del updateRow, updateRows      

NeilAyres

Loading feature data (including the geometry object) into python dictionaries (or lists, whatever) does indeed work and is a great way to do highly flexible geometry manipulation between features classes.

My logic would be something like this...
Using da.SearchCursor get the points into a dictionary. Presumably each has some sort of point ID so this can become the key thus :
{PntID : [PntZipCode, PntGeomObject]}
The zip code polys, the Zip Code itself can be the key pointing to the PolyGeomObject :
{ZipCode : PolyGeomObject}

Then iterate through the PntID's, extract the PntZipCode, use this to find the record in the ZipCode Polys dict, then use the distanceTo geometry method to find the distance between them.

Collect the distances into a list of PntIDs, Distance... Then use this to update the original Fc or a copy.

Cheers,
Neil

ChrisSnyder

Hi Neil,

I will try it again. In previous efforts (see: http://forums.arcgis.com/threads/9555-Modifying-Permanent-Sort-script-by-Chris-Snyder?p=30332&viewfull=1#post30332) I found that the geom objects somehow got corrupted and basically became useless... but interestingly... only if the map units were in non-metric units. So for example, UTM meters worked fine... State Plane Feet (like I use) was totally messed up for some reason. I think Kim Olivier did some deeper digging on this as I recall. I'll try again and let you know.

ChrisSnyder

As I have no experience with pyton yet

Hi Nina - The sooner you learn, the better.

If it helps, here's the documentation for:
search cursor: http://resources.arcgis.com/en/help/main/10.1/index.html#//018w00000011000000
update cursor: http://resources.arcgis.com/en/help/main/10.1/index.html#/UpdateCursor/018w00000014000000/
geometry objects: http://resources.arcgis.com/en/help/main/10.1/index.html#/Geometry/018z00000070000000/

... which are the "bread and butter" of what you are trying to do.

NeilAyres

Chris,
wow, that was a long thread...
Ditto on the "hack-programmer" bit. I go back a long way from aml to avenue and now to python, but learning all the time.

I must admit, all my geom manipulation stuff has been based on either GCS (decimal lat/long) or metric UTM, so I cannot comment on the feet story. But it does sound odd if feet or whatever would make a difference to the geom object.

One of the things I have learned is that if you create a geom object from scratch, you must supply the spatial reference with it, otherwise (esp with GCS based coords) you get weird rounding effects.

Other than that, I am really enthusiastic about the geometry object methods. And its very fast to do this sort of custom spatial analysis using dictionaries.

Currently on v10.2

Cheers,
Neil

ChrisSnyder

Okay, so I guess the geom objects seem to be working as they should in the dictionary. I can't find that post from Kim Oliver about her research into the topic, but in this case it seems to be working fine. Anyone from ESRI remember something about this?

Anyway, here's the fastest way I could think to do it (took 1.6 seconds to process 1000 pnts). I didn't have any zip code data handy, so I hacked together an example that uses counties instead. My randomly placed points were randomly assigned a county code.

Nina, I added some comments in the code here as you suggested, but without some Python experience they probably won't help much.

#assumes the pnts and polys are in the same PRJ, distances are in PRJ's map units...
import arcpy, time
#a pnt FC, with a field that has the claimed county code
pntsFC = r"C:\csny490\temp\test.gdb\random_pnts"
#polygon FC of counties (with a field for COUNTY_CD)
polyFC = r"C:\csny490\temp\test.gdb\counties"
#add a field to store the distance
arcpy.AddField_management(pntsFC, "DIST_TO_CNTY", "DOUBLE")
#put the county geometry into a dictionary for faster lookups
polyGeomDict = {r[0]:(r[1]) for r in arcpy.da.SearchCursor(polyFC, ["COUNTY_CD","SHAPE@"])}
#start a timer
time1 = time.clock()
#initialize an update cursor
updateRows = arcpy.da.UpdateCursor(pntsFC, ["SHAPE@","CLAIMED_CD","DIST_TO_CNTY"])
#a for loop to iterate through the points
for updateRow in updateRows:
   #get a hook to the 1st two updatecursor field values
   pntObj, countyCode = updateRow[0:2]
   #set the DIST_TO_CNTY value to be the distanace of the pnt to the specified polygon 
   updateRow[2] = pntObj.distanceTo(polyGeomDict[countyCode])
   #commit the update to disk
   updateRows.updateRow(updateRow) 
#del the update cursor objects
del updateRow, updateRows 
#stop the timer
time2 = time.clock()
#report the time it took to calc the distances
print "Took " + str(time2 - time1) + " seconds!"

NeilAyres

Very cool Chris.

polyGeomDict = {r[0]:(r[1]) for r in arcpy.da.SearchCursor(polyFC, ["COUNTY_CD","SHAPE@"])}

Those nested thingys still give me the shakes, but I must try to use them more.
And I presume (can't remember), that if the point is inside the correct poly, the distance returned = 0.
If not the distance is to the nearest boundary.
So if the purpose is to identify those points mis-coded, then a simple text field with "Inside" or "Outside" could also be added.
Lets hope that Nina gets back to you and gives you a tick!!!

Cheers,
Neil

ChrisSnyder

Comprehensions (list and dictionary flavors) are pretty darn cool, but in reality they are just a shortened/efficient form of a for loop. And yes, if the point is inside the poly, it has a distance of 0, otherwise, the closest distance is returned.

Seems like what I recall with geom objects in a dictionary (for me at least a few years back) was that certain reported properties, such as .area or .length were incorrect as compared with the geom object itself (not stored in a dictionary)... basically like the ESRI geom object wasn't 100% (but maybe 80%?) compatible with Python. As I recall Jason Scheirer from ESRI also chimed in with some info on this at the time, but again I can't seem to find the post!

Question for any ESRI people in the know (cough... Jason): Are geometry objects stored in Python dictionaries 100% kosher now? Have they always been or am I getting dementia already?

Nimm an der Diskussion teil

Angemeldete Mitglieder können Beiträge verfassen, Updates folgen und mehr. Neu hier? Registriere ein kostenloses Konto.

Schnellzugriffe

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.