How to get the intersecting area of two polygons in MySQL

I was playing around with spatial features of MySQL this weekend and stumbled into a problem where I was looking for the area of two rectangles that overlap. MySQL provides a function to check if they overlap, but no function to extract the region that overlaps.

I’ve never written a stored routine in MySQL before, so I decided it would be a good exercise to try making one. As you can see the function is pretty straightforward and it assumes you are working with rectangles, but other than that it does what it is supposed to.
You pass the function 2 polygons (e.g. Intersection(a.poly,b.poly)), and it returns the intersecting area as a new polygon.

DROP FUNCTION IF EXISTS Intersection;
DELIMITER //
CREATE FUNCTION Intersection(g1 POLYGON, g2 POLYGON) RETURNS POLYGON DETERMINISTIC
BEGIN
  DECLARE minx1, minx2, maxx1, maxx2 DOUBLE;
  DECLARE miny1, miny2, maxy1, maxy2 DOUBLE;
  DECLARE miny,  maxy,  minx,  maxx DOUBLE;
  SET minx1 =    LEAST(X(PointN(ExteriorRing(g1),1)),X(PointN(ExteriorRing(g1),2)),X(PointN(ExteriorRing(g1),3)),X(PointN(ExteriorRing(g1),4)));
  SET maxx1 = GREATEST(X(PointN(ExteriorRing(g1),1)),X(PointN(ExteriorRing(g1),2)),X(PointN(ExteriorRing(g1),3)),X(PointN(ExteriorRing(g1),4)));
  SET miny1 =    LEAST(Y(PointN(ExteriorRing(g1),1)),Y(PointN(ExteriorRing(g1),2)),Y(PointN(ExteriorRing(g1),3)),Y(PointN(ExteriorRing(g1),4)));
  SET maxy1 = GREATEST(Y(PointN(ExteriorRing(g1),1)),Y(PointN(ExteriorRing(g1),2)),Y(PointN(ExteriorRing(g1),3)),Y(PointN(ExteriorRing(g1),4)));
  SET minx2 =    LEAST(X(PointN(ExteriorRing(g2),1)),X(PointN(ExteriorRing(g2),2)),X(PointN(ExteriorRing(g2),3)),X(PointN(ExteriorRing(g2),4)));
  SET maxx2 = GREATEST(X(PointN(ExteriorRing(g2),1)),X(PointN(ExteriorRing(g2),2)),X(PointN(ExteriorRing(g2),3)),X(PointN(ExteriorRing(g2),4)));
  SET miny2 =    LEAST(Y(PointN(ExteriorRing(g2),1)),Y(PointN(ExteriorRing(g2),2)),Y(PointN(ExteriorRing(g2),3)),Y(PointN(ExteriorRing(g2),4)));
  SET maxy2 = GREATEST(Y(PointN(ExteriorRing(g2),1)),Y(PointN(ExteriorRing(g2),2)),Y(PointN(ExteriorRing(g2),3)),Y(PointN(ExteriorRing(g2),4)));
  SET minx = IF ( minx2 > minx1 , minx2 , minx1 );
  SET maxx = IF ( maxx2 > maxx1 , maxx1 , maxx2 );
  SET miny = IF ( miny2 > miny1 , miny2 , miny1 );
  SET maxy = IF ( maxy2 > maxy1 , maxy1 , maxy2 );
  RETURN POLYGON(LINESTRING(POINT(minx,miny),POINT(maxx,miny), POINT(maxx,maxy), POINT(minx,maxy),POINT(minx,miny)));
END //
DELIMITER ;

DROP FUNCTION IF EXISTS Intersection;

DELIMITER //

CREATE FUNCTION Intersection(g1 POLYGON, g2 POLYGON) RETURNS POLYGON DETERMINISTIC

BEGIN

DECLARE minx1, minx2, maxx1, maxx2 DOUBLE;

DECLARE miny1, miny2, maxy1, maxy2 DOUBLE;

DECLARE miny, maxy, minx, maxx DOUBLE;

SET minx1 = LEAST(X(PointN(ExteriorRing(g1),1)),X(PointN(ExteriorRing(g1),2)),X(PointN(ExteriorRing(g1),3)),X(PointN(ExteriorRing(g1),4)));

