Calculate Circle Geometry Based On Point for Map Service

Question

I have a Map Service that is currently just multiple boundary layers. Users query this Map Service with a lat/long to determine whether or not a Point falls within a Boundary. Currently they are using the ArcGIS Identify operation to do this procedure.

We've found cases where some points of data are falling outside of all boundaries. If the Point is X distance away from a boundary, the users would still like it to be returned. In cases like this, I was thinking of having the user calculate a Polygon Circle with X distance (X determined by their specific function, but let's say 5 miles) based upon their Lat/Long.

My thinking was to have them pass in the geometry for this Polygon Circle, rather than passing in the Point. If there is an intersection, return that boundary. I understand that this could lead to returning multiple results, but the users are fine with this. They basically care whether or not the Point falls within, or near, a given boundary.

What are some good methods for calculating the geometry needed for a Polygon Circle, given the Lat/Long, to be used in the Map Service Identify function?

What I've Tried: I have tried using the "Tolerance" parameter. However, I'm not sure I like this given that it's based upon Pixel tolerance (I understand why). I feel like this wouldn't be the best route to go.

I was told by ESRI (am waiting on a link) that there is an SoE that allows for a "Find Nearest" type function. I'm worried that adding in an SoE could cause blocking on future update for ArcGIS Server. E.X. - Map Services get fundamentally changed, which breaks the SoE, which prevents us from upgrading.

Answer 1

Discovered the answer using Shapely and PyProj.

Essentially what I needed to do was transform the input lat/long into a coordinate system with a unit of measurement in meters. I used World Mercator for this.

Once that was done, I was free to use Shapely to buffer the Point with X meters (I converted to miles). Once I was done, I retransformed the resulting Polygon into a geodetic system (NAD 1983 (2011)). From this, I was able to create a dictionary and passed the coordinates key into the Identify function of the Map Service. This returned all polygons that my polygon intersected.

from shapely.geometry import Point, Polygon, LinearRing, mapping
from shapely.ops import transform
from sys import argv
from math import cos, radians
from functools import partial
import pyproj
import json

def partial_projection(from_wkid, to_wkid):
    return partial(
        pyproj.transform,
        pyproj.Proj(init=from_wkid),
        pyproj.Proj(init=to_wkid))

# Args: Longitude Latitude Miles(Radius)
longitude = float(argv[1])
latitude = float(argv[2])
miles = int(argv[3])

point = Point(longitude, latitude)

# Project the inputted Point to Planar Coordinate System before buffering
to_6318 = partial_projection('epsg:6318', 'epsg:3395')

projected_point = transform(to_6318, point)

# Calculation: Convert Miles to Meters, Multiplied by 1 * the Cos(Latitude)
# We provide an additional buffer to accomodate for difference in geodetic vs. planar
buffer_area = (miles * 1609.34) * (1 / cos(radians(latitude)))

# Buffer returns a Polygon object.
buffered_point = projected_point.buffer(buffer_area, 32)

# Project the buffered point (now Polygon) back to a GeoDetic Coordinate System, 
# which is the coordinate system that is used our ArcGIS Feature Classes
to_3395 = partial_projection('epsg:3395', 'epsg:6318') 

# Map out the Point
coordinate_map = mapping(transform(to_3395, buffered_point))

coords = coordinate_map['coordinates'][0]

I took it one step further and requested geometries back from the query. I used these geometries in order to calculate the distance from the original lat/long to the nearest point on each returned polygon. From this I could determine which one was closest.

# Loop through result geometries and calculate distance from lat/long to outer point on boundary
for result in data['results']:
    coordinates = result['geometry']['rings'][0]

    polygon = Polygon(coordinates)

    # Projected Distance expects a LinearRing
    poly_ext = LinearRing(polygon.exterior.coords)

    # Project our Linear Ring to the Lat/Long.
    d = poly_ext.project(point)

    # Based upon the above distance, interpolate to grab the the point at which we meet.
    p = poly_ext.interpolate(d)
    closest_point_coord = list(p.coords)[0]

    # Get distance to closest Point. Use the planar lat/long and project 
    # our closest_point to planar to get a planar distance back.
    # This will return a distance in Meters. 
    # This probably isn't needed as the Geodetic distance is probably sufficient.
    distance = projected_point.distance(transform(to_6318, Point(closest_point_coord)))
    distance_miles = distance / 1609.34

    if distance_miles < closest[0]:
        closest[0] = distance_miles
        closest[1] = result

Credit to https://gis.stackexchange.com/a/142327/62130 for the planar distance calculation