用LISP编写一个“高斯正算”（将经纬度转成平面坐标）函数

zark

;http://surveying.wb.psu.edu/sur351/georef/Ellip3.htm
;将球面坐标转换成平面坐标
(setq a 6378137.0  ;长半轴
      f (/ 1.0 298.257223563) ;扁率
      BB 25.0   ;纬度
      LL 118.0   ;经度
      L0 117.0   ;中央子午线
      H  0.0   ;投影高
      )
;测试 (tt A f BB LL L0 H)
(defun tt (A f BB LL L0 H / B DL E12 E22 F L LB N R2 TE X Y)
  (setq B (/ (* BB PI) 180.0))
  (setq L (/ (* LL PI) 180.0))
  (setq e12 (- (* f 2.0) (* f f)));e1平方
  (setq e22 (/ e12 (- 1.0 e12)));e2平方
  (setq N (+ H (/ a (expt (- 1.0 (* e12 (sin B) (sin B))) 0.5))));卯酉圈半径+投影高
  (setq Lb (* N (cos B) B));从赤道到子午线孤长
  (setq Te (/ (sin B) (cos B)))
  (setq R2 (* e22 (cos B) (cos B)))
  (setq Dl (/ (* (- LL L0) Pi) 180))
  (setq X (+
     Lb
     (* 0.5 Te N (expt (* (cos B) Dl) 2))
     (* (/ 1.0 24) Te N (expt (* (cos B) Dl) 4) (- 5.0 (* Te Te) (* -9.0 R2) (* -4.0 R2 R2)))
     (* (/ 1.0 720) Te N (expt (* (cos B) Dl) 6) (- 61.0 (* 58.0 Te Te) (* -1.0 (expt Te 4)) (* -270.0 R2) (* 330.0 Te Te R2)))
     )
)
  (setq Y (+
     (* N (cos B) Dl)
     (* (/ 1.0 6) N (expt (* (cos B) Dl) 3) (- 1.0 (* Te Te) (* -1.0 R2)))
     (* (/ 1.0 120) N (expt (* (cos B) Dl) 5) (- 5.0 (* 18.0 Te Te) (* -1.0 (expt Te 4)) (* -14.0 R2) (* 58.0 Te Te R2)))
     )
)
  (List (rtos X 2 3) (rtos Y 2 3))
  )

程序执行结果
笑脸转换结果
东坐标都不差，就北坐标差这么多，不知道问题出现在哪？
请高手指点一下。谢谢


该贴已经同步到 zark的微博

004

革天明 发表于 2012-7-6 08:09
笑脸是如何加上去的，不是DCL吧，VB？

笑脸----coordMG是一款应用广泛的测量坐标转换软件，和DCL无关。

satan421

感谢分享，回复学习。

zark

自己先顶一下

VBALISPER

建议不用这个软件对比.用苏某某的一个EXCEL版软件对比

Gu_xl

对比一下！
http://apptools.mjtd.com/apptools.php?mod=view&apptoolsid=36

微博评论

1000[粉] = [8] 元 [羞羞] ①78595⑧333

来自 昔虚宦疙陀 的新浪微博

zark

Gu_xl 发表于 2012-7-4 22:53
对比一下！
http://apptools.mjtd.com/apptools.php?mod=view&apptoolsid=36

版主，解压出问题了

zark

高斯坐标正反matlab算程序：
function [x,y]=zgauss(B,L,B0)
%B是弧度为单位
%L B0都是一度为单位
e2=0.006693421622966;  %第一偏心率
ee2=0.006738525414683;  %第二偏心率
a=6378245;                 %长半轴
X=111134.8611*B0-16036.4803*sin(2*B)+16.8281*sin(4*B)-0.022*sin(6*B); %子午弧长
p=3600*180/pi;  %ρ"
W=sqrt(1-e2*(sin(B)^2)); %W
N=a/W;                    %卯酉圈半径
t=tan(B);
yita=ee2*cos(B)*cos(B); %η的平方
l=L-117;                %L
ll=l*3600;
x=X+(N*sin(B)*cos(B)*ll^2)/(2*p^2)+N*sin(B)*(cos(B)^3)*(5-t^2-9*yita)*(ll^4)/(24*p^4);%x坐标
y=N*cos(B)*ll/p+N*(cos(B)^3)*(1-t^2+yita)*ll^3/(6*p^3)+N*(cos(B)^5)*(5-18*t^2+t^4)*ll^5/(120*p^5);%坐标
function [BB,l]=fgauss(x,y)
%在克拉索夫椭球上
e2=0.006693421622966;  %第一偏心率
ee2=0.006738525414683;  %第二偏心率
a=6378245;                 %长半轴
a0=111134.8611;
Bf(1)=x/a0;
Bf(1)=Bf(1)*pi/180;  %把度换算成弧度
%以下用迭代法解出  Bf
for i=1:5
Bf(i+1)=(x-(-16036.4803*sin(2*Bf(i))+16.8281*sin(4*Bf(i))-0.022*sin(6*Bf(i))))/a0*pi/180;
end
%代入公式计算
BBf=Bf(5)；
tf=tan(BBf);
yitaf2=ee2*(cos(BBf)^2);
Wf=sqrt(1-e2*(sin(BBf)^2));
Mf=a*(1-e2)/Wf;
Nf=a/Wf;
BB=BBf-tf*y^2/(2*Mf*Nf)+tf^3*(5+3*tf^2+yitaf2-9*yitaf2*tf^2)*y^4/(24*Mf*Nf^3);
BB=BB*180/pi;
l=y/(Nf*cos(BBf))-(1+2*tf^2+yitaf2)*y^3/(6*Nf^3*cos(BBf))+(5+28*tf^2+24*tf^4)*y^5/(120*Nf^5*cos(BBf));

这个公式里，计算孤长为什么会用这个公式不太明白

Gu_xl

zark 发表于 2012-7-5 08:25
版主，解压出问题了

现在论坛下载压缩文件好像都有问题！我重新上传了一个不压缩的文件！请去重新下载！
下面附件是 "高斯投影坐标正反算公式.pdf"供参考

游客，本帖隐藏的内容需要发帖数高于 10 才可浏览，你当前发帖数只有 0

zark

Gu_xl 发表于 2012-7-5 09:11
现在论坛下载压缩文件好像都有问题！我重新上传了一个不压缩的文件！请去重新下载！
下面附件是 "高斯投 ...

谢谢，也可以看看我的公式有没有用错

xiabin68

正在找这个，