SET maxx1 = GREATEST(X(PointN(ExteriorRing(g1),1)),X(PointN(ExteriorRing(g1),2)),X(PointN(ExteriorRing(g1),3)),X(PointN(ExteriorRing(g1),4)));

SET miny1 = LEAST(Y(PointN(ExteriorRing(g1),1)),Y(PointN(ExteriorRing(g1),2)),Y(PointN(ExteriorRing(g1),3)),Y(PointN(ExteriorRing(g1),4)));

SET maxy1 = GREATEST(Y(PointN(ExteriorRing(g1),1)),Y(PointN(ExteriorRing(g1),2)),Y(PointN(ExteriorRing(g1),3)),Y(PointN(ExteriorRing(g1),4)));

SET minx2 = LEAST(X(PointN(ExteriorRing(g2),1)),X(PointN(ExteriorRing(g2),2)),X(PointN(ExteriorRing(g2),3)),X(PointN(ExteriorRing(g2),4)));

SET maxx2 = GREATEST(X(PointN(ExteriorRing(g2),1)),X(PointN(ExteriorRing(g2),2)),X(PointN(ExteriorRing(g2),3)),X(PointN(ExteriorRing(g2),4)));

SET miny2 = LEAST(Y(PointN(ExteriorRing(g2),1)),Y(PointN(ExteriorRing(g2),2)),Y(PointN(ExteriorRing(g2),3)),Y(PointN(ExteriorRing(g2),4)));

SET maxy2 = GREATEST(Y(PointN(ExteriorRing(g2),1)),Y(PointN(ExteriorRing(g2),2)),Y(PointN(ExteriorRing(g2),3)),Y(PointN(ExteriorRing(g2),4)));

SET minx = IF ( minx2 > minx1 , minx2 , minx1 );

SET maxx = IF ( maxx2 > maxx1 , maxx1 , maxx2 );

SET miny = IF ( miny2 > miny1 , miny2 , miny1 );

SET maxy = IF ( maxy2 > maxy1 , maxy1 , maxy2 );

RETURN POLYGON(LINESTRING(POINT(minx,miny),POINT(maxx,miny), POINT(maxx,maxy), POINT(minx,maxy),POINT(minx,miny)));

END //

DELIMITER ;

Example comparing some rectangles in 2 tables using the function:

CREATE TABLE `poly1` (
  `id` int(10) unsigned NOT NULL,
  `poly` polygon NOT NULL,
  PRIMARY KEY (`id`),
  SPATIAL KEY `poly` (`poly`)
) ENGINE=MyISAM;

CREATE TABLE `poly2` (
  `id` int(10) unsigned NOT NULL,
  `poly` polygon NOT NULL,
  PRIMARY KEY (`id`),
  SPATIAL KEY `poly` (`poly`)
) ENGINE=MyISAM;

-- Insert some test data into the tables
insert into poly1 (id,poly) VALUES (1, GEOMFROMWKB(POLYGON(LINESTRING(POINT(1,0),POINT(4,0), POINT(4,2), POINT(1,2),POINT(1,0)))));
insert into poly1 (id,poly) VALUES (2, GEOMFROMWKB(POLYGON(LINESTRING(POINT(6,0),POINT(10,0),POINT(10,2),POINT(6,2),POINT(6,0)))));
insert into poly1 (id,poly) VALUES (3, GEOMFROMWKB(POLYGON(LINESTRING(POINT(10,10),POINT(20,10),POINT(20,5),POINT(10,5),POINT(10,10)))));

insert into poly2 (id,poly) VALUES (1, GEOMFROMWKB(POLYGON(LINESTRING(POINT(1,0),POINT(4,0), POINT(4,2), POINT(1,2),POINT(1,0)))));
insert into poly2 (id,poly) VALUES (2, GEOMFROMWKB(POLYGON(LINESTRING(POINT(3,0),POINT(7,0), POINT(7,2), POINT(3,2),POINT(3,0)))));
insert into poly2 (id,poly) VALUES (3, GEOMFROMWKB(POLYGON(LINESTRING(POINT(5,0),POINT(7,0), POINT(7,2), POINT(5,2),POINT(5,0)))));
insert into poly2 (id,poly) VALUES (4, GEOMFROMWKB(POLYGON(LINESTRING(POINT(8,0),POINT(10,0),POINT(10,2),POINT(8,2),POINT(8,0)))));
insert into poly2 (id,poly) VALUES (5, GEOMFROMWKB(POLYGON(LINESTRING(POINT(3,0),POINT(6,0), POINT(6,2), POINT(3,2),POINT(3,0)))));
insert into poly2 (id,poly) VALUES (6, GEOMFROMWKB(POLYGON(LINESTRING(POINT(5,0),POINT(11,0),POINT(11,2),POINT(5,2),POINT(5,0)))));
insert into poly2 (id,poly) VALUES (7, GEOMFROMWKB(POLYGON(LINESTRING(POINT(12,8),POINT(18,8),POINT(18,4),POINT(12,4),POINT(12,8)))));

