Find all your DIY electronics in the MakerShed. 3D Printing, Kits, Arduino, Raspberry Pi, Books & more!

pointpoly_20080302.jpg

update: As readers noted, it’s not the 0 degrees longitude that’s the problem, it’s at 180 degrees where you could encounter issues. I’ve also escaped the gt and lt symbols. Sorry about that.

I spent the weekend participating in the F1 Website Challenge, a coding marathon in which competing teams each produce a mythical man-month’s worth of web site for a worthy non-profit organization—all in the space of 24 hours.

One of the challenges my team faced during development was finding an efficient way for detecting a particular service region for a given address. Our client, Metro Meals on Wheels, has a number of different regions in which they deliver meals, with each region being served by a particular Meals on Wheels organization. These regions are defined by non-overlapping complex polygons. It’s not as simple as a normal vendor search, where you return the nearest location to the requested address. Instead you need to search a database of polygons to find the particular one which intersects the address location.

One of my teammates, Mark Seemann, ended up providing a fairly elegant solution to the problem, and was able to implement it in a simple SQL query. To find out if a point intersects a polygon, it’s as simple as drawing a vector from the point and seeing how many line segments of the polygon it crosses. If the number is even, it’s outside the polygon. If it’s odd, you have an intersection.

So let’s say you have a polygon database which has a row for each line segment of a polygon. You can quickly pull all segments that intersect a vector pointing directly east of your geocoded location like this:

SELECT poly_id, segment_id
    FROM segments
    WHERE ( lnga > thelng OR lngb > thelng )
          AND ( (lata > thelat AND latb < thelat )
              OR (latb > thelat AND lata < thelat ) )

That will return you a list of all line segments that you would cross if you walked directly east from the location at [thelat,thelng] (yes, this assumes you don’t cross 180 degrees longitude). To determine the polygon (or polygons) that intersect our address, it’s as simple as grouping by poly and returning all rows that have an odd number of matches:

SELECT poly_id, COUNT(segment_id) AS segment_count
    FROM segments
    WHERE ( lnga > thelng OR lngb > thelng )
          AND ( (lata > thelat AND latb < thelat )
              OR (latb > thelat AND lata < thelat ) )
          AND segment_count%2 = 1
    GROUP BY poly_id

Of course, the world isn’t flat, though I’ve treated it this way for simplicity. If you wanted this to work for all cases, you’d need to limit your search to a particular distance and translate the coordinates so that the search didn’t cross 180 degrees longitude.


Related

