
-- ---------------------------------------------------------------
-- Authors: Matteo Quintiliani and Stefano Pintore
--          I.N.G.V. - Italy
--          Istituto Nazionale di Geofisica e Vulcanologia
-- ---------------------------------------------------------------

-- ---------------------------------------------------------------


-- http://en.wikipedia.org/wiki/Point_in_polygon
-- http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
-- http://www.sql-statements.com/point-in-polygon.html
-- int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
-- {
--   int i, j, c = 0;
--   for (i = 0, j = nvert-1; i < nvert; j = i++) {
--     if ( ((verty[i]>testy) != (verty[j]>testy)) &&
-- 	 (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
--        c = !c;
--   }
--   return c;
-- }


DELIMITER $$
DROP FUNCTION IF EXISTS `fn_pnpoly`$$
CREATE  DEFINER = CURRENT_USER  FUNCTION `fn_pnpoly`(
in_g    GEOMETRY,
in_X    DOUBLE,
in_Y    DOUBLE
)
RETURNS INT DETERMINISTIC
BEGIN
    DECLARE l_linestring          LineString DEFAULT NULL;
    DECLARE l_num_points          BIGINT DEFAULT NULL;
    DECLARE l_Xj                  DOUBLE DEFAULT 0;
    DECLARE l_Xi                  DOUBLE DEFAULT 0;
    DECLARE l_Yj                  DOUBLE DEFAULT 0;
    DECLARE l_Yi                  DOUBLE DEFAULT 0;
    DECLARE l_in                  INT DEFAULT -1;

    DECLARE i                     INT DEFAULT 0;

    SET l_linestring=ExteriorRing(in_g);
    SET l_num_points=NumPoints(l_linestring);

    SET l_Xj = X(PointN(l_linestring, l_num_points));
    SET l_Yj = Y(PointN(l_linestring, l_num_points));
    SET i = 1;
    WHILE (i <= l_num_points) DO
	SET l_Xi = X(PointN(l_linestring, i));
	SET l_Yi = Y(PointN(l_linestring, i));
	IF( ((l_Yi>in_Y) <> (l_Yj>in_Y)) AND
	    ( in_X <  ( (l_Xj-l_Xi) * ((in_Y-l_Yi) / (l_Yj-l_Yi)) + l_Xi) ) ) THEN
		SET l_in = l_in * (-1);
	END IF;
	SET l_Xj = l_Xi;
	SET l_Yj = l_Yi;
	SET i = i + 1;
    END WHILE;

    RETURN l_in;
END $$
DELIMITER ;

-- ---------------------------------------------------------------

DELIMITER $$
DROP FUNCTION IF EXISTS `fn_find_region_name`$$
CREATE  DEFINER = CURRENT_USER  FUNCTION `fn_find_region_name`(
in_lat                   DOUBLE    /*latitude (North=positive)*/,
in_lon                   DOUBLE    /*longitude (East=positive)*/
) RETURNS VARCHAR(1024) NOT DETERMINISTIC READS SQL DATA
 BEGIN
    DECLARE l_region    VARCHAR(1024);
    DECLARE l_sp        VARCHAR(1024);
    DECLARE l_p         GEOMETRY;

    SET l_region = NULL;
    -- WARNING: bug or problem setting l_p directly using the following instruction
    --          SET l_p = GeomFromText( CONCAT('POINT', '(',CAST(in_lat AS CHAR) ,' ',  CAST(in_lon AS CHAR), ')' ) );
    --  Possible solution: set an intermediate string variable l_sp and successively pass it to GeomFromText()
    --  Calling Contains(g, l_p) is the fastest way I found.
    SET l_sp = CONCAT('POINT', '(',CAST(in_lat AS CHAR) ,' ',  CAST(in_lon AS CHAR), ')' );
    SET l_p = GeomFromText(l_sp);

    SELECT
    region
    /*
    -- For debugging
    GROUP_CONCAT(
	region, '{', kind, '}', IF(fn_pnpoly(g, in_lat, in_lon)>0, '', '[?]')
	ORDER BY fn_pnpoly(g, in_lat, in_lon)  DESC, kind, Area(g)
	SEPARATOR '; '
    )
    */
    INTO l_region
    FROM (
	SELECT *
	FROM ew_region_geom
	WHERE Contains(g, l_p)
    ) tmp_table1
    ORDER BY fn_pnpoly(g, in_lat, in_lon)  DESC, kind, Area(g)
    LIMIT 1
    ;

    if l_region IS NULL then 
	SET l_region = 'Region not found!'; 
    END IF;

    RETURN l_region;
END $$
DELIMITER ;

-- ---------------------------------------------------------------

-- DROP TABLE IF EXISTS `ew_region_kind`;
-- DROP TABLE IF EXISTS `ew_region_geom`;

CREATE TABLE IF NOT EXISTS `ew_region_kind` (
    `kind`               INT               NOT NULL DEFAULT 0 COMMENT 'identifier of kind of the regions',
    `name`               VARCHAR(255)      NOT NULL COMMENT 'Name of the kind of the regions',
    PRIMARY KEY (`kind`)
)
ENGINE = MYISAM COMMENT 'Kind of the regions';

CREATE TABLE IF NOT EXISTS `ew_region_geom` (
    `kind`               INT               NOT NULL DEFAULT 0 COMMENT 'identifier of kind of the region',
    `g`                  GEOMETRY          NOT NULL COMMENT 'Region geometry',
    `region`             VARCHAR(255)      NOT NULL COMMENT 'Region name',
    INDEX (`kind`),
    FOREIGN KEY (`kind`) REFERENCES ew_region_kind(`kind`) ON DELETE CASCADE ON UPDATE CASCADE
)
ENGINE = MYISAM COMMENT 'Region geometry and names';

-- ---------------------------------------------------------------
-- Counterclockwise for polygon definition.
-- Suggestion: first and last point be the same.
-- ---------------------------------------------------------------

-- SET @g = 'POLYGON((44.0 7.0, 47.0 19.0, 47.0 19.0, 47.0 7.0 ,44.0 7.0))';
-- INSERT INTO ew_region_geom VALUES (PolygonFromText(@g),'Nord Italia');
-- 
-- SET @g = 'POLYGON((41.5 7.0, 41.5 19, 44.0 19.0, 44.0 7.0,  41.5 7.0))';
-- INSERT INTO ew_region_geom VALUES (PolygonFromText(@g),'Centro Italia');
-- 
-- SET @g = 'POLYGON((36.0 7.0, 36.0 19.0, 41.5 19.0, 41.5 7.0,  36.0 7.0))';
-- INSERT INTO ew_region_geom VALUES (PolygonFromText(@g),'Sud Italia');

-- SELECT region FROM ew_region_geom WHERE Contains(g, GeomFromText('POINT(9.0 37.0)'));

-- ---------------------------------------------------------------

