SQL Server Puzzle #5: Points Within Polygon

[gmmastros]

For my application, I need to create reports that show which points are contained within a polygon.

The coordinate system for this example is Latitude/Longitude. To define the polygon, there is a table MapcountiesPoints. In this table, there is a FeatureId (to identify the county), a point number field to identify the order in which points are put together, a latitude and longitude to identify the location.

CREATE TABLE [MapCountiesPoints] (
	[FeatureId] [int] NOT NULL ,
	[PointNumber] [int] NOT NULL ,
	[Latitude] [decimal](8, 5) NOT NULL ,
	[Longitude] [decimal](8, 5) NOT NULL ,
	CONSTRAINT [PK_MapCountiesPoints] PRIMARY KEY  CLUSTERED 
	(
		[FeatureId],
		[PointNumber]
	)  ON [PRIMARY] 
) ON [PRIMARY]

Sample data looks like this....

FeatureId   PointNumber Latitude   Longitude  
----------- ----------- ---------- ---------- 
6           1           40.22638   -75.60358
6           2           40.22615   -75.60386
6           3           40.22580   -75.60429
6           4           40.22574   -75.60436
6           5           40.22562   -75.60452
6           6           40.22544   -75.60476
6           7           40.22538   -75.60488
6           8           40.22534   -75.60493

A polygon is made up from a series of lines. To visualize the boundary, imagine a line drawn from point 1 to point 2, then point 2 to point 3, then point 3 to point 4, and so on. To close the boundary, the last point is connected to the first point.

Download the county data

There is another table, Location that identifies points we want to test. In this table, there is a unique identifier LocationId, and Latitude & longitude to identify the location.

Create TABLE [Location] (
	[LocationId] [int] NOT NULL ,
	[Longitude] [decimal](8, 5) NULL ,
	[Latitude] [decimal](8, 5) NULL ,
	CONSTRAINT [PK_Location_LocationId] PRIMARY KEY  CLUSTERED 
	(
		[LocationId]
	)  ON [PRIMARY] 
) ON [PRIMARY]

Download the location data

Here's the task....

Write a stored procedure that
1. Accepts FeatureId as an argument
2. Returns a list of LocationId's that are contained within the polygon identifed in the MapCountiesPoints table.

This functionality already works within my application, so I will humbly exclude myself from this puzzle.

Remember, fastest execution time wins. If you have any questions, just ask.

-George

Strong and bitter words indicate a weak cause. - Fortune cookie wisdom

 

[pwilson]

von, you write some quick stuff

   CREATE PROCEDURE County_pwilson
AS

Declare @StartTime DateTime
Declare @MidTime DateTime
Declare @EndTime DateTime
Declare @EndTime2 DateTime
Set @StartTime = GetDate()

DECLARE @maxpoint int
SELECT @maxpoint = MAX(PointNumber) FROM MapPoints

Set @EndTime = GetDate()


    SELECT COUNT(*) AS Total FROM (
    SELECT  LocationId  FROM 
    MapCountiesPoints MP1 INNER JOIN MapCountiesPoints MP2 ON MP1.PointNumber = (MP2.PointNumber + 1) %@maxpoint
    CROSS JOIN Location L 
    WHERE 
                    (
                        (MP1.Longitude <= L.Longitude AND L.Longitude < MP2.Longitude) OR
                        (MP2.Longitude <= L.Longitude AND L.Longitude < MP1.Longitude)
                   )
                   AND
                   (
                        L.Latitude < (MP2.Latitude - MP1.Latitude) * ((L.Longitude - MP1.Longitude) / (MP2.Longitude - MP1.Longitude)) + MP1.Latitude OR 
                        (
                            L.Latitude = (MP2.Latitude - MP1.Latitude) * ((L.Longitude - MP1.Longitude) / (MP2.Longitude - MP1.Longitude)) + MP1.Latitude AND (L.Latitude > MP1.Latitude AND L.Latitude < MP2.Latitude) OR (L.Latitude < MP1.Latitude AND L.Latitude > MP2.Latitude)
                        )
                    )
    GROUP BY LocationId 
    HAVING Count(*) % 2 = 1
    ) AS S1
Select @maxpoint, Convert(Decimal(10,5), DateDiff(Millisecond, @StartTime, @MidTime)) / 1000 AS Insert1, Convert(Decimal(10,5), DateDiff(Millisecond, @MidTime, @EndTime)) / 1000 AS Delete1, Convert(Decimal(10,5), DateDiff(Millisecond, @EndTime, Getdate())) / 1000 AS Select1, Convert(Decimal(10,5), DateDiff(Millisecond, @startTime, Getdate())) / 1000 AS Total1