Comments

  1. Ed Davies says:

    Some test data for you to consider:

    lata = lnga = 1.0
    latb = lngb = 10.0
    thelat = 3.0
    thelng = 8.0

    The point is south and east of the line yet:

    (lngb > thelng) AND
    (latb > thelat AND lata

    Also, why do you think this would break as you cross the prime meridian (0° longitude) assuming a straightforward convention such as negative longitudes to represent points to the west? As somebody born and brought up within sight of the Greenwich meridian I’m a bit sensitive about this sort of thing, though it hasn’t stopped me releasing code with (at least) one bug in it which was triggered by data points on the meridian (and would have been by points on the equator, too).

  2. NoseyNick says:

    It might be confused by crossing 180degrees longitude – approximately the international dateline, but I don’t think it’ll be confused by crossing the 0 meridian.

  3. Steve says:

    Just use Oracle Spatial, SQL Server 2008, or PostGIS…it’s all been done before.

  4. Ed Davies says:

    Ug. I previewed and edited so that the less-than symbol was HTML escaped properly and it’s still broken in the display. The less-thans were also broken in my feed reader (Liferea), for what it’s worth.

    I agree with NoseyNick on the 180° and 0° thing.

  5. Ed Davies says:

    Now I think of it, you need to consider the case where one of the polygon nodes is at exactly the same latitude as the point under test.

  6. Jason Striegel says:

    Sorry about that. The gt and lt signs are encoded correctly now.

    Ed and NoseyNick are correct. The issue is when you as you near 180 degrees. There are segments that you’d want to test against, but as longitude increases past 180, it drops to -180 and you fail the “WHERE ( lnga > thelng OR lngb > thelng )” part of the test.

  7. valery says:

    Well, it’s a good effort … at re-inventing the wheel. PostGIS is a good example what has been done.

  8. Joel Byrnes says:

    Sure this might be available in PostGIS, Oracle Spatial, etc, but what if you don’t use Postgres and don’t have Spatial (only available in Oracle Enterprise)? This SQL looks to be compatible across all SQL databases, and is very easy to implement. I might find a use for this, and I thank the poster for their time and effort, rather than calling it a waste of time.

  9. Brainy Smurf says:

    I agree with Joel, you haters. “Just use something else” isn’t always an option. A nice, easy to implement solution that works across several DB’s.

  10. JB says:

    Great Idea! How could you loop this thru several points?

  11. G Moon says:

    This approach only tells you if the point is inside or outside the rectangle formed by the vector coordinates and not if it is E or W of the actual vector. You have to check with something like:

    SELECT poly_id, segment_id FROM segments WHERE same as before with the additional check
    AND ( thelng

    Also it does not work for all polygons as you do not check when edges (or vertices) align with your ‘walking’ vector so your counts need to reflect other cases to give you correct results. Assuming your points always line inside or outside and not on a polygon edge or vertices, I expect you will get ok results for your purpose. Spatial DBMS’s will have the additional checks as well as some form of spatial indices to allow large numbers of polygons to be efficiently checked.

  12. G Moon says:

    Lovely, messed up the posting even though it looked okay in preview. The additional check can be done with something like:

    … AND (thelng <= ((thelat – lata) * (lngb – lnga) / (latb – lata))) …

    Of course there are special cases for more complex situations so the approach you have works only for simple polygons and points definitely inside or outside.

  13. JB says:

    G Moon,

    Thanks for the addition to the post. I have found that several cases where one of the polygon nodes is at exactly the same latitude as the point under test. How could I adjust the code to account for this problem?

  14. ryan beckham says:

    i noticed you are using a polygon database.
    what field type are you using for the polygon points, and when you say each row contains a line segment that does mean each row contains two points of the polygon?
    thanks

  15. ryan beckham says:

    i noticed you are using a polygon database.
    what field type are you using for the polygon points,……… and when you say each row contains a line segment does that mean each row contains two points of the polygon?
    thanks

  16. Jason Striegel says:

    Ryan,

    That’s correct. Each row contains the id of the polygon and a lat/lon for each segment. In my pseudo-code above, there’s a segment table that has: segment_id int, poly_id int, lata float, lona float, latb float, lonb float.

    If you have a triangle with points a, b and c. There will be a segment entry for ab, bc, and ca. The poly id would be the same for all, and would point to an entry in a polygon table. The polygon table would then have the information about that feature, like its name and anything else associated with the polygon.

  17. Borna Segulin says:

    After digging through web I didn’t found complete solution so I’m posting my function here:

    suppose we have table with points of polygon:
    CREATE TABLE [tblArea](
    [pointID] [int] IDENTITY(1,1) NOT NULL,
    [polyID] [int] NOT NULL,
    [pointOrderNum] [smallint] NOT NULL,
    [latitude] [decimal](12, 8) NOT NULL,
    [longitude] [decimal](12, 8) NOT NULL
    )

    This function returns whether point is in poly:

    CREATE FUNCTION [dbo].[PointInPoly]
    (
    @latitude DECIMAL(12,8),
    @longitude DECIMAL(12,8),
    @polyID INT
    )
    RETURNS BIT
    AS
    BEGIN
    DECLARE @pointCount INT;
    SET @pointCount = (SELECT COUNT(*) FROM tblArea WHERE polyID=@polyID);

    RETURN (
    SELECT COUNT(*)
    FROM tblArea p1, tblArea p2
    WHERE p1.polyID=@polyID
    AND p2.polyID=@polyID
    AND p1.pointOrderNum = (CASE WHEN p2.pointOrderNum+1 >= @pointCount THEN 0 ELSE p2.pointOrderNum+1 END)
    AND ( (p1.longitude > @longitude AND p2.longitude @longitude) )
    AND ( @longitude * (P1.latitude – P2.latitude)/(P1.longitude – P2.longitude) ) + ( P1.latitude – ((P1.latitude – P2.latitude)/(P1.longitude – P2.longitude) * P1.longitude) ) > @latitude
    ) % 2;

    END

    and it works ok on all types of polygons.

  18. Borna Segulin says:

    After digging through web I didn’t found complete solution so I’m posting my function here:

    suppose we have table with points of polygon:
    CREATE TABLE [tblArea](
    [pointID] [int] IDENTITY(1,1) NOT NULL,
    [polyID] [int] NOT NULL,
    [pointOrderNum] [smallint] NOT NULL,
    [latitude] [decimal](12, 8) NOT NULL,
    [longitude] [decimal](12, 8) NOT NULL
    )

    This function returns whether point is in poly:

    CREATE FUNCTION [dbo].[PointInPoly]
    (
    @latitude DECIMAL(12,8),
    @longitude DECIMAL(12,8),
    @polyID INT
    )
    RETURNS BIT
    AS
    BEGIN
    DECLARE @pointCount INT;
    SET @pointCount = (SELECT COUNT(*) FROM tblArea WHERE polyID=@polyID);

    RETURN (
    SELECT COUNT(*)
    FROM tblArea p1, tblArea p2
    WHERE p1.polyID=@polyID
    AND p2.polyID=@polyID
    AND p1.pointOrderNum = (CASE WHEN p2.pointOrderNum+1 >= @pointCount THEN 0 ELSE p2.pointOrderNum+1 END)
    AND ( (p1.longitude > @longitude AND p2.longitude @longitude) )
    AND ( @longitude * (P1.latitude – P2.latitude)/(P1.longitude – P2.longitude) ) + ( P1.latitude – ((P1.latitude – P2.latitude)/(P1.longitude – P2.longitude) * P1.longitude) ) > @latitude
    ) % 2;

    END

    and it works ok on all types of polygons.

  19. Borna Segulin says:

    After digging through web I didn’t found complete solution so I’m posting my function here:

    suppose we have table with points of polygon:
    CREATE TABLE [tblArea](
    [pointID] [int] IDENTITY(1,1) NOT NULL,
    [polyID] [int] NOT NULL,
    [pointOrderNum] [smallint] NOT NULL,
    [latitude] [decimal](12, 8) NOT NULL,
    [longitude] [decimal](12, 8) NOT NULL
    )

    This function returns whether point is in poly:

    CREATE FUNCTION [dbo].[PointInPoly]
    (
    @latitude DECIMAL(12,8),
    @longitude DECIMAL(12,8),
    @polyID INT
    )
    RETURNS BIT
    AS
    BEGIN
    DECLARE @pointCount INT;
    SET @pointCount = (SELECT COUNT(*) FROM tblArea WHERE polyID=@polyID);

    RETURN (
    SELECT COUNT(*)
    FROM tblArea p1, tblArea p2
    WHERE p1.polyID=@polyID
    AND p2.polyID=@polyID
    AND p1.pointOrderNum = (CASE WHEN p2.pointOrderNum+1 >= @pointCount THEN 0 ELSE p2.pointOrderNum+1 END)
    AND ( (p1.longitude > @longitude AND p2.longitude @longitude) )
    AND ( @longitude * (P1.latitude – P2.latitude)/(P1.longitude – P2.longitude) ) + ( P1.latitude – ((P1.latitude – P2.latitude)/(P1.longitude – P2.longitude) * P1.longitude) ) > @latitude
    ) % 2;

    END

    and it works ok on all types of polygons.

  20. Stan says:

    in the following line :
    AND ( (p1.longitude > @longitude AND p2.longitude @longitude) )

    Is there a symbol between p2.longitude and @longitude?

    I m looking for the same query.

    Stan