14

We need to figure out a quick and fairly accurate method for point-in-polygon for lat/long values and polygons over google maps. After some research - came across some posts about mysql geometric extensions, and did implement that too - 

SELECT id, Contains( PolyFromText( 'POLYGON(".$polygonpath.")' ) , PointFromText( concat( \"POINT(\", latitude, \" \", longitude, \")\" ) ) ) AS
            CONTAINS
FROM tbl_points

That did not however work with polygons made up of a large number of points :(

After some more research - came across a standard algorithm called the Ray-casting algorithm but before trying developing a query for that in MySQL, wanted to take my chances if someone had already been through that or came across a useful link which shows how to implement the algorithm in MySQL / SQL-server.

So, cutting it short - question is:

Can anyone please provide the MySQL/SQL-server implementation of Ray casting algorithm?

Additional detail:

  • Polygons are either of concave, convex or complex.
  • Targeting quick execution over 100% accuracy.
Johan
  • 74,508
  • 24
  • 191
  • 319
zarun
  • 910
  • 9
  • 26
  • 4
    When I looked at geospatial extensions in MySQL about a year ago, they were beta quality at best. I ended up using PostgreSQL for my geo database. I have never heard of anyone trying to implement ray casting in a database, however... – Eric J. Sep 27 '11 at 18:26
  • @EricJ. - thank you for your response. I wish I could use postGIS for this problem .. but can't do so because its only a (small probably) part of a huge system. :( – zarun Sep 27 '11 at 18:32
  • MySQL has cos/sin/tan functions, does that help you? link: http://dev.mysql.com/doc/refman/5.0/en/mathematical-functions.html – Johan Oct 03 '11 at 20:49

7 Answers7

18

The following function (MYSQL version of Raycasting algorithm) rocked my world :

CREATE FUNCTION myWithin(p POINT, poly POLYGON) RETURNS INT(1) DETERMINISTIC 
BEGIN 
DECLARE n INT DEFAULT 0; 
DECLARE pX DECIMAL(9,6); 
DECLARE pY DECIMAL(9,6); 
DECLARE ls LINESTRING; 
DECLARE poly1 POINT; 
DECLARE poly1X DECIMAL(9,6); 
DECLARE poly1Y DECIMAL(9,6); 
DECLARE poly2 POINT; 
DECLARE poly2X DECIMAL(9,6); 
DECLARE poly2Y DECIMAL(9,6); 
DECLARE i INT DEFAULT 0; 
DECLARE result INT(1) DEFAULT 0; 
SET pX = X(p); 
SET pY = Y(p); 
SET ls = ExteriorRing(poly); 
SET poly2 = EndPoint(ls); 
SET poly2X = X(poly2); 
SET poly2Y = Y(poly2); 
SET n = NumPoints(ls); 
WHILE i<n DO 
SET poly1 = PointN(ls, (i+1)); 
SET poly1X = X(poly1); 
SET poly1Y = Y(poly1); 
IF ( ( ( ( poly1X <= pX ) && ( pX < poly2X ) ) || ( ( poly2X <= pX ) && ( pX < poly1X ) ) ) && ( pY > ( poly2Y - poly1Y ) * ( pX - poly1X ) / ( poly2X - poly1X ) + poly1Y ) ) THEN 
SET result = !result; 
END IF; 
SET poly2X = poly1X; 
SET poly2Y = poly1Y; 
SET i = i + 1; 
END WHILE; 
RETURN result; 
End;

Add 

  DELIMITER ;;

before the function as required. The usage for the function is:

 SELECT myWithin(point, polygon) as result;

where

 point  = Point(lat,lng) 
 polygon = Polygon(lat1 lng1, lat2 lng2, lat3 lng3, .... latn lngn, lat1 lng1)

Please note that the polygon ought to be closed (normally it is closed if you're retrieving a standard kml or googlemap data but just make sure it is - note lat1 lng1 set is repeated at the end)

I did not have points and polygons in my database as geometric fields, so I had to do something like: 

 Select myWithin(PointFromText( concat( "POINT(", latitude, " ", longitude, ")" ) ),PolyFromText( 'POLYGON((lat1 lng1, ..... latn lngn, lat1 lng1))' ) ) as result

I hope this might help someone.

zarun
  • 910
  • 9
  • 26
5

I would write a custom UDF that implements the ray-casting algorithm in C or Delphi or whatever high level tool you use:

Links for writing a UDF
Here's sourcecode for a MySQL gis implementation that looks up point on a sphere (use it as a template to see how to interact with the spatial datatypes in MySQL).
http://www.lenzg.net/archives/220-New-UDF-for-MySQL-5.1-provides-GIS-functions-distance_sphere-and-distance_spheroid.html

From the MySQL manual:
http://dev.mysql.com/doc/refman/5.0/en/adding-functions.html

UDF tutorial for MS Visual C++
http://rpbouman.blogspot.com/2007/09/creating-mysql-udfs-with-microsoft.html

UDF tutorial in Delphi:
Creating a UDF for MySQL in Delphi

Source-code regarding the ray-casting algorithm
Pseudo-code: http://rosettacode.org/wiki/Ray-casting_algorithm
Article in drDobbs (note the link to code at the top of the article): http://drdobbs.com/cpp/184409586
Delphi (actually FreePascal): http://www.cabiatl.com/mricro/raycast/

Community
  • 1
  • 1
Johan
  • 74,508
  • 24
  • 191
  • 319
  • +1 Just be aware that custom UDFs may break when upgrading MySQL – cEz Oct 03 '11 at 21:08
  • Thanks! A lot of useful links there. Going through each of them. Also, I was stuck in the problem when it came to this notification system which should notify users in a boundary when something of importance is posted within. So - users(points), event(points) and boundary(polygon) .. I figure it would be dead slow to pass hundreds of such boundary information to a function and then check if the incident point is within or not the boundary - then check all users again within such a valid boundary at a time and notify them ;( – zarun Oct 07 '11 at 20:36
  • 1
    @zarun the trick is to make sure you can use indexes. Index allow you to eliminate 99.99% of the search space. Thats how they speed things up. – Johan Oct 07 '11 at 20:48
  • @Johan - thanks man .. I acknowledge the input. But what I meant in my reply above was regarding "UDF that implements the ray-casting algorithm in C or Delphi or whatever high level tool you use" ---------- i.e., it would be slow to pass on values into a php function (I'm doing mysql/php for the problem above). But nevermind thanks. :) – zarun Oct 09 '11 at 06:41
  • You cannot write the UDF in php, it must **compile** into a DLL. That dll you attach to MySQL. This only leaves tools that generate native x86 code. – Johan Oct 09 '11 at 18:15
  • thx got it. I was just replyin to the original comment actually. Also, I've a solution now - got Raycasting implemented in MYSQL! Thanks for the resources! – zarun Oct 10 '11 at 06:20
2

Just in case, one MySQL function which accepts MULTIPOLYGON as an input: http://forums.mysql.com/read.php?23,286574,286574

DELIMITER $$

CREATE DEFINER=`root`@`localhost` FUNCTION `GISWithin`(pt POINT, mp MULTIPOLYGON) RETURNS int(1)
    DETERMINISTIC
BEGIN 

DECLARE str, xy TEXT; 
DECLARE x, y, p1x, p1y, p2x, p2y, m, xinters DECIMAL(16, 13) DEFAULT 0; 
DECLARE counter INT DEFAULT 0; 
DECLARE p, pb, pe INT DEFAULT 0; 

SELECT MBRWithin(pt, mp) INTO p; 

IF p != 1 OR ISNULL(p) THEN 
    RETURN p; 
END IF; 

SELECT X(pt), Y(pt), ASTEXT(mp) INTO x, y, str; 
SET str = REPLACE(str, 'POLYGON((',''); 
SET str = REPLACE(str, '))', ''); 
SET str = CONCAT(str, ','); 

SET pb = 1; 
SET pe = LOCATE(',', str); 
SET xy = SUBSTRING(str, pb, pe - pb); 
SET p = INSTR(xy, ' '); 
SET p1x = SUBSTRING(xy, 1, p - 1); 
SET p1y = SUBSTRING(xy, p + 1); 
SET str = CONCAT(str, xy, ','); 

WHILE pe > 0 DO 
    SET xy = SUBSTRING(str, pb, pe - pb); 
    SET p = INSTR(xy, ' '); 
    SET p2x = SUBSTRING(xy, 1, p - 1); 
    SET p2y = SUBSTRING(xy, p + 1); 

    IF p1y < p2y THEN SET m = p1y; ELSE SET m = p2y; END IF; 

    IF y > m THEN 
        IF p1y > p2y THEN SET m = p1y; ELSE SET m = p2y; END IF; 
        IF y <= m THEN 
            IF p1x > p2x THEN SET m = p1x; ELSE SET m = p2x; END IF; 
            IF x <= m THEN 
                IF p1y != p2y THEN 
                    SET xinters = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x; 
                END IF; 
                IF p1x = p2x OR x <= xinters THEN 
                    SET counter = counter + 1; 
                END IF; 
            END IF; 
        END IF; 
    END IF; 

    SET p1x = p2x; 
    SET p1y = p2y; 
    SET pb = pe + 1; 
    SET pe = LOCATE(',', str, pb); 
END WHILE; 

RETURN counter % 2; 

END
Igor
  • 2,619
  • 6
  • 24
  • 36
2

In reply to zarun function for finding lat/long within polygon.

I had a property table having lat/long information. So I had to get the records whose lat/long lies within polygon lats/longs (which I got from Google API). At first I was dumb how to use the Zarun function. So here is the solution query for it.

  • Table: properties
  • Fields: id, latitude, longitude, beds etc...
  • Query:
SELECT id
FROM properties
WHERE myWithin(
    PointFromText(concat( "POINT(", latitude, " ", longitude, ")")),
    PolyFromText('POLYGON((37.628134 -77.458334,37.629867 -77.449021,37.62324 -77.445416,37.622424 -77.457819,37.628134 -77.458334))' )
) = 1 limit 0,50;

Hope it will save time for dumbs like me ;)

M. Mimpen
  • 1,212
  • 16
  • 21
Nasir Hussain
  • 31
  • 2
1

I wanted to use the above mywithin stored procedure on a table of polygons so here are the commands to do just that.

After importing a shapefile containing polygons into mysql using ogr2ogr as follows

ogr2ogr -f "mysql" MYSQL:"mydbname,host=localhost,user=root,password=mypassword,port=3306" -nln "mytablename" -a_srs "EPSG:4326" /path/to/shapefile.shp

you can then use MBRwithin to prefilter your table and mywithin to finish as follows

DROP TEMPORARY TABLE IF EXISTS POSSIBLE_POLYS;
CREATE TEMPORARY TABLE POSSIBLE_POLYS(OGR_FID INT,SHAPE POLYGON);
INSERT INTO POSSIBLE_POLYS (OGR_FID, SHAPE) SELECT mytablename.OGR_FID,mytablename.SHAPE FROM mytablename WHERE MBRwithin(@testpoint,mytablename.SHAPE);

DROP TEMPORARY TABLE IF EXISTS DEFINITE_POLY;
CREATE TEMPORARY TABLE DEFINITE_POLY(OGR_FID INT,SHAPE POLYGON);
INSERT INTO DEFINITE_POLY (OGR_FID, SHAPE) SELECT POSSIBLE_POLYS.OGR_FID,POSSIBLE_POLYS.SHAPE FROM POSSIBLE_POLYS WHERE mywithin(@testpoint,POSSIBLE_POLYS.SHAPE);

where @testpoint is created, for example, from

SET @longitude=120;
SET @latitude=-30;
SET @testpoint =(PointFromText( concat( "POINT(", @longitude, " ", @latitude, ")" ) ));
Mr Purple
  • 2,325
  • 1
  • 18
  • 15
1

It is now a Spatial Extension as of MySQL5.6.1 and above. See function_st-contains at Docs.

Nasser Al-Wohaibi
  • 4,562
  • 2
  • 36
  • 28
  • Simply one line with a link is not an answer on StackOverflow. Please incorporate information - such as a summary and code examples - from the link into this post. – M. Mimpen Jan 21 '14 at 12:41
1

Here is a version that works with MULTIPOLYGONS (an adaptation of Zarun's one which only works for POLYGONS).

CREATE FUNCTION GISWithin(p POINT, multipoly MULTIPOLYGON) RETURNS INT(1) DETERMINISTIC 
BEGIN 
DECLARE n,i,m,x INT DEFAULT 0; 
DECLARE pX,pY,poly1X,poly1Y,poly2X,poly2Y DECIMAL(13,10); 
DECLARE ls LINESTRING; 
DECLARE poly MULTIPOLYGON;
DECLARE poly1,poly2 POINT; 
DECLARE result INT(1) DEFAULT 0; 
SET pX = X(p); 
SET pY = Y(p); 
SET m = NumGeometries(multipoly);
WHILE x<m DO 
SET poly = GeometryN(multipoly,x);
SET ls = ExteriorRing(poly); 
SET poly2 = EndPoint(ls); 
SET poly2X = X(poly2); 
SET poly2Y = Y(poly2); 
SET n = NumPoints(ls); 
WHILE i<n DO 
SET poly1 = PointN(ls, (i+1)); 
SET poly1X = X(poly1); 
SET poly1Y = Y(poly1); 
IF ( ( ( ( poly1X <= pX ) && ( pX < poly2X ) ) || ( ( poly2X <= pX ) && ( pX < poly1X ) ) ) && ( pY > ( poly2Y - poly1Y ) * ( pX - poly1X ) / ( poly2X - poly1X ) + poly1Y ) ) THEN 
SET result = !result; 
END IF; 
SET poly2X = poly1X; 
SET poly2Y = poly1Y; 
SET i = i + 1; 
END WHILE; 
SET x = x + 1;
END WHILE;
RETURN result; 
End;
daamsie
  • 1,539
  • 10
  • 15