原创：（网格化地图之初）高斯大地座标（经度、纬度）→3°投影带平面直角坐标（X、Y）换算

好像不能上传Excel，暂时把结果用表格显示出来吧

高斯大地座标（经度、纬度）→3°投影带平面直角坐标（X、Y）换算表											数     据      预      处        理
投影点编号	纬度(B)			经度(L)			投影带号	纵座标（X,m）	横座标（Y,m）	子午线收敛角    （γ°）	B弧度值	L弧度值	X0	l=经差    （L－L0）	卯酉圈曲率半径，N	t=tgB	ρ’	η2	座标带编号	5-t2+9η2+4η4	61-58t2+t4
	度	分	秒	度	分	秒
1	35	32	7	107	4	36	35	3936032.1591	35688348.2838	1.206980965	0.62020735	1.868840253	3934047.827	0.036244538	6385467.934	0.7142	57.2958	0.0045	35	4.5301	31.6737

以下是Oracle中创建的计算函数，可以根据经纬度计算出平面直角坐标的X、Y的距离，进而在网格化的地图中得知该点位于哪个方格或者哪个象限！

CREATE OR REPLACE FUNCTION Alpha(N number) RETURN NUMBER
 is
 --η2
 begin
     return (0.006738525415)*power(COS(Radian(N)),2);
 end ;
CREATE OR REPLACE FUNCTION Beta(N number) RETURN NUMBER
 is
 --ρ’
 begin
      return 180/acos(-1);
 end ;
CREATE OR REPLACE FUNCTION Coordinate(Lon number) RETURN NUMBER
 is
 --座标带编号/投影带号
 begin
     return trunc(Lon/3,0)+1;
 end ;
CREATE OR REPLACE FUNCTION Fun_Meridian_C_angle(Lon number,Lat number) RETURN NUMBER
 is
 --子午线收敛角    （γ°）
 begin
     return Beta(Lat)*Lon_Diff(Lon)*SIN(RADIAN(Lat))+SIN(RADIAN(Lat))*POWER(COS(RADIAN(Lat)),2)*(1+3*Alpha(Lat)+2*POWER(Alpha(Lat),2))*POWER(Lon_Diff(Lon),3)+SIN(RADIAN(Lat))*power(COS(RADIAN(Lat)),4)*(2-power(Tange(Lat),2))*power(Lon_Diff(Lon),5);
 end ;
CREATE OR REPLACE FUNCTION Fun_X0(N number) RETURN NUMBER
 is
 --得出范围X0
 begin
 return  6378245*(1-0.006693421623)*(1.0050517739*Radian(N)-0.00506237764*SIN(2*Radian(N))/2+0.00001062451*SIN(4*Radian(N))/4-0.0000002081*SIN(6*Radian(N))/6);
 end ;
CREATE OR REPLACE FUNCTION GETBUSLINEFIRSTP(lineverid number, theupdown number, firstlastpointtag integer) RETURN VARCHAR2
 is
 Result VARCHAR2(4000); --定义变量
 --公交线第一个点
 begin
  if firstlastpointtag = 1 then
  select 'POINT(' || a.longitude_offset || ' ' || a.latitude_offset || ')' into Result from TB_BSINFO_BUSLINE_ROUTE a where a.order_num=(select min(b.order_num) from TB_BSINFO_BUSLINE_ROUTE b where b.line_ver_id=lineverid and b.updown=theupdown) and a.line_ver_id=lineverid and a.updown=theupdown;
  end if;
  if firstlastpointtag = 2 then
  select 'POINT(' || a.longitude_offset || ' ' || a.latitude_offset || ')' into Result from TB_BSINFO_BUSLINE_ROUTE a where a.order_num=(select max(b.order_num) from TB_BSINFO_BUSLINE_ROUTE b where b.line_ver_id=lineverid and b.updown=theupdown) and a.line_ver_id=lineverid and a.updown=theupdown;
  end if;
  return Result;
 end;
