Select to view content in your preferred language

Test for XY inside Polygon

4416
25
Jump to solution
09-01-2020 03:05 PM
jaykapalczynski
Honored Contributor

Im trying to do this and stuck

  1. Take a JSON string with address and distance set those to variables
  2. Geocode the Address to get an XY and set those to Variables
  3. Buffer a polygon of the state using the distance variable created from the DataInput JSON 
  4. Define the State Boundary
  5. Buffer the State Boundary and set into memory
  6. Test ot see if the point is within that Buffered State Boundary
  7. Return True if its within and False if Outside the Buffered State Boundary

ERROR:

Traceback (most recent call last):
File "E:\ArcGIS Server\ArcGIS\xxxxxxxx.py", line 87, in <module>
results = arcpy.SelectLayerByLocation_management(ptGeometry, 'INTERSECT', 'IN_MEMORY/BufferGeom')
File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 7891, in SelectLayerByLocation
raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000368: Invalid input data.
Failed to execute (SelectLayerByLocation).

Can Anyone see anything obvious as to why I am tripping this error?

import arcpy
import requests
import json
import uuid

#dataInputParameter1 = arcpy.GetParameterAsText(0)
dataInputParameter1 = '{"Location":[{"address":"6188 Maddox Blvd, Chincoteague, VA, 23336","distance":"12" }]}'

        
# ADD ADDITIONAL VARIABLE FOR RETURN RESULT  
varResult = "true"
    
textDict = json.loads(dataInputParameter1)

langs = []
       
for SpecificLocation in textDict['Location']:
    varsearchAddress =  SpecificLocation["address"]
    varsearchDistance =  SpecificLocation["distance"]

    print ("address: " + varsearchAddress)
    print ("distance: " + varsearchDistance)
    print ""

    # LOCAL VARIABLES
    coordinates = ""
    countList = []   
    countyField = "FIPS2"  
    countList2 = []
    finalList = []
    txt_list = ""
    txt_list2 = ""
    data = ""

    # GEOCODERS
    geoCodeUrl = r"https://gismaps.vita.virginia.gov/arcgis/rest/services/Geocoding/VGIN_Composite_Locator/GeocodeServer/findAddressCandidates"

    searchAddress = varsearchAddress
    searchDistance = varsearchDistance

    # CONVERT MILE TO METERS FOR SPATIAL BUFFER
    distanceMeters = int(varsearchDistance) * 1609.34

    def singleAdressGeocode(address, geoCodeUrl, outSR = "4326"):
    #def singleAdressGeocode(address, geoCodeUrl, outSR = "26917"):
        ##clean up the address for url encoding
        address = address.replace(" ", "+")
        address = address.replace(",", "%3B")

        #send address to geocode service
        lookup = requests.get(geoCodeUrl + "?SingleLine=" + address + "&outSR=" + outSR + "&maxLocations=1&f=pjson")
        data = lookup.json()
        if data["candidates"]:
          ##coords = data["candidates"][0]["location"]
          ##return coords
          return data
        else:
          ##no results
          return "Address not geocoded: " + address

    # DEFINE VARIABLES FOR THE ADDRESS AND X AND Y VALUES
    geocodeResult = singleAdressGeocode(searchAddress, geoCodeUrl)
    addressResult = geocodeResult["candidates"][0]["address"]
    coordinateX = geocodeResult["candidates"][0]["location"]["x"]
    coordinateY = geocodeResult["candidates"][0]["location"]["y"]
    strcoordinateX = str(coordinateX)
    strcoordinateY = str(coordinateY)


    arcpy.env.workspace = r"E:/ArcGISProjects/CountySelection/xxxxx.sde"

    # CREATE THE POINT FROM THE X AND Y
    point = arcpy.Point(coordinateX, coordinateY)
    ptGeometry = arcpy.PointGeometry(point)
    
    # DEFINE THE STATE FEATURE CLASS
    arcpy.MakeFeatureLayer_management(r"E:/ArcGISProjects/CountySelection/xxxxx.sde/Boundary", "state_lyr")

    # BUFFER THE STATE LAYER AND PUT IT IN MEMORY
    arcpy.Buffer_analysis('state_lyr', 'IN_MEMORY/BufferGeom', distanceMeters)

    # TEST TO SEE IF POINT IS INSIDE THE BUFFERED GEOMETRY
    results = arcpy.SelectLayerByLocation_management(ptGeometry, 'INTERSECT', 'IN_MEMORY/BufferGeom')

    #Test if there are any polygons points:  
    if int(arcpy.GetCount_management('IN_MEMORY/BufferGeom').getOutput(0))==0:  
        raise Exception("Selected points do not overlap any polygons.")  
        print "Selected points do not overlap any polygons."
        varResult = "true"‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
25 Replies
jaykapalczynski
Honored Contributor

