Test for XY inside Polygon

3362
25
Jump to solution
09-01-2020 03:05 PM
jaykapalczynski
Frequent 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
1 Solution

Accepted Solutions
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'>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

View solution in original post

0 Kudos
25 Replies
DanPatterson
MVP Esteemed Contributor

Select Layer By Location (Data Management)—ArcGIS Pro | Documentation 

you need to make a layer for your pntGeometry


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

OK that removed the error and an oversight by me.....thanks for the eyes....

For some reason I am getting true for any address....I think my code to test if its in or not is not working correctly.

If I put in an address for Wyoming it returns true....I set the variable to false at the beginning but for some reason I get a count from the results variable either way....

Any ideas where this test for the result is wrong?

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

    #Test if there are any polygons points:  
    if int(arcpy.GetCount_management(results).getOutput(0))==0:  
        #raise Exception("Selected points do not overlap any polygons.")  
        print "Selected points do not overlap any polygons."
        varResult = "true"
0 Kudos
jaykapalczynski
Frequent Contributor

Woke up at 3am and it dawned on my that I might be having issues with the projection of the Feature Class and the Point created from the Geocode...think once I fix that it might work....testing that this morning.

0 Kudos
jaykapalczynski
Frequent Contributor

NAH that did not work....but at least the polygon state layer is WGS84 and the point being geocoded is output in WGS 84 so thats on the same page but I never return any features.....

0 Kudos
jaykapalczynski
Frequent Contributor

dan patterson‌ OK I got it working but it seems the issue is around the point feature that I am creating from the Geocoded results...

Any thoughts? The data is in WGS84

If I point this to a point Feature Class as such this WORKS

arcpy.MakeFeatureLayer_management(r"E:/ArcGISProjects/CountySelection/xx.sde/Test", "testPoint")

results = arcpy.SelectLayerByLocation_management('testPoint', 'WITHIN', 'IN_MEMORY/BufferGeom') ‍‍‍‍‍‍

But I need to create a Point reference from the Geocoded XY this FAILS

# CREATE THE POINT FROM THE X AND Y
point = arcpy.Point(coordinateX, coordinateY)
ptGeometry = arcpy.PointGeometry(point)

# DEFINE THE POINT FEATURE CLASS
arcpy.MakeFeatureLayer_management(ptGeometry, 'point_lyr')

results = arcpy.SelectLayerByLocation_management('point_lyr', 'WITHIN', 'IN_MEMORY/BufferGeom') ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
RandyBurton
MVP Alum

I assume BufferGeom is the polygon layer, so perhaps something like:

x = -75.369027 # longitude
y = 37.935822 # latitude

# WGS 1984: 4326
# GCS North American 1983: 4269 ## use the factory/authority code for BufferGeom layer
ptGeometry = arcpy.PointGeometry(arcpy.Point(x,y),arcpy.SpatialReference(4326)).projectAs(arcpy.SpatialReference(4269))

results = arcpy.SelectLayerByLocation_management('IN_MEMORY/BufferGeom', 'CONTAINS', ptGeometry, selection_type="NEW_SELECTION")‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
jaykapalczynski
Frequent Contributor

randy burton‌  Thanks for your input...

If I do the below I am getting this error thrown back at me....

Maybe I should project all my data and use a different coord system?

Note:  My polygon layer (state of VA) is:  GCS_WGS_1984 WKID: 4326 Authority: EPSG

ERROR:

Traceback (most recent call last):
File "E:\ArcGIS Server\ArcGIS Published MXDs\Projects\NWTL\Python\NWTL_In_Virginia.py", line 109, in <module>
results = arcpy.SelectLayerByLocation_management('IN_MEMORY/BufferGeom', 'CONTAINS', ptGeometry, 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).

varsearchDistance = This value comes from a Parameter being sent to the .py script via a JSON string.  It is in miles and I then convert ot meters.  

## CONVERT MILE TO METERS FOR SPATIAL BUFFER
distanceMeters = int(varsearchDistance) * 1609.34
    
