Webservice to determine if a Latitude / Longitude Point is inside a polygon or building

I spent some time trying to find an easy, do-it-yourself method to build a web-service that would take a Lat/Long as input and tell me if that point is contained within a polygon.  For my purposes I wanted to know if a point was within a building or not and a polygon will easily satisfy this request.

The Ray-casting algorithm is a method to determine just that.  Google maps and Esri have function calls that one can use within their API's, but those are client-side as you can only call such functions through javascript, JQuery or the like.  Thus if you want something on the server side you must write it yourself.

Luckily I found this rosettacode side that had the ray casting algorithm coded in several languages!  Their example uses pre-canned polygons and points.  I modified it to accept a lat/long from a user and return true if the point is within any of my polygons/buildings and false otherwise. The list of polys contain the Lat/Long points for each building.



My python code:

# Code was taken from http://rosettacode.org/wiki/Ray-casting_algorithm#Python
# uses simple ray-casting algoritm to see if point lies within a polygon
# 
# Adapted by Nick Gramsky to be used as a service to check if LatLong was inside  
#    a GIS polygon
from collections import namedtuple
from pprint import pprint as pp
import sys
 
Pt = namedtuple('Pt', 'x, y')               # Point
Edge = namedtuple('Edge', 'a, b')           # Polygon edge from a to b
Poly = namedtuple('Poly', 'name, edges')    # Polygon
 
_eps = 0.00001
_huge = sys.float_info.max
_tiny = sys.float_info.min
 
def rayintersectseg(p, edge):
    ''' takes a point p=Pt() and an edge of two endpoints a,b=Pt() of a line segment returns boolean
    '''
    a,b = edge
    if a.y > b.y:
        a,b = b,a
    if p.y == a.y or p.y == b.y:
        p = Pt(p.x, p.y + _eps)
 
    intersect = False
 
    if (p.y > b.y or p.y < a.y) or (
        p.x > max(a.x, b.x)):
        return False
 
    if p.x < min(a.x, b.x):
        intersect = True
    else:
        if abs(a.x - b.x) > _tiny:
            m_red = (b.y - a.y) / float(b.x - a.x)
        else:
            m_red = _huge
        if abs(a.x - p.x) > _tiny:
            m_blue = (p.y - a.y) / float(p.x - a.x)
        else:
            m_blue = _huge
        intersect = m_blue >= m_red
    return intersect
 
def _odd(x): return x%2 == 1
 
def ispointinside(p, poly):
    ln = len(poly)
    return _odd(sum(rayintersectseg(p, edge)
                    for edge in poly.edges ))
 
def polypp(poly):
    print "\n  Polygon(name='%s', edges=(" % poly.name
    print '   ', ',\n    '.join(str(e) for e in poly.edges) + '\n    ))'
 
if __name__ == '__main__':
    polys = [
      Poly(name='avwilliams', edges=(
        Edge(a=Pt(x=38.9913160, y=-76.937079), b=Pt(x=38.991333, y=-76.936119)),
        Edge(a=Pt(x=38.991333, y=-76.936119), b=Pt(x=38.990287, y=-76.936108)),
        Edge(a=Pt(x=38.990287, y=-76.936108), b=Pt(x=38.990278, y=-76.937057)),
        Edge(a=Pt(x=38.990278, y=-76.937057), b=Pt(x=38.990495,y=-76.937052)),
        Edge(a=Pt(x=38.990495,y=-76.937052), b=Pt(x=38.990499,y=-76.936424)),
        Edge(a=Pt(x=38.990499,y=-76.936424), b=Pt(x=38.991091,y=-76.93643)),
        Edge(a=Pt(x=38.991091,y=-76.93643), b=Pt(x=38.991104,y=-76.937079)),
        Edge(a=Pt(x=38.991104,y=-76.937079), b=Pt(x=38.9913160, y=-76.937079))
        )),
    ]  

    if len(sys.argv) != 3:
       print "Incorrect number of arguments.  Please submit a lat and a long...."

    userpoint = (Pt(x=float(sys.argv[1]), y=float(sys.argv[2])))

    testpoints = (Pt(x=float(sys.argv[1]), y=sys.argv[2]), Pt(x=38.990842, y=-76.93625),
                  Pt(x=38.9021466, y=-77), Pt(x=0, y=5),
                  Pt(x=10, y=5), Pt(x=8, y=5),
                  Pt(x=10, y=10))
 
    #print "\n TESTING WHETHER POINTS ARE WITHIN POLYGONS"
    inside = False
    for poly in polys:
        #polypp(poly)
        #print '   ', '\t'.join("%s: %s" % (p, ispointinside(p, poly))
        #                       for p in testpoints[:3])
        #if ispointinside(testpoints[1], poly):
        if ispointinside(userpoint, poly):
           inside = True

    print inside

And the corresponding php code:

Oh and this was done on a mac, hence the need to call out 'python' in the system call. Syntax for this call will vary by OS. Of course this example is using a hard-coded polygon within the python script. That's fine if you want to check only a few polygons/buildings but not a great way to code if you want to check many buildings. One could simply create a small database that stores polygons and iterate through each polygon. Using PostGres you can perform a bounds call to retrieve only polygons that contain the point within their minimin bounding rectangle and iterate through this python code. Additionally this could very well be coded in php but I just didn't have the time :) It would reduce the need to perform a system call and reduce the software dependencies....

 <PUBLIC SERVICE ANNOUNCEMENT> Yes, my php code IS very light and simple. One should ALWAYS clean user input prior to working with the parameters they provide. My example/prototype did not include that (BAD DEVELOPER) but this was an internal prototype only. </PUBLIC SERVICE ANNOUNCEMENT>