说明

之前一个帖子给出了正反算的C语言代码
链接如下：
link

我把它移植到了matlab中
支持向量输入

高斯投影正算

高斯投影正算matlab代码

function [X, Y] = GaussProjCal(Lon, Lat)
%%  高斯投影正算： 经纬度 算 直角坐标
% Longitude  Latitude 均为角度制
Lon = deg2rad(Lon);
Lat = deg2rad(Lat);

ZoneWide = deg2rad(6); %带宽为6度
a=6378140.0; f=1/298.257;  %80年西安坐标系参数

ProjNo = floor(Lon./ZoneWide);

lon0 = ProjNo.*ZoneWide+ZoneWide./2;

e2 = 2*f-f*f;
ee = e2*(1.0-e2);
NN = a./sqrt(1.0-e2.*sin(Lat).*sin(Lat));
T = tan(Lat).*tan(Lat);
C = ee.*cos(Lat).*cos(Lat);
A = (Lon-lon0).*cos(Lat);
M = a.*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256).*Lat-(3*e2/8+3*e2*e2/32+45*e2*e2*e2/1024).*sin(2.*Lat)...
  +(15*e2*e2/256+45*e2*e2*e2/1024).*sin(4.*Lat)-(35*e2*e2*e2/3072).*sin(6.*Lat));

xval = NN.*(A+(1-T+C).*A.*A.*A./6+(5-18.*T+T.*T+72.*C-58.*ee).*A.*A.*A.*A.*A/120);
yval = M+NN.*tan(Lat).*(A.*A./2+(5-T+9.*C+4.*C.*C).*A.*A.*A.*A./24 ...
       +(61-58.*T+T.*T+600.*C-330.*ee).*A.*A.*A.*A.*A.*A./720);

X0 = 1000000.*(ProjNo+1)+500000; 
Y0 = 0; 

xval = xval+X0; yval = yval+Y0; 
X = xval;
Y = yval;

高斯投影反算

高斯投影反算matlab代码。

function [Lon, Lat] = GaussProjInvCal(X,Y)
%%   高斯投影反算：  直角坐标 反算 经纬度
a=6378140.0; f=1/298.257; %80年西安坐标系参数
ZoneWide = deg2rad(6); %6度带宽
ProjNo = floor(X./1000000);

longitude0 = (ProjNo-1).*ZoneWide+ZoneWide./2;

X0 = ProjNo.*1000000+500000;
Y0=0;
xval = X-X0; yval = Y-Y0;

e2 = 2.*f-f.*f;
e1 = (1.0-sqrt(1-e2))./(1.0+sqrt(1-e2));
ee = e2./(1-e2);
M = yval;
u = M./(a.*(1-e2./4-3.*e2.*e2/64-5.*e2.*e2.*e2/256));
fai = u+(3.*e1/2-27.*e1.*e1.*e1./32).*sin(2.*u)+(21.*e1.*e1/16-55.*e1.*e1.*e1.*e1./32)*sin(4.*u)...
        +(151.*e1.*e1.*e1./96).*sin(6.*u)+(1097.*e1.*e1.*e1.*e1./512).*sin(8.*u);
C = ee.*cos(fai).*cos(fai);
T = tan(fai).*tan(fai);
NN = a./sqrt(1.0-e2.*sin(fai).*sin(fai));
R = a.*(1-e2)./sqrt((1-e2.*sin(fai).*sin(fai)).*(1-e2.*sin(fai).*sin(fai)).*(1-e2.*sin(fai).*sin(fai)));
D = xval./NN;

%计算经度（Longitude） 纬度（Latitude）
Lon = longitude0+(D-(1+2.*T+C).*D.*D.*D./6+(5-2.*C+28.*T-3.*C.*C+8.*ee+24.*T.*T).*D.*D.*D.*D.*D./120)./cos(fai);
Lon = rad2deg(Lon);
Lat = fai -(NN.*tan(fai)./R).*(D.*D./2-(5+3.*T+10.*C-4.*C.*C-9.*ee).*D.*D.*D.*D./24+(61+90.*T+298.*C+45.*T.*T-256.*ee-3.*C.*C).*D.*D.*D.*D.*D.*D./720);
Lat = rad2deg(Lat);