SELECT
 a.id as 'poly1',
 b.id as 'poly2',
 cast(100/Area(a.poly)*Area(Intersection(a.poly,b.poly)) as decimal(4,2)) as '% of poly2 in poly1'
FROM
 poly1 a,
 poly2 b
WHERE
 MBROverlaps(b.poly,a.poly);

CREATE TABLE `poly1` (

`id` int(10) unsigned NOT NULL,

`poly` polygon NOT NULL,

PRIMARY KEY (`id`),

SPATIAL KEY `poly` (`poly`)

) ENGINE=MyISAM;

CREATE TABLE `poly2` (

`id` int(10) unsigned NOT NULL,

`poly` polygon NOT NULL,

PRIMARY KEY (`id`),

SPATIAL KEY `poly` (`poly`)

) ENGINE=MyISAM;

-- Insert some test data into the tables

insert into poly1 (id,poly) VALUES (1, GEOMFROMWKB(POLYGON(LINESTRING(POINT(1,0),POINT(4,0), POINT(4,2), POINT(1,2),POINT(1,0)))));

insert into poly1 (id,poly) VALUES (2, GEOMFROMWKB(POLYGON(LINESTRING(POINT(6,0),POINT(10,0),POINT(10,2),POINT(6,2),POINT(6,0)))));

insert into poly1 (id,poly) VALUES (3, GEOMFROMWKB(POLYGON(LINESTRING(POINT(10,10),POINT(20,10),POINT(20,5),POINT(10,5),POINT(10,10)))));

insert into poly2 (id,poly) VALUES (1, GEOMFROMWKB(POLYGON(LINESTRING(POINT(1,0),POINT(4,0), POINT(4,2), POINT(1,2),POINT(1,0)))));

insert into poly2 (id,poly) VALUES (2, GEOMFROMWKB(POLYGON(LINESTRING(POINT(3,0),POINT(7,0), POINT(7,2), POINT(3,2),POINT(3,0)))));

insert into poly2 (id,poly) VALUES (3, GEOMFROMWKB(POLYGON(LINESTRING(POINT(5,0),POINT(7,0), POINT(7,2), POINT(5,2),POINT(5,0)))));

insert into poly2 (id,poly) VALUES (4, GEOMFROMWKB(POLYGON(LINESTRING(POINT(8,0),POINT(10,0),POINT(10,2),POINT(8,2),POINT(8,0)))));

insert into poly2 (id,poly) VALUES (5, GEOMFROMWKB(POLYGON(LINESTRING(POINT(3,0),POINT(6,0), POINT(6,2), POINT(3,2),POINT(3,0)))));

insert into poly2 (id,poly) VALUES (6, GEOMFROMWKB(POLYGON(LINESTRING(POINT(5,0),POINT(11,0),POINT(11,2),POINT(5,2),POINT(5,0)))));

insert into poly2 (id,poly) VALUES (7, GEOMFROMWKB(POLYGON(LINESTRING(POINT(12,8),POINT(18,8),POINT(18,4),POINT(12,4),POINT(12,8)))));

SELECT

a.id as 'poly1',

b.id as 'poly2',

cast(100/Area(a.poly)*Area(Intersection(a.poly,b.poly)) as decimal(4,2)) as '% of poly2 in poly1'

FROM

poly1 a,

poly2 b

WHERE

MBROverlaps(b.poly,a.poly);

Result:

poly1 poly2 % of poly2 in poly1

1 2 33.33

1 5 33.33

2 2 25.00

2 3 25.00

3 7 36.00