I was looking for a "point in polygon" function which doens't use simple MBR.

Finally I found a post from Adam Smith with a working MySQL stored function, based on the idea from

http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
But, the function was too slow, so I optimized it a little bit.

If someone does have a better alternative or a faster MySQL stored function, please let me know!!!

Here is the code...

DELIMITER //

CREATE 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;

DELIMITER ;

-- usage example:

CREATE TABLE `polygons` (

`id` int(11) NOT NULL auto_increment,

`polygon_data` polygon NOT NULL,

PRIMARY KEY (`id`),

SPATIAL KEY `polygon_data` (`polygon_data`(32))

) ENGINE=MyISAM DEFAULT CHARSET=latin1;

SET @x = -79.12345678901234;

SET @y = 43.12345678901234;

SET @point = CONCAT('POINT(',@x,' ',@y,')');

SELECT id FROM polygons WHERE GISWithin(GeomFromText(@point), polygon_data);

Edited 3 time(s). Last edit at 10/16/2009 11:46AM by Gert Serneels.