I have city boundaries from the census loaded into an oracle database. I am trying to use the data to resolve the city and state for any latitude or longitude. I have sql statement to reduce the cities I need to check for an area, however, I had to write my own pl/SQL function to test if a point is with in city limits. Its a ray casting algorithm I got from the ADA version listed here http://rosettacode.org/wiki/Ray-casting_algorithm . The problem is when the poloygon is complex (i.e it has shapes with in the shape) the algorithm is returning invalid results. Can anyone suggest a way to get around this problem ?

Here is the Pl/SQL Polygon Function:

CREATE OR REPLACE PACKAGE  Polygons as
    type Point is record
    (
        X number, 
        Y number
    );
    type Point_List is table of Point;
    type poly_Segment is VARRAY(2) of Point;
    type Polygon is TABLE of poly_Segment;
    function Get_Point(x number,y number) return Point;
    function Create_Polygon (listOfPoints Point_List) return Polygon;
    function Is_Inside (pt2Check Point, shape2Check Polygon) return varchar2;
end Polygon;

/

CREATE OR REPLACE package body Polygons as

    EPSILON CONSTANT number := 0.00001;

    function Get_Point(x number,y number) return Point is
        pnt2Rtn Point;
    BEGIN
        pnt2Rtn.x := x;
        pnt2Rtn.y := y;
        return pnt2Rtn;
    END Get_Point;

    function Ray_Intersects_Segment(Who Point, locSec poly_Segment)  return  Boolean
    is
        The_Point Point := Who;
        Above Point;
        Below Point;
        M_Red   number;
        Red_Is_Infinity Boolean := False;
        M_Blue number;
        Blue_Is_Infinity Boolean := False;
    begin
        if locSec (1).Y < locSec (2).Y then
            Above := locSec (2);
            Below := locSec (1);
        else
            Above := locSec (1);
            Below := locSec (2);
        end if;
        if The_Point.Y = Above.Y or The_Point.Y = Below.Y then
            The_Point.Y := The_Point.Y + EPSILON;
        end if;
        if The_Point.Y < Below.Y or The_Point.Y > Above.Y then
            return False;
        elsif The_Point.X > Above.X and The_Point.X > Below.X then
            return False;
        elsif The_Point.X < Above.X and The_Point.X < Below.X then
            return True;
        else
            if Above.X <> Below.X then
                M_Red := (Above.Y - Below.Y) / (Above.X - Below.X);
            else
                Red_Is_Infinity := True;
            end if;
            if Below.X <> The_Point.X then
                M_Blue := (The_Point.Y - Below.Y) / (The_Point.X - Below.X);
            else
                Blue_Is_Infinity := True;
            end if;
            if Blue_Is_Infinity then
                return True;
            elsif Red_Is_Infinity then
                return False;
            elsif M_Blue >= M_Red then
                return True;
            else
                return False;
            end if;
        end if;
    end Ray_Intersects_Segment;

    function Create_Polygon (listOfPoints Point_List) return Polygon is
        polyRes Polygon := new Polygon();
        Side poly_Segment;
    begin
        polyRes.extend(listOfPoints.Count);
        for I in 1..listOfPoints.Count loop
             Side := new poly_Segment();
             Side.extend(2);
             Side (1) := listOfPoints (I);
             --connect connect the lines if its the last item connect back to the first
             if I = listOfPoints.COUNT then
                Side (2) := listOfPoints (listOfPoints.FIRST);
             else
                Side (2) := listOfPoints (I + 1);
             end if;

             polyRes(I) := Side;
        end loop;

        return polyRes;
    end Create_Polygon;

    function Is_Inside (pt2Check Point, shape2Check Polygon) return varchar2 is
        cnt4Intersect number := 0;
    begin
        for Side in 1..shape2Check.count loop
            if Ray_Intersects_Segment (pt2Check, shape2Check (Side)) then
            cnt4Intersect := cnt4Intersect + 1;
            end if;
        end loop;

        if cnt4Intersect mod 2 = 0 then
            return 'N';
        else
            return 'Y';
        end if;
    end Is_Inside;

end Polygons;
/

Here is a link to the shape data http://pastebin.com/yxrftXQn

I was testing to see if 33.66840,-87.59806 was with in the city (it is, but according to my function its not)

EDIT So after some research I have found the issue is actually how I am defining the polygons shapes. I needed to adjust the data feed.

EDIT Is it possible to keep the Polygon Package and have the body use SDO_RELATE and still gain efficiency (can some one provide example code ) ?

有帮助吗?

解决方案

