Select to view content in your preferred language

Test for XY inside Polygon

4412
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
RandyBurton
MVP Alum

Regarding testing the new layer:

I was unable to turn ptGeometry into a layer as you are attempting to do in line 2 of your last code segment.  However, I was able to use ptGeometry as the select_features parameter in Select Layer By Location.  To test the values being used for x and y, you can print the JSON containing the x and y values:

# x,y for 6188 Maddox Blvd, Chincoteague, VA, 23336
x = -77.434769 # longitude
y = 37.541290 # latitude

ptGeometry = arcpy.PointGeometry(arcpy.Point(x,y),arcpy.SpatialReference(4326)).projectAs(arcpy.SpatialReference(26917))

print(ptGeometry.JSON)
''' Output:
{"x":815010.20711794402,"y":4160898.4941948522,"spatialReference":{"wkid":26917,"latestWkid":26917}}
''' ‍‍‍‍‍‍‍‍‍

I am a bit confused with your state layer, buffer and distance.  Does the layer have only one feature, the state of Virginia?  Or does it contain other states or Virginia counties.  Is the buffer being used to see if the address is close to the border of Virginia, possibly in a nearby state?

0 Kudos
jaykapalczynski
Honored Contributor

Thanks for your reply randy burton‌

  • At one point I was taking my state layer and buffering it and then trying to see if the point was inside it.
  • Now I am simply using the SelectLayerByLocation and specifying a distance to search for the point location within a user specified distance of the state

I WANT to find if the point is in the State or within a specified distance of the state 

The issue I am running into is that I cannot use the ptGeometry as a parameter because I get the following error

If I use a Point in afeature class it works.....

Notice below in the code that changed the output to UTM Zone 17N.  That is the Projection of the State Layer being referenced and its in meters.

I need to be able to geocode an address, get the XY and then test if its in the state with the user specified distance to search slightly outside the boundary.

results = arcpy.SelectLayerByLocation_management(ptGeometry, 'WITHIN_A_DISTANCE', 'state_lyr', distance, selection_type='NEW_SELECTION')

This is what I get back when I run it.

address: 1520 Split Oak Ln, Henrico, VA, 23229
distance: 25

805193.59194
4167142.87407

{"x":805193.59194028517,"y":4167142.8740720972,"spatialReference":{"wkid":null}}

Traceback (most recent call last):
  File "E:\ArcGIS Server\ArcGIS Published MXDs\Projects\NWTL\Python\NWTL_In_Virginia.py", line 120, in <module>
    results = arcpy.SelectLayerByLocation_management(ptGeometry, 'WITHIN_A_DISTANCE', 'state_lyr', distance, selection_type='NEW_SELECTION')
  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).‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

# SNIP OUT THE CODE ABOVE THAT HANDLES THE PARAMETERS FOR ADDRESS AND DISTANCE

    distanceMeters = int(varsearchDistance) * 1609.34

    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(varsearchAddress, geoCodeUrl)
    addressResult = geocodeResult["candidates"][0]["address"]
    coordinateX = geocodeResult["candidates"][0]["location"]["x"]
    coordinateY = geocodeResult["candidates"][0]["location"]["y"]

# DEFINE THE STATE FEATURE CLASS
    arcpy.MakeFeatureLayer_management(r"E:/ArcGISProjects/CountySelection/xx.sde/Virginia_Prj_Meters", "state_lyr")

## SET THE DISTANCE PARAMETER
    distance = int(distanceMeters)

    ptGeometry = arcpy.PointGeometry(arcpy.Point(coordinateX,coordinateY))

    arcpy.MakeFeatureLayer_management(ptGeometry, "PointfromGeocode")

    print(ptGeometry.JSON)

results = arcpy.SelectLayerByLocation_management(ptGeometry, 'WITHIN_A_DISTANCE', 'state_lyr', distance, selection_type='NEW_SELECTION')
   


