0

I am trying to join points to polygons based on whether their addresses match as well as if the distance between them is less than 2500 feet. I defined a haversine function that computes the distance between the point and the polygon using their coordinates. However I get an error that the sql is invalid.

import arcpy
import harversine
arcpy.env.overwriteOutput=True
arcpy.env.workspace=r"C:\mygdb.gdb"        
where="pointshp.CAdd = polyshp.P_Address and harversine.myhaversine(pointshp.NewLong,pointshp.NewLat,polyshp.polyX,polyshp.polyY) < 2500"
keyfield="pointshp.OBJECTID"
arcpy.MakeQueryTable_management(["pointshp","polyshp"],"joined","USE_KEY_FIELDS",keyfield,"",where)
arcpy.CopyFeatures_management("joined","joined")

The error I get:

An invalid SQL statement was used. [joined]

The haversine function I use which is saved in a separate script called harversine.py i.e. it is not a spelling mistake.

import math
def myhaversine(lon1, lat1, lon2, lat2):
    # convert decimal degrees to radians 
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    c = 2 * math.asin(math.sqrt(a)) 
    r = 20887680 # Radius of earth in kilometers. Use 3956 for miles, 20,887,680 feet
    return math.ceil(c * r)
  • You simply cannot use formulae in your query, e.g. try this query ([FID] mod 2) = 0 – FelixIP Nov 3 '15 at 22:10
  • Do you know if there's another way to accomplish what I am trying to do? – ketar Nov 3 '15 at 22:13
  • Spatial join polygons one to many. Populate new field with your distance thing between pairs. Similar to make query table, the only difference is field – FelixIP Nov 3 '15 at 22:27
  • Unfortunately the spatial join won't cater to the condition of point address being equal to the polygon address. pointshp.CAdd = polyshp.P_Address – ketar Nov 3 '15 at 22:49
  • Use select by location, could be very slow though – FelixIP Nov 3 '15 at 22:54
2

I think your best bet will be to determine all of the addresses that fit the distance criteria first, and then filter for those. Here's one way to do it. Determine the point X and Y values from the fields, and add them to a dictionary with the key being the address. Then perform the calculations per address while iterating through the polygons, and if the address meets the distance criteria add the address to a 'good address' list. Finally incorporate those addresses into the where clause. It's all tabular, so it should be relatively quick.

The code (untested):

import arcpy
import harversine
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\mygdb.gdb"

print "creating point dictionary"
#Dictionary with XYs by address from point
pntDi = dict ([(addr, (x, y)) for addr, x, y
                in arcpy.da.SearchCursor
                ("pointshp", ["CAdd", "NewLong", "NewLat"])
                if addr and x and y])

#empty list for good addresses
goodAddresses = []

#iterate polygon table
print "iterating polygons and calculating"
with arcpy.da.SearchCursor ("polyshp", ["P_Address", "polyX", "polyY"]) as cursor:
    for addr, polyX, polyY in cursor:
        if not addr or not polyX or not polyY:
            continue
        if not addr in pntDi:
            continue

        pntX = pntDi[addr][0]
        pntY = pntDi[addr][1]

        #Calculate value
        value = harversine.myhaversine(pntX,
                                       pntY,
                                       polyX,
                                       polyY)

        #Add good value to addresslist
        if value < 2500:
            goodAddresses += [addr]


#New where clause
addressStr = "', '".join (goodAddresses)
where = "pointshp.CAdd = polyshp.P_Address and pointshp.CAdd in ('{0}')".format (addressStr)
  • Let me test the code and see what results I get. – ketar Nov 4 '15 at 15:15
  • fixed the line to if not addr or not polyX or not polyY: – Emil Brundage Nov 4 '15 at 15:17
  • Your code definitely works. My only problem is that the dictionary leaves out some addresses because they are duplicates yet i would like for each point address that fulfills the conditions to be assigned a polygon even if it's a duplicate address. – ketar Nov 4 '15 at 16:03
  • Try swapping - creating the dictionary from the polygons and iterating the points. – Emil Brundage Nov 4 '15 at 17:13
  • I run into the same problem with the duplicate addresses being dropped. I am going to try using a list of tuples instead of a dictionary and see how that works...What I hadn't mentioned is that if a point address matches 2 polygon addresses, then the point shouldn't be joined to any of the polygons. – ketar Nov 4 '15 at 18:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.