# Calculate Circle Geometry Based On Point for Map Service

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.

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 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))

longitude = float(argv)
latitude = float(argv)
miles = int(argv)

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']
``````

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']

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)

# 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:
closest = distance_miles
closest = result
``````

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