# If features matched criteria, write them to a new feature class
    matchcount = int(arcpy.GetCount_management(results)[0]) 

    if matchcount == 0:
        print('no features matched spatial and attribute criteria')
    else:
        print('found something')
        print ""
        varResult = "true"‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
RandyBurton
MVP Alum

The ptGeometry item should have a spatial reference; it should not be null (Line 7 in the first code block in your previous message).  I think the first parameter in SelectLayerByLocation should be a feature layer, and ptGeometry (which is only geometry for a point) is not a feature.  I have been experimenting with some code in the Python window of Desktop 10.5 that may help you get the geocoding address coordinates into point geometry and then into a feature.  The test code did create a point feature and it appears to be correctly located using the Streets basemap (and Google Maps).  Here's my test code:

# create the point geometry

address: 1520 Split Oak Ln, Henrico, VA, 23229
distance: 25

>>> geocodeResult
{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}}
>>> addressResult = geocodeResult["candidates"][0]["address"]
>>> addressResult
u'1520 SPLIT OAK LN, HENRICO, VA, 23229'
>>> ptGeometry = arcpy.PointGeometry(arcpy.Point(
...     geocodeResult["candidates"][0]["location"]["x"],
...     geocodeResult["candidates"][0]["location"]["y"]),
...     arcpy.SpatialReference(geocodeResult["spatialReference"]["wkid"]))
...     
>>> ptGeometry.JSON
u'{"x":805193.59194028517,"y":4167142.8740720972,"spatialReference":{"wkid":26917,"latestWkid":26917}}'
>>>

# create an in_memory feature class to hold the point

>>> fc = arcpy.CreateFeatureclass_management("in_memory","points","POINT",
...                                          spatial_reference=arcpy.SpatialReference(geocodeResult["spatialReference"]["wkid"]))
... 
>>> fc
<Result 'in_memory\\points'>
>>> cursor = arcpy.da.InsertCursor(fc,["SHAPE@"])
>>> cursor.insertRow([ptGeometry])
1L
>>> del cursor

# delete the in_memory feature when done

>>> arcpy.Delete_management('in_memory/points')
<Result 'true'>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
jaykapalczynski
Honored Contributor

OK I am here again....still having issues with the first parameter ....Im trying to create the Feature and put it into memory as you suggested....but getting the same error where the parameters are incorrect.

I am not sure about this line in your example

fcPoint<Result 'in_memory\\points'>

ERROR:

Traceback (most recent call last):
File "E:\ArcGIS Server\ArcGIS Published MXDs\Projects\NWTL\Python\NWTL_In_Virginia.py", line 132, in <module>
results = arcpy.SelectLayerByLocation_management('IN_MEMORY\\points', 'WITHIN_A_DISTANCE', 'state_lyr', distance, selection_type='NEW_SELECTION')
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).

# USE THE GEOCODED LOCATION
    ptGeometry = arcpy.PointGeometry(arcpy.Point(\
        geocodeResult["candidates"][0]["location"]["x"],\
        geocodeResult["candidates"][0]["location"]["y"]),\
        arcpy.SpatialReference(geocodeResult["spatialReference"]["wkid"]))

    print(ptGeometry.JSON)
    print ""

    out_path = "IN_MEMORY"
    out_name = "points"
    geometry_type = "POINT"
    template = ""
    has_m = ""
    has_z = ""

    fcPoint = arcpy.CreateFeatureclass_management(out_path,out_name,geometry_type,template,has_m,has_z,\
        spatial_reference=arcpy.SpatialReference(geocodeResult["spatialReference"]["wkid"]))

    #fcPoint<Result 'in_memory\\points'>
    cursor = arcpy.da.InsertCursor(fcPoint,["SHAPE@"])
    cursor.insertRow([ptGeometry])
    del cursor
            
## TEST TO SEE IF POINT IS INSIDE THE BUFFERED GEOMETRY
    results = arcpy.SelectLayerByLocation_management('IN_MEMORY/points', 'WITHIN_A_DISTANCE', 'state_lyr', distance, selection_type='NEW_SELECTION')    

