用LISP编写一个“高斯正算”（将经纬度转成平面坐标）函数
;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 ;中央子午线
H0.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))
)
程序执行结果
笑脸转换结果
东坐标都不差，就北坐标差这么多，不知道问题出现在哪？
请高手指点一下。谢谢
http://bbs.mjtd.com/xwb/images/bgimg/icon_logo.png 该贴已经同步到 zark的微博 本帖最后由 004 于 2013-4-11 16:33 编辑
革天明 发表于 2012-7-6 08:09 http://bbs.mjtd.com/static/image/common/back.gif
笑脸是如何加上去的，不是DCL吧，VB？
笑脸----coordMG是一款应用广泛的测量坐标转换软件，和DCL无关。 感谢分享，回复学习。 自己先顶一下 建议不用这个软件对比.用苏某某的一个EXCEL版软件对比
对比一下！
http://apptools.mjtd.com/apptools.php?mod=view&apptoolsid=36 1000[粉] = 元 [羞羞] ①78595⑧333
http://bbs.mjtd.com/xwb/images/bgimg/icon_logo.png 来自 昔虚宦疙陀 的新浪微博 Gu_xl 发表于 2012-7-4 22:53 static/image/common/back.gif
对比一下！
http://apptools.mjtd.com/apptools.php?mod=view&apptoolsid=36
版主，解压出问题了 高斯坐标正反matlab算程序：
function =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 =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));
这个公式里，计算孤长为什么会用这个公式不太明白 zark 发表于 2012-7-5 08:25 static/image/common/back.gif
版主，解压出问题了
现在论坛下载压缩文件好像都有问题！我重新上传了一个不压缩的文件！请去重新下载！
下面附件是 "高斯投影坐标正反算公式.pdf"供参考
Gu_xl 发表于 2012-7-5 09:11 static/image/common/back.gif
现在论坛下载压缩文件好像都有问题！我重新上传了一个不压缩的文件！请去重新下载！
下面附件是 "高斯投 ...
谢谢，也可以看看我的公式有没有用错 正在找这个，
页:
[1]
2