find minimum bounding coordinate

find_minimum_bounding_coordinates.sql

-- find minimum bounding coordinate
-- ref:
-- [1] http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
-- check the ref[1]
--
-- find minimum bounding coordinate,
-- if lon_min or lon_max outside [-π, π],
-- then the lon_min or lat_min will translate to lon_min + 2π or lon_max - 2π,
-- when check condition for this situation, need to check lon <= lon_min or lon >= lon_max


delimiter $$
drop function if exists geo_bounding_coordinates_contains$$
create function geo_bounding_coordinates_contains (bounding_min POINT, bounding_max POINT, pt point)
       returns boolean
       deterministic
       comment 'POINT(longitude, latitude)'
       begin

	IF x(bounding_min) < x(bounding_max) THEN -- normal situation
	   return (x(pt) >= x(bounding_min) AND x(pt) <= x(bounding_max)) AND (y(pt) >= y(bounding_min) AND y(pt) <= y(bounding_max));
	ELSE
	   return (x(pt) <= x(bounding_max) OR x(pt) >= x(bounding_min)) AND (y(pt) >= y(bounding_min) AND y(pt) <= y(bounding_max));
	END IF;
       end$$

drop procedure if exists find_minimum_bounding_coordinates_degree$$
create procedure find_minimum_bounding_coordinates_degree (center POINT, distance DOUBLE, OUT coord_min POINT, OUT coord_max POINT)
       deterministic
       begin
	call find_minimum_bounding_coordinates_degree_R(center, distance, coord_min, coord_max, 6370986);
       end$$
       
drop procedure if exists find_minimum_bounding_coordinates_radian$$
create procedure find_minimum_bounding_coordinates_radian (center POINT, distance DOUBLE, OUT coord_min POINT, OUT coord_max POINT)
       deterministic
       begin
       	call find_minimum_bounding_coordinates_radian_R(center, distance, coord_min, coord_max, 6370986);
       end$$

drop procedure if exists find_minimum_bounding_coordinates_degree_R$$
create procedure find_minimum_bounding_coordinates_degree_R (center POINT, distance DOUBLE, OUT coord_min POINT, OUT coord_max POINT, R DOUBLE)
       deterministic
       begin
       declare center_rad POINT;
       declare coord_min_rad POINT;
       declare coord_max_rad POINT;

       set center_rad = POINT(RADIANS(x(center)), RADIANS(y(center)));
       
       call find_minimum_bounding_coordinates_radian_R(center_rad, distance, coord_min_rad, coord_max_rad, R);

       set coord_min = POINT(DEGREES(x(coord_min_rad)) , DEGREES(y(coord_min_rad)));
       set coord_max = POINT(DEGREES(x(coord_max_rad)) , DEGREES(y(coord_max_rad)));
       end$$
       
drop procedure if exists find_minimum_bounding_coordinates_radian_R$$
create procedure find_minimum_bounding_coordinates_radian_R (center POINT, distance DOUBLE, OUT coord_min POINT, OUT coord_max POINT,  R DOUBLE)
       deterministic
       begin
	
	declare lat_min DOUBLE;
	declare lat_max DOUBLE;
	declare lon_min DOUBLE;
	declare lon_max DOUBLE;

	declare `α` DOUBLE; -- angular radius
	declare lat DOUBLE;
	declare lon DOUBLE;

	declare lat_t DOUBLE;
	declare `Δlon` DOUBLE;

	declare `π` DOUBLE default pi();
	declare `π_2` DOUBLE;
	set `π_2` = `π` / 2;

	set `α` = distance / R;

	set lat = y(center);
	set lon = x(center);
	set lat_min = lat - `α`;
	set lat_max = lat + `α`;

	
	IF lat_max > `π_2` OR lat_min < -`π_2` THEN
	   set coord_min = POINT(-`π`, GREATEST(lat_min, -`π_2`));
	   set coord_max = POINT(`π`, LEAST(lat_max, `π_2`));
	ELSE
	   set lat_t = asin(sin(lat)/cos(`α`));
           -- set `Δlon` = asin(sin(`α`) / cos(lat));
	   set `Δlon` = acos((cos(`α`) - sin(lat_t)*sin(lat)) / (cos(lat_t) * cos(lat)));
	   set lon_min = lon - `Δlon`;
	   set lon_max = lon + `Δlon`;

	   if lon_min < -`π` THEN
	      set lon_min = lon_min + 2 * `π`;
	   END IF;
	   if lon_max > `π` THEN
	      set lon_max = lon_max - 2 * `π`;
	   END IF;
	   
	   set coord_min = POINT(lon_min, lat_min);
	   set coord_max = POINT(lon_max, lat_max);
	END IF;

	-- set coord_min = POINT(DEGREES(x(coord_min)), DEGREES(y(coord_min)));
	-- set coord_max = POINT(DEGREES(x(coord_max)), DEGREES(y(coord_max)));
	-- select astext(coord_min), astext(coord_max), `Δlon`;
	-- select lat_min, lat_max, lat_t, lon_min, lon_max,`Δlon`, `α`, sin(`α`), cos(lat), sin(`α`) / cos(lat);

       end$$

delimiter ;