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

2664
11
01-16-2014 11:36 AM
NinaKruse1
New Contributor
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!!
Tags (2)
0 Kudos
11 Replies
ChrisSnyder
Regular Contributor III
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
0 Kudos
ChrisSnyder
Regular Contributor III
"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      
0 Kudos
NinaKruse1
New Contributor
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
0 Kudos
NinaKruse1
New Contributor
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      
0 Kudos
NeilAyres
MVP Alum
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
0 Kudos
ChrisSnyder
Regular Contributor III
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&viewfu...) 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.
0 Kudos
ChrisSnyder
Regular Contributor III
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.
0 Kudos
NeilAyres
MVP Alum
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
0 Kudos
ChrisSnyder
Regular Contributor III
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!" 
0 Kudos