arcpy.MakeFeatureLayer_management(r"E:/ArcGISProjects/CountySelection/xx.sde/Virginia", "state_lyr")

    x = -75.369027 # longitude
    y = 37.935822 # latitude

    # WGS 1984: 4326
    # GCS North American 1983: 4269 ## use the factory/authority code for BufferGeom layer
    ptGeometry = arcpy.PointGeometry(arcpy.Point(x,y),arcpy.SpatialReference(4326)).projectAs(arcpy.SpatialReference(4269))
    
    ## BUFFER THE STATE LAYER AND PUT IT IN MEMORY
    arcpy.Buffer_analysis('state_lyr', 'IN_MEMORY/BufferGeom', distanceMeters)
    results = arcpy.SelectLayerByLocation_management('IN_MEMORY/BufferGeom', 'CONTAINS', ptGeometry, selection_type="NEW_SELECTION")‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
jaykapalczynski
Frequent Contributor

Think I might have it....Using the MakeFeatureLayer to take the ptGeometry and make a feature layer....

ptGeometry = arcpy.PointGeometry(arcpy.Point(coordinateX,coordinateY),arcpy.SpatialReference(4326)).projectAs(arcpy.SpatialReference(4269))

arcpy.MakeFeatureLayer_management(ptGeometry, "PointfromGeocode")

Also changed the Coord system to match that of my State Layer (webaux sphere)

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

For some reason not finding the XY coordinates keep telling me no matches...

    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
jaykapalczynski
Frequent Contributor

OK ... Puzzled here....I can get this to work if I reference the point that is in a Feature Class but

I cannot get it to work with the result from my Geocoded address and resulting XY

My original address that is geocoded is WGS84 I convert that to UTM Nad83 Zone 17N (26917)

The state State layer is stored in UTM Nad83 Zone17N (26917) Meters

The input from the parameters is Miles and I am converting that distance to Meters.

I know the distance thing is working because I can test that against a Point that is outside the state in a Feature Class and increase that distance and it eventually finds it...

Any thoughts as to why the Geocoded result is not being found...I am converting this to the Appropriate Spatial Reference.....ugggggg

This Does NOT work

    ## CONVERT MILE TO METERS FOR SPATIAL BUFFER
    # Value received from input parameter in JSON string
    distanceMeters = int(varsearchDistance) * 1609.34
   arcpy.MakeFeatureLayer_management(r"E:/ArcGISProjects/Selection/xx.sde/Virginia_UTM_Meters", "state_lyr")

    x = -77.434769 # longitude
    y = 37.541290 # latitude

    # convert from WGS84 to UTM zone 17N
    ptGeometry = arcpy.PointGeometry(arcpy.Point(x,y),arcpy.SpatialReference(4326)).projectAs(arcpy.SpatialReference(26917))

    arcpy.MakeFeatureLayer_management(ptGeometry, "PointfromGeocode")

    distance = int(distanceMeters)
       
    results = arcpy.SelectLayerByLocation_management('PointfromGeocode', 'WITHIN_A_DISTANCE', 'state_lyr', distance, selection_type='NEW_SELECTION')  ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

This works - only real difference is that I am pointing to a Feature Class point instead of usign the XY to create the point geometry

    ## CONVERT MILE TO METERS FOR SPATIAL BUFFER
    # Value received from input parameter in JSON string
    distanceMeters = int(varsearchDistance) * 1609.34
   arcpy.MakeFeatureLayer_management(r"E:/ArcGISProjects/Selection/xx.sde/Virginia_UTM_Meters", "state_lyr")

arcpy.MakeFeatureLayer_management(r"E:/ArcGISProjects/Selection/xx.sde/TestPoint", "testPoint")

    distance = int(distanceMeters)
       
    results = arcpy.SelectLayerByLocation_management('testPoint', 'WITHIN_A_DISTANCE', 'state_lyr', distance, selection_type='NEW_SELECTION')  ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

HOW CAN I Test this new Layer?  I tried the below but I dont get anything returned.

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

    print "start"

    fields = ['SHAPE@XY']

    with arcpy.da.SearchCursor('PointfromGeocode', fields) as cursor:
        for row in cursor:
            print(u'{0}'.format(row[0]))

    print "end of first"
0 Kudos