The ray casting algorithm is not the problem it works ! The issue was I needed to separate the data points using the parts column which was not included(because it wasn't clear that it would be useful) . In short the data was bad, and I needed to correctly format the shapes (i also had WAY too much shape data, because of the way I was exporting the .shp file..) ..

Each shape should be its own Polygons, and checking if a shape is inside is as simple as checking that all the points are with in the shape...

I added a few new functions to the polygon in order to support what I needed to do..ideally, I'd like to switch to using the upstream functions but they aren't installed on the server I am using..

Here is the complete Polygon Package:

CREATE OR REPLACE PACKAGE  Polygons as

    type Point is record
    (
        X number, 
        Y number
    );

    type Point_List is table of Point;

    type poly_Segment is VARRAY(2) of Point;

    type Polygon is TABLE of poly_Segment;

    type Polygons is TABLE of Polygon;

    function Get_Point(x number,y number) return Point;

    function Create_Polygon (listOfPoints Point_List) return Polygon;

    function Is_Inside (pt2Check Point, shape2Check Polygon) return varchar2;


    function Is_Inside_Polygon (pt2Check Point, shape2Check Polygon,shapes2Include Polygons) return varchar2;

    function Is_Inside_Polygon (pt2Check Point, shape2Check Polygon,shapes2Exclude Polygons,shapes2Include Polygons) return varchar2;

end Polygons;


/
CREATE OR REPLACE package body Polygons as

    EPSILON CONSTANT number := 0.00001;

    function Get_Point(x number,y number) return Point is
        pnt2Rtn Point;
    BEGIN
        pnt2Rtn.x := x;
        pnt2Rtn.y := y;
        return pnt2Rtn;
    END Get_Point;

    function Ray_Intersects_Segment(Who Point, locSec poly_Segment)  return  Boolean
    is
        The_Point Point := Who;
        Above Point;
        Below Point;
        M_Red   number;
        Red_Is_Infinity Boolean := False;
        M_Blue number;
        Blue_Is_Infinity Boolean := False;
    begin
        if locSec (1).Y < locSec (2).Y then
            Above := locSec (2);
            Below := locSec (1);
        else
            Above := locSec (1);
            Below := locSec (2);
        end if;
        if The_Point.Y = Above.Y or The_Point.Y = Below.Y then
            The_Point.Y := The_Point.Y + EPSILON;
        end if;
        if The_Point.Y < Below.Y or The_Point.Y > Above.Y then
            return False;
        elsif The_Point.X > Above.X and The_Point.X > Below.X then
            return False;
        elsif The_Point.X < Above.X and The_Point.X < Below.X then
            return True;
        else
            if Above.X <> Below.X then
                M_Red := (Above.Y - Below.Y) / (Above.X - Below.X);
            else
                Red_Is_Infinity := True;
            end if;
            if Below.X <> The_Point.X then
                M_Blue := (The_Point.Y - Below.Y) / (The_Point.X - Below.X);
            else
                Blue_Is_Infinity := True;
            end if;
            if Blue_Is_Infinity then
                return True;
            elsif Red_Is_Infinity then
                return False;
            elsif M_Blue >= M_Red then
                return True;
            else
                return False;
            end if;
        end if;
    end Ray_Intersects_Segment;

    function Create_Polygon (listOfPoints Point_List) return Polygon is
        polyRes Polygon := new Polygon();
        Side poly_Segment;
    begin
        polyRes.extend(listOfPoints.Count);
        for I in 1..listOfPoints.Count loop
             Side := new poly_Segment();
             Side.extend(2);
             Side (1) := listOfPoints (I);
             --connect connect the lines if its the last item connect back to the first
             if I = listOfPoints.COUNT then
                Side (2) := listOfPoints (listOfPoints.FIRST);
             else
                Side (2) := listOfPoints (I + 1);
             end if;

             polyRes(I) := Side;
        end loop;

        return polyRes;
    end Create_Polygon;

    function Is_Inside (pt2Check Point, shape2Check Polygon) return varchar2 is
        cnt4Intersect number := 0;
    begin
        for Side in 1..shape2Check.count loop
            if Ray_Intersects_Segment (pt2Check, shape2Check (Side)) then
            cnt4Intersect := cnt4Intersect + 1;
            end if;
        end loop;

        if  cnt4Intersect mod 2 = 0  then
            --return to_char(shape2Check.count);
            return 'N';
        else
            return 'Y';
        end if;
    end Is_Inside;

    function Is_Inside_Polygon (pt2Check Point, shape2Check Polygon,shapes2Include Polygons) return varchar2 is
  begin
        return Is_Inside_Polygon(pt2Check,shape2Check,new Polygons(),shapes2Include);
    end Is_Inside_Polygon;

    function Is_Inside_Polygon (pt2Check Point, shape2Check Polygon,shapes2Exclude Polygons,shapes2Include Polygons) return varchar2 is
        foundInShape varchar2(4);
        currentSet number := 1;
        shape2BeExcluded varchar2(4);
    begin
        foundInShape := Is_Inside(pt2Check,shape2Check);
        IF foundInShape = 'N' THEN
            IF shapes2Include is not null THEN
                for currentSet in 1..shapes2Include.count loop
                    foundInShape := Is_Inside(pt2Check,shapes2Include(currentSet));
                    exit when foundInShape = 'Y';
                END Loop;
            END IF;
        END IF;
        IF foundInShape = 'Y' THEN
            IF shapes2Exclude is not null THEN
                for currentSet in 1..shapes2Exclude.count loop
                    shape2BeExcluded := Is_Inside(pt2Check,shapes2Exclude(currentSet));
                    IF shape2BeExcluded = 'Y' THEN
                        foundInShape := 'N';
                        exit;
                    END IF;
                END Loop;
            END IF;
        END IF;
        return foundInShape;
    end Is_Inside_Polygon;

end Polygons;
/
许可以下: CC-BY-SA归因
不隶属于 StackOverflow
scroll top