MySQL – Ryan Schulze

How to get the intersecting area of two polygons in MySQL

December 3, 2012September 27, 2014RyanGeometric Objects, Intersection, MySQL, polygon, spatial, stored procedure

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

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

How to build an efficient GeoIP SQL table

August 22, 2011April 10, 2012Ryan1 Commentbash, GeoIP, MySQL

This here is a very handy little script I threw together to generate a geoip.sql table for quickly determining which country a IP is from. I already hear you saying “Just convert the IP to an INT and use BETWEEN, how hard can it be”. And you are right, that works. And it may even be your easiest solution, but it just isn’t fast. And if you are planning on hammering the table with thousands of queries you are going to end up looking for something fast.

A while back I found a very interesting posting at www.jcole.us that described how to use Spacial Indexes together with MySQL’s GIS to speed up the queries. The posting has been online for a while and both it and the replies are worth reading.

All I did was make a small bash script to download the current “lite” version of GeoIP CSV file from maxmind.com, use the information from the posting to throw/transform it into a local database table and dump out a .sql file that can be easily imported into any other database. The script isn’t failproof though, it expects your user to be able to use mysql and have permission to create databases/tables and “load data local infile”.

generate_geoip_sql.sh

#!/bin/bash
#===============================================================================
#
#          FILE: generate_geoip_sql.sh
#
#         USAGE: ./generate_geoip_sql.sh
#
#   DESCRIPTION: Generates a optimized sql dump with geoipdata from maxmind.com
#
#         USAGE: import sql into a database
#
#  REQUIREMENTS: mysql database rights to create databases and LOAD DATA LOCAL INFILE
#
#       EXAMPLE:
#                 SELECT country_code,country_name
#                 FROM   geo_ip
#                 WHERE  MBRCONTAINS(ip_poly, POINTFROMWKB(POINT(INET_ATON('1.2.3.4'), 0)))
#
#                 SELECT country_code,country_name
#                 FROM   geo_ip
#                 WHERE  INET_ATON('1.2.3.4') 
#                 BETWEEN ip_from AND ip_to;
#
#        AUTHOR:  Ryan Schulze (rs), ryan@dopefish.de
#
#===============================================================================

set -o nounset                              # Treat unset variables as an error

CSVFile="GeoIPCountryWhois.csv"

if [[ -e ${CSVFile} ]]
  then
  echo -en "\n${CSVFile} already exists, download a newer version? [Y/n]: "
  read answer
  if [[ "$(echo $answer|tr "N" "n")" != "n" ]]
  then
    rm ${CSVFile}
  fi
fi

if [[ ! -e ${CSVFile} ]]
then
  tmpfile=$(mktemp)
  wget -O ${tmpfile} http://geolite.maxmind.com/download/geoip/database/GeoIPCountryCSV.zip
  unzip ${tmpfile} &amp;&amp; rm ${tmpfile}
fi

echo -e "\ngenerating geoip.sql"

mysql -e "create database if not exists temp"
mysql -D temp -e "drop table if exists geo_ip"

echo '
CREATE TABLE geo_ip
(
  id           INT UNSIGNED  NOT NULL auto_increment,
  ip_poly      POLYGON       NOT NULL,
  ip_from      INT UNSIGNED  NOT NULL,
  ip_to        INT UNSIGNED  NOT NULL,
  country_code CHAR(2)       NOT NULL,
  country_name CHAR(50)      NOT NULL,
  PRIMARY KEY (id),
  SPATIAL INDEX (ip_poly)
) ENGINE=MyISAM DEFAULT CHARSET=latin1;
LOAD DATA LOCAL INFILE "GeoIPCountryWhois.csv"
INTO TABLE geo_ip
FIELDS
  TERMINATED BY ","
  ENCLOSED BY "\""
LINES
  TERMINATED BY "\n"