CREATE OR REPLACE FUNCTION GETBUSLINEOVERLAPLENGTH(aline_ver_id number,aupdown number,bline_ver_id number,bupdown number) RETURN VARCHAR2
 is
 a_line_ver_id number;
 a_updown number;
 b_line_ver_id number;
 b_updown number;
 overlapbuslinelenght VARCHAR2(4000); --定义变量
 a_buslinelenght VARCHAR2(4000);
 b_buslinelenght VARCHAR2(4000);
 begin

    a_line_ver_id := aline_ver_id;
    a_updown := aupdown;
    b_line_ver_id := bline_ver_id;
    b_updown := bupdown;
    if aline_ver_id > bline_ver_id then
      a_line_ver_id := bline_ver_id;
      a_updown := bupdown;
      b_line_ver_id := aline_ver_id;
      b_updown := aupdown;
    end if;

    select
      sdo_geom.sdo_length(a.shape, 0.1, 'unit=M') ,
      sdo_geom.sdo_length(c.shape, 0.1, 'unit=M'),
      sdo_geom.sdo_length(SDO_GEOM.sdo_intersection(a.shape,SDO_GEOM.SDO_BUFFER(c.shape, 15, 0.2, 'unit= M'),0.001), 0.1, 'unit=M')
      into a_buslinelenght,b_buslinelenght,overlapbuslinelenght
      from
      (select b.shape as shape from TB_BSINFO_BUSLINE_SHAPE b
      where b.line_ver_id=b_line_ver_id and b.updown=b_updown) c,TB_BSINFO_BUSLINE_SHAPE a
      where a.line_ver_id=a_line_ver_id and a.updown=a_updown;
    if overlapbuslinelenght > 15 then
      overlapbuslinelenght := overlapbuslinelenght - 15;
      if aline_ver_id > bline_ver_id then
         overlapbuslinelenght := to_char(b_buslinelenght) || '~' || to_char(overlapbuslinelenght) || '~' || to_char(a_buslinelenght);
      else
        overlapbuslinelenght := to_char(a_buslinelenght) || '~' || to_char(overlapbuslinelenght) || '~' || to_char(b_buslinelenght);
      end if;
      return overlapbuslinelenght;
    else
      select
      sdo_geom.sdo_length(a.shape, 0.1, 'unit=M') ,
      sdo_geom.sdo_length(c.shape, 0.1, 'unit=M'),
      sdo_geom.sdo_length(SDO_GEOM.sdo_intersection(a.shape,SDO_GEOM.SDO_BUFFER(c.shape, 15, 0.2, 'unit= M'),0.001), 0.1, 'unit=M')
      into a_buslinelenght,b_buslinelenght,overlapbuslinelenght
      from
      (select b.shape as shape from TB_BSINFO_BUSLINE_SHAPE b
      where b.line_ver_id=a_line_ver_id and b.updown=a_updown) c,TB_BSINFO_BUSLINE_SHAPE a
      where a.line_ver_id=b_line_ver_id and a.updown=b_updown;
      overlapbuslinelenght := overlapbuslinelenght - 15;
      if aline_ver_id > bline_ver_id then
         overlapbuslinelenght := to_char(b_buslinelenght) || '~' || to_char(overlapbuslinelenght) || '~' || to_char(a_buslinelenght);
      else
        overlapbuslinelenght := to_char(a_buslinelenght) || '~' || to_char(overlapbuslinelenght) || '~' || to_char(b_buslinelenght);
      end if;
      return overlapbuslinelenght;
    end if;

 end;
CREATE OR REPLACE FUNCTION GetDistance(lng1 number,
                                        lat1 number,
                                        lng2 number,
                                        lat2 number) RETURN NUMBER is
   earth_padius number := 6378.137;
   radLat1      number := Radian(lat1);
   radLat2      number := Radian(lat2);
   a            number := radLat1 - radLat2;
   b            number := Radian(lng1) - Radian(lng2);
   s            number := 0;
 begin
   s := 2 *
        Asin(Sqrt(power(sin(a / 2), 2) +
                  cos(radLat1) * cos(radLat2) * power(sin(b / 2), 2)));
   s := s * earth_padius;
   s := Round(s * 10000) / 10000;
   return s*1000;

 end;
create or replace function getGrid_Code(Lon in number, lat in number) return varchar2 is
  Load_X       number(10,2);
  Load_Y       number(10,2);
  Addr_X       number(10,2);
  Addr_Y       number(10,2);
  NumOne       number;
  NumY         number;
  NumFour      number;
  NumX         number;
  NumYcode     number;
  NumXcode     number;
 -- Result       varchar2;
begin
  --人民广场上海国际饭店121.46711830,31.23546313横纵坐标
  select X_Coordinate(121.46711830,31.23546313),Y_Coordinate(121.46711830,31.23546313) into Load_X,Load_Y from dual;
  select X_Coordinate(lon,lat),Y_Coordinate(lon,lat) into Addr_X,Addr_Y from dual;
  if (Addr_Y - Load_Y) > 0 THEN
    NumOne := 1;
    else
      NumOne := 3;
    end if;
   if (Addr_X - Load_X) > 0 then
     NumFour := 4;
     ELSE
       NumFour := 2;
    END IF;
    select round(abs((Addr_Y - Load_Y)/1000),0),round(abs((Addr_X - Load_X)/1000),0) into NumY, NumX from dual;
    if NumY < 10 then
       NumYcode :=NumOne||'0'||NumY;
    ELSE
      NumYcode:=NumOne||NumY;
    end if;
    if NumX < 10 then
       NumXcode := NumFour||'0'||NumX;
    ELSE
      NumXcode:= NumFour||NumX;
    end if;
  return(NumYcode||NumXcode);
