MySQL 转换GPS84为百度BD09的函数

-- 使用示例
SELECT ConvertPolygonToBD("POLYGON((109.04876 31.325246,109.048966 31.326721,109.049111 31.328609,109.048909 31.33047,109.04807 31.33334))");

-- 创建转换函数
DELIMITER $
DROP FUNCTION ConvertPolygonToBD$
CREATE FUNCTION ConvertPolygonToBD(GpsPolygon LONGTEXT)
RETURNS LONGTEXT
BEGIN
 -- 转换 GPS 84坐标到 BD09格式
 SET @pi = 3.14159265358979324;
 SET @x_pi = 3.14159265358979324 * 3000.0 / 180.0;
 SET @a = 6378245.0;
 SET @ee = 0.00669342162296594323;
 
 IF GpsPolygon IS NULL THEN RETURN ""; END IF;
 
 SET @R = "POLYGON((";
 SET @Y = REPLACE(GpsPolygon, @R, "");
 SET @Y = REPLACE(@Y, "))", "");
 
 REPEAT
   SET @L = SUBSTRING_INDEX(@Y, ",", 1);
   SET @Y = SUBSTRING(@Y, POSITION("," IN @Y) + 1);
   
   SET @LON = SUBSTRING_INDEX(@L, " ", 1);
   SET @LAT = SUBSTRING_INDEX(@L, " ", -1);
   
   -- 转换 WGS 84 到 GCJ 02（国内坐标系）
   SET @X1 = @LON - 105.0; SET @Y1 = @LAT - 35.0;
   
   SET @DLAT = -100.0 + 2.0 * @X1 + 3.0 * @Y1 + 0.2 * @Y1 * @Y1 + 0.1 * @X1 * @Y1 + 0.2 * SQRT(ABS(@X1));
   SET @DLAT = @DLAT + (20.0 * SIN(6.0 * @X1 * @PI) + 20.0 * SIN(2.0 * @x1 * @PI)) * 2.0 / 3.0;
   SET @DLAT = @DLAT + (20.0 * SIN(@y1 * @PI) + 40.0 * SIN(@y1 / 3.0 * @PI)) * 2.0 / 3.0;
   SET @DLAT = @DLAT + (160.0 * SIN(@y1 / 12.0 * @PI) + 320 * SIN(@y1 * @PI / 30.0)) * 2.0 / 3.0;
   
   SET @DLON = 300.0 + @x1 + 2.0 * @y1 + 0.1 * @x1 * @x1 + 0.1 * @x1 * @y1 + 0.1 * SQRT(ABS(@x1));
   SET @DLON = @DLON + (20.0 * SIN(6.0 * @x1 * @PI) + 20.0 * SIN(2.0 * @x1 * @PI)) * 2.0 / 3.0;
   SET @DLON = @DLON + (20.0 * SIN(@x1 * @PI) + 40.0 * SIN(@x1 / 3.0 * @PI)) * 2.0 / 3.0;
   SET @DLON = @DLON + (150.0 * SIN(@x1 / 12.0 * @PI) + 300.0 * SIN(@x1 / 30.0 * @PI)) * 2.0 / 3.0;
   
   SET @radLat = @LAT / 180.0 * @PI;
   SET @MAGIC = SIN(@radLat);
   SET @MAGIC = 1 - @ee * @MAGIC * @MAGIC;
   SET @sqrtMagic = SQRT(@magic);
   
   SET @DLAT = (@dLat * 180.0) / ((@a * (1 - @ee)) / (@magic * @sqrtMagic) * @PI);
   SET @dLon = (@dLon * 180.0) / (@a / @sqrtMagic * COS(@radLat) * @PI);
   -- 转换 GCJ 02 到 BD 09
   SET @Z = SQRT(@DLON * @DLON + @DLAT * @DLAT) + 0.00002 * SIN(@DLAT * @x_pi);
   SET @theta = ATAN2(@DLAT, @DLON) + 0.000003 * COS(@DLON * @x_pi);  
   
   SET @DLON = CAST(@z * COS(@theta) + 0.0065 AS DECIMAL(11,7));
   SET @DLAT = CAST(@z * SIN(@theta) + 0.006 AS DECIMAL(11,7));
   
   SET @L = CONCAT(@LON + @DLON, " ", @LAT + @DLAT);
   SET @R = CONCAT(@R, @L, ",");
 UNTIL (POSITION("," IN @Y) = 0)
 END REPEAT;
 
 -- 处理最后一个数据
   SET @LON = SUBSTRING_INDEX(@Y, " ", 1);          SET @LAT = SUBSTRING_INDEX(@Y, " ", -1);
   SET @X1 = @LON - 105.0;                        SET @Y1 = @LAT - 35.0;
   
   SET @DLAT = -100.0 + 2.0 * @X1 + 3.0 * @Y1 + 0.2 * @Y1 * @Y1 + 0.1 * @X1 * @Y1 + 0.2 * SQRT(ABS(@X1));
   SET @DLAT = @DLAT + (20.0 * SIN(6.0 * @X1 * @PI) + 20.0 * SIN(2.0 * @x1 * @PI)) * 2.0 / 3.0;
   SET @DLAT = @DLAT + (20.0 * SIN(@y1 * @PI) + 40.0 * SIN(@y1 / 3.0 * @PI)) * 2.0 / 3.0;
   SET @DLAT = @DLAT + (160.0 * SIN(@y1 / 12.0 * @PI) + 320 * SIN(@y1 * @PI / 30.0)) * 2.0 / 3.0;
   
   SET @DLON = 300.0 + @x1 + 2.0 * @y1 + 0.1 * @x1 * @x1 + 0.1 * @x1 * @y1 + 0.1 * SQRT(ABS(@x1));
   SET @DLON = @DLON + (20.0 * SIN(6.0 * @x1 * @PI) + 20.0 * SIN(2.0 * @x1 * @PI)) * 2.0 / 3.0;
   SET @DLON = @DLON + (20.0 * SIN(@x1 * @PI) + 40.0 * SIN(@x1 / 3.0 * @PI)) * 2.0 / 3.0;
   SET @DLON = @DLON + (150.0 * SIN(@x1 / 12.0 * @PI) + 300.0 * SIN(@x1 / 30.0 * @PI)) * 2.0 / 3.0;
   
   SET @radLat = @LAT / 180.0 * @PI;
   SET @MAGIC = SIN(@radLat);
   SET @MAGIC = 1 - @ee * @MAGIC * @MAGIC;
   SET @sqrtMagic = SQRT(@magic);
   
   SET @DLAT = (@dLat * 180.0) / ((@a * (1 - @ee)) / (@magic * @sqrtMagic) * @PI);
   SET @dLon = (@dLon * 180.0) / (@a / @sqrtMagic * COS(@radLat) * @PI);
   SET @Z = SQRT(@DLON * @DLON + @DLAT * @DLAT) + 0.00002 * SIN(@DLAT * @x_pi);
   SET @theta = ATAN2(@DLAT, @DLON) + 0.000003 * COS(@DLON * @x_pi);  
   
   SET @DLON = CAST(@z * COS(@theta) + 0.0065 AS DECIMAL(11,7));
   SET @DLAT = CAST(@z * SIN(@theta) + 0.006 AS DECIMAL(11,7));
   SET @L = CONCAT(@LON + @DLON, " ", @LAT + @DLAT);
 
 RETURN CONCAT(@R, @L, "))");
END$
DELIMITER ;

/*
-- 使用示例
UPDATE SYS_REGION SET PGN_BD = POLYGONFROMTEXT(ConvertPolygonToBD(ASTEXT(pgn)))
WHERE city = "奉节县"
*/