(
  @ip_from_string, @ip_to_string,
  @ip_from, @ip_to,
  @country_code, @country_string
)
SET
  id      := NULL,
  ip_from := @ip_from,
  ip_to   := @ip_to,
  ip_poly := GEOMFROMWKB(POLYGON(LINESTRING(
    /* clockwise, 4 points and back to 0 */
    POINT(@ip_from, -1), /* 0, top left */
    POINT(@ip_to,   -1), /* 1, top right */
    POINT(@ip_to,    1), /* 2, bottom right */
    POINT(@ip_from,  1), /* 3, bottom left */
    POINT(@ip_from, -1)  /* 0, back to start */
  ))),
  country_code := @country_code,
  country_name := @country_string
;
' | mysql --local-infile=1 -D temp

mysqldump --opt temp geo_ip &gt; geoip.sql

echo -e "\ndone generating geoip.sql"

#!/bin/bash

#===============================================================================

# FILE: generate_geoip_sql.sh

# USAGE: ./generate_geoip_sql.sh

# DESCRIPTION: Generates a optimized sql dump with geoipdata from maxmind.com

# USAGE: import sql into a database

# REQUIREMENTS: mysql database rights to create databases and LOAD DATA LOCAL INFILE

# EXAMPLE:

# SELECT country_code,country_name

# FROM geo_ip

# WHERE MBRCONTAINS(ip_poly, POINTFROMWKB(POINT(INET_ATON('1.2.3.4'), 0)))

# SELECT country_code,country_name

# FROM geo_ip

# WHERE INET_ATON('1.2.3.4')

# BETWEEN ip_from AND ip_to;

# AUTHOR: Ryan Schulze (rs), ryan@dopefish.de

#===============================================================================

set -o nounset # Treat unset variables as an error

CSVFile="GeoIPCountryWhois.csv"

if [[ -e ${CSVFile} ]]

then

echo -en "\n${CSVFile} already exists, download a newer version? [Y/n]: "

read answer

if [[ "$(echo $answer|tr "N" "n")" != "n" ]]

then

rm ${CSVFile}

if [[ ! -e ${CSVFile} ]]

then

tmpfile=$(mktemp)

wget -O ${tmpfile} http://geolite.maxmind.com/download/geoip/database/GeoIPCountryCSV.zip

unzip ${tmpfile} && rm ${tmpfile}

echo -e "\ngenerating geoip.sql"

mysql -e "create database if not exists temp"

mysql -D temp -e "drop table if exists geo_ip"

echo '

CREATE TABLE geo_ip

(

id INT UNSIGNED NOT NULL auto_increment,

ip_poly POLYGON NOT NULL,

ip_from INT UNSIGNED NOT NULL,

ip_to INT UNSIGNED NOT NULL,

country_code CHAR(2) NOT NULL,

country_name CHAR(50) NOT NULL,

PRIMARY KEY (id),

SPATIAL INDEX (ip_poly)

) ENGINE=MyISAM DEFAULT CHARSET=latin1;

LOAD DATA LOCAL INFILE "GeoIPCountryWhois.csv"

INTO TABLE geo_ip

FIELDS

TERMINATED BY ","

ENCLOSED BY "\""

LINES

TERMINATED BY "\n"

(

@ip_from_string, @ip_to_string,

@ip_from, @ip_to,

@country_code, @country_string

)

SET

id := NULL,

ip_from := @ip_from,

ip_to := @ip_to,

ip_poly := GEOMFROMWKB(POLYGON(LINESTRING(

/* clockwise, 4 points and back to 0 */

POINT(@ip_from, -1), /* 0, top left */

POINT(@ip_to, -1), /* 1, top right */

POINT(@ip_to, 1), /* 2, bottom right */

POINT(@ip_from, 1), /* 3, bottom left */

POINT(@ip_from, -1) /* 0, back to start */

))),

country_code := @country_code,

country_name := @country_string

;

' | mysql --local-infile=1 -D temp

mysqldump --opt temp geo_ip > geoip.sql

echo -e "\ndone generating geoip.sql"