Last questions.....Can you do?

distance = 25 + " Mile"
point.with-distance-of(polygon, distance)‍‍
0 Kudos
DanPatterson
MVP Esteemed Contributor

Just explore the arcpy geometry objects.  If you are using arcmap (see Blake's link) then you will see the generic geometry properties and methods.  For other geometry objects, you can examine those as well.  Similar information is available for ArcGIS Pro which uses python 3 (currently 3.6.12 )


... sort of retired...
jaykapalczynski
Honored Contributor

Dan Patterson

was I correct in the logic about reading in my geometry for the polygon....

Yea I am using python here....trying to publish this out to a service but its not letting me.....keeps blowing up saying Error in creating the definition file

0 Kudos
RandyBurton
MVP Alum

Here is some test code for geometry methods .within and .distanceTo.

import arcpy

# IMPORTANT: all geometry should be in same projection

KMile = 1609.34

def ptGeo(geoResult):
    ptGeometry = arcpy.PointGeometry(arcpy.Point(
        geocodeResult["candidates"][0]["location"]["x"],
        geocodeResult["candidates"][0]["location"]["y"]),
        arcpy.SpatialReference(geocodeResult["spatialReference"]["wkid"]))
    return ptGeometry

# feature class contains only one feature, polygon outline of Virginia in UTM Zone 17N (spatial ref: 26917)
va = r'C:\Path\to\Test.gdb\Virginia'

# easy way to get Virginia geometry 
vaGeo = arcpy.da.SearchCursor(va,"SHAPE@").next()[0]

# paired address point and geocoding results
points = [['1520 Split Oak Ln, Henrico, VA, 23229',
           {u'candidates': [{u'attributes': {}, u'score': 100, u'location': {u'y': 4167142.8740720972, u'x': 805193.5919402852}, u'address': u'1520 SPLIT OAK LN, HENRICO, VA, 23229'}], u'spatialReference': {u'wkid': 26917, u'latestWkid': 26917}}],
          ['439 Rockford St, Mt Airy, NC 27030',
           {u'candidates': [{u'attributes': {}, u'score': 100, u'location': {u'y': 4038943.928205234, u'x': 533895.533587685}, u'address': u'27030'}], u'spatialReference': {u'wkid': 26917, u'latestWkid': 26917}}]
          ]

for address, geocodeResult in points:
    pt = ptGeo(geocodeResult)
    if pt.within(vaGeo):
        print("Location is within Virginia: {}".format(address))
    else:
        miles = pt.distanceTo(vaGeo)/KMile
        print("Location is {:.1f} miles outside Virginia: {}".format(miles, address))

# output
# Location is within Virginia: 1520 Split Oak Ln, Henrico, VA, 23229
# Location is 4.3 miles outside Virginia: 439 Rockford St, Mt Airy, NC 27030

''' methods:
within (second_geometry, {relation})
Indicates if the base geometry is within the comparison geometry. 

distanceTo (other)	
Returns the minimum distance between two geometries. If the geometries intersect, the minimum distance is 0.
Both geometries must have the same projection.
'''
jaykapalczynski
Honored Contributor

Randy Burton‌ Thanks Randy...I really do like this approach and learning as I go....much appreciated.  Think this is the approach I am going to take...modified your example a bit and added some more validation and it is working great..

0 Kudos
RandyBurton
MVP Alum

Another approach...

The Spatial Join (Analysis) will compute distances in certain cases.  As a test, a Virginia layer (containing only one feature - an outline of the state in UTM Zone 17N) was opened in Desktop 10.5.  I geocoded the two addresses previously used with some others that lie in and outside the state.  I pasted code into the Python window, ran the spatial join and had the following results:

The DistanceTo field was generated by the spatial join.  A distance of 0 indicates the point is inside the state.  If it lies outside, a  distance in meters is calculated.  This distance can be compared to the buffer in the Distance field.  If the DistanceTo is less than the Distance buffer, the point is within the buffer.  A select layer by attribute can use this in a where clause.

Here's the test code, which may need modification depending on your settings or if run outside ArcMap. In_memory was used in creating the points and spatialJoin layers, so they should be deleted at some point.

import arcpy, requests, json

#              ['searchAddress', searchDistance (in miles)],
addressList = [['1520 Split Oak Ln, Henrico, VA, 23229', 25],
               ['6188 Maddox Blvd, Chincoteague, VA, 23336', 12],
               ['318 Market St, Peterstown, WV 24963', 1],
               ['109 Federal St, Rich Creek, VA 24147',1],
               ['439 Rockford St, Mt Airy, NC 27030',2] ]

KMile = 1609.34

def singleAdressGeocode(address):
    geoCodeUrl = r"https://gismaps.vita.virginia.gov/arcgis/rest/services/Geocoding/VGIN_Composite_Locator/GeocodeServer/findAddressCandidates"
    outSR = "26917"
    r = requests.get(geoCodeUrl, params={'SingleLine': address, 'outSR' : outSR, 'maxLocations': 1, 'f': 'pjson'}) # requests will url encode params
    data = r.json()
    if data["candidates"]:
        return data
    else:
        return None # indicates error
    
def insertResult(address, distance, result, cursor):
    geoAddress = geocodeResult["candidates"][0]["address"]
    score = geocodeResult["candidates"][0]["score"]
    pointX = geocodeResult["candidates"][0]["location"]["x"]
    pointY = geocodeResult["candidates"][0]["location"]["y"]
    ptGeometry = arcpy.PointGeometry(arcpy.Point(
        geocodeResult["candidates"][0]["location"]["x"],
        geocodeResult["candidates"][0]["location"]["y"]),
        arcpy.SpatialReference(geocodeResult["spatialReference"]["wkid"]))
    cursor.insertRow([ptGeometry, address, geoAddress, score, pointX, pointY, distance])

# create in_memory feature to save address points
fc = arcpy.CreateFeatureclass_management("in_memory","points","POINT", spatial_reference=arcpy.SpatialReference(26917))
arcpy.AddField_management(fc,'Address','TEXT',"#","#",100,'Address',"NULLABLE","NON_REQUIRED")
arcpy.AddField_management(fc,'geoAddress','TEXT',"#","#",100,'Geocoded Address',"NULLABLE","NON_REQUIRED")
arcpy.AddField_management(fc,'geoScore','Double',"#","#","#",'Geocode Score',"NULLABLE","NON_REQUIRED")
arcpy.AddField_management(fc,'POINT_X','Double',"#","#","#",'Point X',"NULLABLE","NON_REQUIRED")
arcpy.AddField_management(fc,'POINT_Y','Double',"#","#","#",'Point Y',"NULLABLE","NON_REQUIRED")
arcpy.AddField_management(fc,'Distance','Double',"#","#","#",'Distance',"NULLABLE","NON_REQUIRED")

# start an insert cursor
cursor = arcpy.da.InsertCursor(fc,["SHAPE@", 'Address', 'geoAddress', 'geoScore', 'POINT_X', 'POINT_Y', 'Distance'])

# loop through addresses and insert in points feature
for address, distance in addressList:
    geocodeResult = singleAdressGeocode(address)
    if geocodeResult:
        insertResult(address, int(distance)*KMile, geocodeResult, cursor)
    else:
        print('ERROR: Address ({}) did not geocode'.format(address))

del cursor # done with insert cursor

# a spatial join can be used to compute distances
# Virginia layer contains only 1 feature, an outline of VA using NAD_1983_UTM_Zone_17N
# field_mapping is optional, only fields in points feature were used
arcpy.SpatialJoin_analysis(target_features="points",
                           join_features="Virginia",
                           out_feature_class="in_memory/spatialJoin",
                           join_operation="JOIN_ONE_TO_ONE",
                           join_type="KEEP_ALL",
                           field_mapping='''
Address "Address" true true false 100 Text 0 0 ,First,#,points,Address,-1,-1;
geoAddress "Geocoded Address" true true false 100 Text 0 0 ,First,#,points,geoAddress,-1,-1;
geoScore "Geocode Score" true true false 0 Double 0 0 ,First,#,points,geoScore,-1,-1;
POINT_X "Point X" true true false 0 Double 0 0 ,First,#,points,POINT_X,-1,-1;
POINT_Y "Point Y" true true false 0 Double 0 0 ,First,#,points,POINT_Y,-1,-1;
Distance "Distance" true true false 0 Double 0 0 ,First,#,points,Distance,-1,-1''',
                           match_option="CLOSEST_GEODESIC",
                           distance_field_name="DistanceTo")

# a search cursor can print a results table
fields = ['Address', 'geoAddress', 'geoScore', 'POINT_X', 'POINT_Y', 'Distance', 'DistanceTo']
print('\t'.join(fields)+'\tLocation')

with arcpy.da.SearchCursor("spatialJoin", fields) as cursor:
    for a, g, s, x, y, d, dt in cursor:
        if dt == 0:
            inout = 'Inside Virginia'
        elif dt < d:
            inout = 'Within Buffer'
        else:
            inout = 'Outside Buffer'
        print('{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}'.format(a, g, s, x, y, d, dt, inout))

# delete the in_memory feature when done
# arcpy.Delete_management('in_memory/points')
# arcpy.Delete_management('in_memory/spatialJoin')

''' # address points can be selected where the distance is greater than the buffer
arcpy.SelectLayerByAttribute_management(in_layer_or_view="spatialJoin",
                                        selection_type="NEW_SELECTION",
                                        where_clause="DistanceTo > Distance")
'''

Hope this helps.