近日工程项目涉及基本比例尺分幅,用了CASS11的相关功能,发现有些缺陷,就打算自己整一个,期间遇到些问题,记载下来备查。
cass11见上图,不妥之处有以下几点:
1.界面比较乱,条理不清。
2.“设置中央线”没有必要,按规范即可。
3.缺失10万、25万、50万比例尺。
4.无大地坐标输出方案。
重新设计后的程序界面见下图:
分为参数设置、数据录入两个区,参数设置全部用选择框,界面简洁、清晰。程序编写过程中遇到的问题及处理方式:
1.当输入的生成范围跨多个带区时,采用大地坐标没有影响,采用投影坐标需分别计算每个带区的中央子午线及其经区范围,注意处理范围左边缘带区的最小经度、右边缘带区的最大经度范围。另外3度带的50万,100万比例尺的中央线也必须按6度带进行计算,这样中央子午线两边的投影变形一致。
2.高斯投影正算就不想浪费时间了,在CSDN检索到“氧气女神”的代码,程序结构很棒,BUT,计算结果是错的,私信了也不回话,这妹子估计比较忙。这是用她的数据及她的程序计算结果。“https://blog.csdn.net/qq_44718128/article/details/117874196”
后来检索到“卤雅少年”的代码“https://blog.csdn.net/Zj1638/article/details/125380593”,将两人的代码整合了一下,得到高斯投影正算代码,精度为0.0005,先将就用吧,哪位兄台有更高精度的代码请发一个,谢谢!我的邮箱:huaheroys@126.com
#region 高斯投影正算 纬度,经度,中央子午线,代号,是否加带号大数 单位度
public static double[] BLtoXY(double B,double L,double L0,double dh,bool jdh) {
Editor ed = Autodesk.AutoCAD.ApplicationServices.Application.DocumentManager.MdiActiveDocument.Editor;
//CGCS2000椭球参数
double a = 6378137.0;
double f = 0.00335281066475;
double e1 = 0.00669437999014631;
double e2 = 0.00673949674228149;
double m0, m2, m4, m6, m8;
double a0, a2, a4, a6, a8;
//m0-m8的值
m0 = a * (1 - e1);
m2 = 1.5 * m0 * e1;
m4 = 1.25 * m2 * e1;
m6 = (7.0 / 6) * m4 * e1;
m8 = (9.0 / 8) * m6 * e1;
//a0-a8的值
a0 = m0 + m2 / 2.0 + (3.0 / 8) * m4 + (5.0 / 16) * m6 + (35.0 / 128) * m8;
a2 = (1.0 / 2) * (m2 + m4) + (15.0 / 32) * m6 + (7.0 / 16) * m8;
a4 = (1.0 / 8) * m4 + (3.0 / 16) * m6 + (7.0 / 32) * m8;
a6 = (1.0 / 32) * m6 + (1.0 / 16) * m8;
a8 = (1.0 / 128) * m8;
//求椭球面上某点到赤道的子午线弧长,精度至0.001m
double radB = (Math.PI / 180) * B; //将纬度转换为弧度
double Δ1, Δ2, Δ3, Δ4, Δ5;
Δ1 = a0 * radB;
Δ2 = (a2 / 2) * Math.Sin(2 * radB);
Δ3 = (a4 / 4) * Math.Sin(4 * radB);
Δ4 = (a6 / 6) * Math.Sin(6 * radB);
Δ5 = (a8 / 8) * Math.Sin(8 * radB);
double X = Math.Round(Δ1 - Δ2 + Δ3 - Δ4 + Δ5, 8);
//卯酉圈曲率半径N
double N = a / Math.Sqrt(1 - e1 * Math.Sin(radB) * Math.Sin(radB));
// tan(B)的值
double t = Math.Tan(radB);
// η的值
double η = Math.Sqrt(e2) * Math.Cos(radB);
//l的值
double ρm = 180 * 60 * 60 / Math.PI;
double l = ((L - L0) * 3600) / ρm;
//高斯坐标中的x值和y值
double x = Math.Round(X + (N / 2) * t * Math.Cos(radB) * Math.Cos(radB) * l * l + (N / 24) * t * (5 - t * t + 9 * η * η + 4 * Math.Pow(η, 4)) * Math.Pow(Math.Cos(radB), 4) * Math.Pow(l, 4) + (N / 720) * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(l, 6) * (Math.Pow(Math.Cos(radB), 6)), 4);
double y = 500000 + Math.Round(N * Math.Cos(radB) * l + (N / 6.0) * (1 - t * t + η * η) * Math.Pow(Math.Cos(radB), 3) * Math.Pow(l, 3) + (N / 120.0) * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * η * η - 58 * η * η * t * t) * Math.Pow(Math.Cos(radB), 5) * Math.Pow(l, 5), 4);
if (jdh) { y = y + dh * 1000000; }
return new double[] { y, x };
}
#endregion
程序运行效果-大地坐标
程序运行效果-投影坐标
欢迎批评指导,QQ群:861650360.