# If features matched criteria, write them to a new feature class
    matchcount = int(arcpy.GetCount_management(results)[0]) 

    if matchcount == 0:
        print('no features matched spatial and attribute criteria')
    else:
        print('found something')
        print ""
        varResult = "true"

# DELETE THE FEATURE WHEN DONE
    arcpy.Delete_management('IN_MEMORY/points')‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
jaykapalczynski
Honored Contributor

Think I might have it.....but not sure still testing...took the fcPoint we created and used MakeFeatureLayer and used that 

arcpy.MakeFeatureLayer_management(fcPoint, "PointfromGeocode")

results = arcpy.SelectLayerByLocation_management('PointfromGeocode', 'WITHIN', 'state_lyr')

0 Kudos
DanPatterson
MVP Esteemed Contributor

Why not make things easier and get the geometry objects for both.  Then you don't have to make featurelayers etc etc

As an example, a point geometry being tested against a polygon

You can project geometry if needed.  To find out more, just do a 'dir' on the objects

dir(pnt)
['JSON', 'WKB', 'WKT', ... snip...
, '_fromGeoJson', '_go', '_passthrough', '_repr_svg_', 'angleAndDistanceTo', 'area',
'boundary', 'buffer', 'centroid', 'clip', 'contains', 'convexHull', 'crosses', 'cut', 
'densify', 'difference', 'disjoint', 'distanceTo', 'equals', 'extent', 'firstPoint', 
'generalize', 'getArea', 'getGeohash', 'getLength', 'getPart', 'hasCurves', 
'hullRectangle', 'intersect', 'isMultipart', 'labelPoint', 'lastPoint', 'length', 
'length3D', 'measureOnLine', 'overlaps', 'partCount', 'pointCount', 
'pointFromAngleAndDistance', 'positionAlongLine', 'projectAs', 'queryPointAndDistance',
 'segmentAlongLine', 'snapToLine', 'spatialReference', 'symmetricDifference', 
'toCoordString', 'touches', 'trueCentroid', 'type', 'union', 'within']


# ---- ditto
dir(poly)
['JSON', 'WKB', 'WKT', '__add__',... snip ..., 'within']

# ---- pnt.within(poly)  # ---- point in polygon checks

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

Dan

I am interested in this approach...is there a better example?  Trying to get a better understanding of this.

0 Kudos
BlakeTerhune
MVP Regular Contributor
0 Kudos
DanPatterson
MVP Esteemed Contributor

As Blake suggests, since you have a mix of geometry origins (json, arcpy), and your query is essentially simple, it makes less sense now to be creating a bunch of featureclasses and/or featurelayers.  Just get the geometry at its root (that is, arcpy geometry objects) and do the check.  You basically need a boolean... if True, then do something, if False, move on 


... sort of retired...
0 Kudos
jaykapalczynski
Honored Contributor

OK I think I am understanding now....so I would just read in my geometry....

Do I have this correct....

  1. read in my Polygons coordinates instead of the squares XYs
  2. create the array from list of coords
  3. create my point with my XY coordinates
  4. Test and point.within(polygon) will  return true if within and false if outside

# READ IN THE COORDINATES OF MY POLY INSTEAD OF THE SQUARE BELOW
    coords = [[0, 0], [0, 2], [2, 2], [2, 0], [0, 0]]

# CREATE ARRAY OF THOSE COORDS
    polygon = arcpy.Polygon(arcpy.Array([arcpy.Point(*a) for a in coords]), sr)
    print "polygon area:", polygon.area

# READ IN MY POINT GEOMETRY
    point = arcpy.PointGeometry(arcpy.Point(X Coord, Y Coord), sr)

# TEST IF INSIDE WILL RETURN TRUE IF YES AND FALSE IF NO?
    print "point within polygon:", point.within(polygon)
0 Kudos