;

Runs in about 11-13. 3754 records. Algorithm taken from

http://astronomy.swin.edu.au/~pbourke/geometry/insidepoly/

I ran some other tests on data similar to
Polygon

INSERT [MapCountiesPoints] ( FeatureId, PointNumber, Latitude, Longitude) VALUES (6,1,10.00000,10.00000)
INSERT [MapCountiesPoints] ( FeatureId, PointNumber, Latitude, Longitude) VALUES (6,2,10.00000,20.00000)
INSERT [MapCountiesPoints] ( FeatureId, PointNumber, Latitude, Longitude) VALUES (6,3,20.00000,20.00000)
INSERT [MapCountiesPoints] ( FeatureId, PointNumber, Latitude, Longitude) VALUES (6,4,25.00000,15.00000)
INSERT [MapCountiesPoints] ( FeatureId, PointNumber, Latitude, Longitude) VALUES (6,5,20.00000,10.00000)
INSERT [MapCountiesPoints] ( FeatureId, PointNumber, Latitude, Longitude) VALUES (6,6,10.00000,10.00000)



INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (1,9.00000,10.00000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (2,10.00000,9.00000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (3,10.00000,10.00000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (4,11.00000,10.00000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (5,10.00000,11.00000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (6,11.00000,11.00000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (7,12.50000,22.50000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (8,12.00000,22.50000)
INSERT [Location] ( LocationId, Longitude, Latitude) VALUES (9,13.00000,22.50000)

It finds points 3,4,5,6,7,9


NOTE: if the SP doesnt work, it is probably because I used X,Y instead of Long/Lat for my readability. I did a replace to get it back to long/lat

 

[vongrunt]

With that data my sproc crashes while attempting to ignore horizontal polygon lines... one of two exceptions I haven't covered yet.

------
"There's a man... He's bald and wears a short-sleeved shirt, and somehow he's very important to me. I think his name is Homer."
(Jack O'Neill, Stargate)

 

[gmmastros]

Here's my code. It produces the same count as vongrunt an pwilson.

Using the data supplied by pwilson, it produces the correct results.

ALTER  Procedure GetPointInPolygon_gmmastros
            @FeatureId Integer
As
SET NOCOUNT ON

Declare @MaxLongitude Decimal(8,5),
		@MinLongitude Decimal(8,5),
		@MaxLatitude Decimal(8,5),
		@MinLatitude Decimal(8,5)

Select	@MinLongitude = Min(Longitude),
		@MaxLongitude = Max(Longitude),
		@MinLatitude = Min(Latitude),
		@MaxLatitude = Max(Latitude)
From	MapCountiesPoints
Where	FeatureId = @FeatureId

Create
Table	#Candidates
		(
		LocationId Integer, 
		Longitude Decimal(8,5), 
		Latitude Decimal(8,5),
		Maybe Bit default 0,
		IntersectLatitude Decimal(8,5)
		)

Insert
Into	#Candidates(LocationId, Longitude, Latitude)
Select	LocationId,
		Longitude,
		Latitude
From	Location
Where	Longitude Between @MinLongitude and @MaxLongitude
		And Latitude Between @MinLatitude And @MaxLatitude

Declare @Lines
Table 	(
		LineNumber Integer,
		MinLongitude Decimal(8,5),
		MaxLongitude Decimal(8,5),
		MinLatitude Decimal(8,5),
		MaxLatitude Decimal(8,5)
		)

Insert
Into	@Lines
		(
		LineNumber,
		MinLongitude,
		MaxLongitude,
		MinLatitude,
		MaxLatitude
		)

Select  A.PointNumber,
        Case    When A.Longitude < B.Longitude
	            Then A.Longitude
	            Else B.Longitude
	            End As MinLongitude,
        Case    When A.Longitude > B.Longitude
                Then A.Longitude
                Else B.Longitude
                End As MaxLongitude,
        Case    When A.Longitude < B.Longitude
                Then A.Latitude
                Else B.Latitude
                End As MinLatitude,
        Case    When A.Longitude > B.Longitude
                Then A.Latitude
                Else B.Latitude
                End As MaxLatitude
From	MapCountiesPoints A
		Inner Join MapCountiesPoints B 
			On	A.FeatureId = B.FeatureId
			And A.PointNumber = B.PointNumber - 1
Where   A.FeatureId = @FeatureId

Update	C
Set		Maybe = 1
From 	#Candidates C
		Inner Join @Lines L
			On 	C.Longitude > L.MinLongitude 
			And C.Longitude < L.MaxLongitude
			And C.Latitude > L.MinLatitude 
			And C.Latitude < L.MaxLatitude

Declare @Locations Table(LocationId Integer)

-- 
-- This next part will get the points that are within the polygon, 
-- but only for points where the Latitude > Maxlatitude
--

Insert
Into	@Locations(LocationId)
Select	C.LocationId
From    #Candidates C
		Inner Join @Lines L
            On  C.Longitude > L.MinLongitude
            And C.Longitude <= L.MaxLongitude
            And C.Latitude > L.MaxLatitude
Where	C.Maybe = 0
Group By C.LocationId
Having Count(L.LineNumber) % 2 = 1

Insert
Into	@Locations(LocationId)
Select	C.LocationId
From    #Candidates C
		Inner Join @Lines L
            On  C.Longitude >= L.MinLongitude
            And C.Longitude < L.MaxLongitude
            And C.Latitude < L.MaxLatitude    --- I set this to equal because I think this is all we need to check
Where	C.Maybe = 0
Group By C.LocationId
Having Count(L.LineNumber) % 2 = 1

Insert 	Into @Locations(LocationId)
Select	A.LocationId
From	(
		Select	C.LocationId,
				L.LineNumber,
				C.Latitude,
				MinLatitude + (Longitude-MinLongitude) * (MaxLatitude-MinLatitude) / (MaxLongitude-MinLongitude) As IntersectLatitude
		From	#Candidates C
				Inner Join @Lines L
					On 	C.Longitude > L.MinLongitude 
					And C.Longitude < L.MaxLongitude
					And C.Latitude > L.MinLatitude 
					And C.Latitude < L.MaxLatitude
		Where	L.MaxLongitude <> L.MinLongitude
				And C.Maybe = 1
		) A
Where	A.Latitude <= A.IntersectLatitude
Group By A.LocationId
Having Count(A.LineNumber) %2 = 1

Select Distinct LocationId From @Locations

GO
SET QUOTED_IDENTIFIER OFF 
GO
SET ANSI_NULLS ON 
GO

Declare @S DateTime
Set @S = GetDate()

Exec GetPointInPolygon_gmmastros 6

Select DateDiff(Millisecond, @S, GetDate())

Running on the large dataset, it takes approximately 2 seconds.

-George

Strong and bitter words indicate a weak cause. - Fortune cookie wisdom

 

[SQLDenis]

I have an idea for a new puzzle, wanna hear it?

Denis The SQL Menace
SQL blog:

http://sqlservercode.blogspot.com/

Personal Blog:

http://otherthingsnow.blogspot.com/

 

[gmmastros]

I was planning on giving this one a couple more hours before declaring a winner. Right now, I'm leaning toward pwilson because his SP produces the correct results.

-George

Strong and bitter words indicate a weak cause. - Fortune cookie wisdom

 

[vongrunt]

So far I'd also vote for pwilson's code - if not for speed then for clear code and accuracy.

Btw. I'm just finishing Ultimate Sample Data (20-line polygon looking like Lego monster from Loch Ness) so let's delay results for #5 for day or two and proceed with #6 in the meantime.

------
"There's a man... He's bald and wears a short-sleeved shirt, and somehow he's very important to me. I think his name is Homer."
(Jack O'Neill, Stargate)

 

[vongrunt]

OK, here is sample data:

INSERT INTO [MapCountiesPoints] VALUES(6, 1, 1.00, 9.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 2, 5.00, 9.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 3, 6.00, 7.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 4, 4.00, 7.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 5, 2.00, 5.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 6, 4.00, 5.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 7, 2.00, 3.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 8, 3.00, 3.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 9, 3.00, 1.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 10, 2.00, 1.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 11, 2.00, 0.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 12, 6.00, 0.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 13, 6.00, -2.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 14, 5.00, -1.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 15, 4.00, -1.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 16, 3.00, -1.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 17, 0.00, -1.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 18, 0.00, 2.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 19, 0.00, 9.00)
INSERT INTO [MapCountiesPoints] VALUES(6, 20, 1.00, 9.00)

INSERT INTO [Location] VALUES(1, 1.00, 7.00)
INSERT INTO [Location] VALUES(2, 1.00, 10.00)
INSERT INTO [Location] VALUES(3, 1.00, 9.00)
INSERT INTO [Location] VALUES(4, 2.00, 7.00)
INSERT INTO [Location] VALUES(5, 2.00, 8.00)
INSERT INTO [Location] VALUES(6, 3.00, 10.00)
INSERT INTO [Location] VALUES(7, 5.00, 8.00)
INSERT INTO [Location] VALUES(8, 4.50000, 4.00)
INSERT INTO [Location] VALUES(9, 6.00, 4.00)
INSERT INTO [Location] VALUES(10, 3.00, 4.00)
INSERT INTO [Location] VALUES(11, 2.00, 4.00)
INSERT INTO [Location] VALUES(12, 3.00, 0.50000)
INSERT INTO [Location] VALUES(13, 2.00, -0.50000)
INSERT INTO [Location] VALUES(14, -2.00, -1.00)
INSERT INTO [Location] VALUES(15, 3.00, -1.00)
INSERT INTO [Location] VALUES(16, 4.00, -2.00)
INSERT INTO [Location] VALUES(17, 6.00, -2.00)
INSERT INTO [Location] VALUES(18, 6.00, -3.00)
INSERT INTO [Location] VALUES(19, 7.00, 0.00)
INSERT INTO [Location] VALUES(20, 6.00, 0.00)
INSERT INTO [Location] VALUES(21, 0.00, 0.00)

I guess it covers all extreme situations/exceptions. Polygon itself is not convex and lines are defined counter-clockwise (ccw) - some algorithms are sensitive to that order so beware. Wanna reverse order, simply do

UPDATE MapCountiesPoints SET PointNumber = 22-PointNumber WHERE FeatureID = 6

Locations (points) within polygon are:

1 3 4 5 7 10 11 13 15 17 20 21

Of those, the following locations fall exactly at polygon's boundaries (lines):

3 10 15 17 20 21

Needless to say, my sproc fails... what about yours?

------
"There's a man... He's bald and wears a short-sleeved shirt, and somehow he's very important to me. I think his name is Homer."
(Jack O'Neill, Stargate)

 

[vongrunt]

Aww... cr@p.

I forgot that Location table has swapped lat/long... use code from above but swap these values with:

update Location set Longitude = Latitude, Latitude = Longitude

Mea culpa .

------
"There's a man... He's bald and wears a short-sleeved shirt, and somehow he's very important to me. I think his name is Homer."
(Jack O'Neill, Stargate)

 

[gmmastros]

vongrunt, is this what your 'map' looks like? I wrote this to help me figure out why I am missing 1 point.

When you move the mouse, you can see the latitude & longitude for the position of the mouse.

When you hover over a point, you can see the FeatureId for the point. To run this, you will have to have the vb6 runtimes.

Point Within Polygon

I am missing point 17. If you view this as a 'loch ness monster' travelling from right to left, it's the point at the nose.


-George

Strong and bitter words indicate a weak cause. - Fortune cookie wisdom

 

[vongrunt]

Yup, that's the map. And even when we pass this test, there are two other simple tests to perform... and then we can announce code 100% safe.

Your sproc misses Location #17 (data precision problems?), pwilson's #3, #17 and #20.

My sproc cannot get the job done (runtime error), plus there are two fundamental flaws I'll have to fix... probably tomorrow.

------
"There's a man... He's bald and wears a short-sleeved shirt, and somehow he's very important to me. I think his name is Homer."
(Jack O'Neill, Stargate)

 

[gmmastros]

>> And even when we pass this test, there are two other simple tests to perform... and then we can announce code 100% safe.

Ok, I'm getting the right results now. What are the other 2 simple tests?



-George

Strong and bitter words indicate a weak cause. - Fortune cookie wisdom

 

[vongrunt]

1. Turn entire poly set clockwise (see UPDATE query below sample data, I think there should be value 21 instead of 22). Some routines are sensitive to that order but maybe shouldn't (depending on algorithm used).

2. Because polygon is represented as cyclical ordered set of points (first and last points are same), some routines are sensitive about initial poly point (PointNum=1). Try to put this point:

a) on "spike". Say, move poly point 2 from long./lat. = (9, 5) to (8, 0.5). That should also move points #5 and #7 out of polygon.

b) on beginning of horizontal line: move point 2 from (9, 5) to (8, 1). #5 and #7 should get excluded again.

c) in the middle of horizontal line, by making PointNum=18 (2, 0) initial point and shifting all other PointNums accordingly.

Of course - a) thru c) must give identical results.

------
"There's a man... He's bald and wears a short-sleeved shirt, and somehow he's very important to me. I think his name is Homer."
(Jack O'Neill, Stargate)