end getGrid_Code;
CREATE OR REPLACE FUNCTION GETNEXTSTATION(lineverid number, theupdown number, currentordernum number) RETURN VARCHAR2
 is
 nextstationname VARCHAR2(4000); --定义变量
 nextordernum     integer;
 begin
  nextordernum := currentordernum + 1;
  select a.name into nextstationname from TB_BSINFO_BUSLINE_STATION t,TB_BSINFO_STATION_DICT a
         where t.line_ver_id=lineverid and t.updown=theupdown and t.order_num=nextordernum and t.station_id=a.station_id;
  return nextstationname;
 end;
CREATE OR REPLACE FUNCTION GETSTATIONBUSLINE(stationid number) RETURN VARCHAR2
 is
 buslines VARCHAR2(4000); --定义变量
 begin
   select WMSYS.Wm_Concat(b.name || '(版本：' || t.line_ver_id || '-上下行:' || t.updown || ')') into buslines from TB_BSINFO_BUSLINE_STATION t, TB_BSINFO_BUSLINE_VER a,TB_BSINFO_BUSLINE_DICT b
         where t.station_id = stationid and t.line_ver_id=a.id and a.line_id=b.line_id;

  return buslines;
 end;
CREATE OR REPLACE FUNCTION GETSTATIONBUSLINENUM(stationid number) RETURN VARCHAR2
 is
 buslinenum VARCHAR2(4000); --定义变量
 begin
  select count(t.station_id) into buslinenum from TB_BSINFO_BUSLINE_STATION t
         where t.station_id = stationid;
  return buslinenum;
 end;
create or replace function Get_Gps_time(gps_time in number) return varchar2 is
  Result varchar2(50);
begin
  select SUBSTR(TO_CHAR(gps_time,'FM000000'),1,2)||':'||SUBSTR(TO_CHAR(gps_time,'FM000000'),3,2)||':'||SUBSTR(TO_CHAR(gps_time,'FM000000'),-2) into Result from dual;
  return(Result);
end Get_Gps_time;
CREATE OR REPLACE FUNCTION Lon_Diff(Lon number) RETURN NUMBER
 is
 --l=经差 （L－L0）
 begin
     return Radian(Lon)-3*Coordinate(Lon)*acos(-1)/180;
 end ;
CREATE OR REPLACE FUNCTION Radian(N number) RETURN NUMBER
 is
 --得出弧度
 begin
 return  N* acos(-1)/180.0;
 end ;
CREATE OR REPLACE FUNCTION ROCIPV(N number) RETURN NUMBER
 is
 --卯酉圈曲率半径，N
 begin
      return 6378245/SQRT(1-0.006693421623*power((SIN(Radian(N))),2));
 end ;
CREATE OR REPLACE FUNCTION Tange(N number) RETURN NUMBER
 is
 --t=tgB
 begin
      return TAN(Radian(N));
 end ;
CREATE OR REPLACE FUNCTION U_QSS(N number) RETURN NUMBER
 is
 --5-t2+9η2+4η4
 begin
     return 5-power(Tange(N),2)+9*Alpha(N)+4*power(Alpha(N),2);
 end ;
CREATE OR REPLACE FUNCTION V_QQ(N number) RETURN NUMBER
 is
 --61-58t2+t4
 begin
     return 61-58*power(Tange(N),2)+power(Tange(N),4);
 end ;
CREATE OR REPLACE FUNCTION X_Coordinate(Lon number,lat number) RETURN NUMBER
 is
 --纵座标（X,m）
 begin
     return Fun_X0(lat)+POWER(Lon_Diff(Lon),2)*ROCIPV(lat)*SIN(radian(lat))*COS(radian(lat))/2+POWER(Lon_Diff(Lon),4)*ROCIPV(lat)*SIN(radian(lat))*POWER(COS(radian(lat)),3)*(U_QSS(lat))/24+POWER(Lon_Diff(Lon),6)*ROCIPV(lat)*SIN(radian(lat))*POWER(COS(radian(lat)),5)*V_QQ(lat)/720;
 end ;
CREATE OR REPLACE FUNCTION Y_Coordinate(Lon number,lat number) RETURN NUMBER
 is
 --横座标（Y,m）
 begin
     return Coordinate(Lon)*1000000+500000+Lon_Diff(Lon)*ROCIPV(lat)*COS(RADIAN(lat))*(1+POWER(Lon_Diff(Lon)*COS(RADIAN(lat)),2)*(1-POWER(Tange(lat),2)+Alpha(lat))/6+POWER(Lon_Diff(Lon)*COS(RADIAN(lat)),4)*(5-18*POWER(Tange(lat),2)+POWER(Tange(lat),4)+14*Alpha(lat)-58*Alpha(lat)*power(Tange(lat),2))/120);
 end ;