高斯投影正算公式以及java代码
2011-12-16 11:04
260 查看
正算公式(将经纬度转化为坐标):
![](http://blog.51cto.com/attachment/201112/110706360.png)
![](http://blog.51cto.com/attachment/201112/110751678.png)
![](http://blog.51cto.com/attachment/201112/110835917.png)
java代码(附有源代码和修改后的代码):
源代码:
/**
* 由经纬度反算成高斯投影坐标
*
* @param longitude
* @param latitude
* @return
*/
public static double[] GaussToBLToGauss(double longitude, double latitude) {
int ProjNo = 0;
int ZoneWide; // //带宽
double[] output = new double[2];
double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
double a, f, e2, ee, NN, T, C, A, M, iPI;
iPI = 0.0174532925199433; // //3.1415926535898/180.0;
ZoneWide = 6; // //6度带宽
a = 6378245.0;
f = 1.0 / 298.3; // 54年北京坐标系参数
// //a=6378140.0; f=1/298.257; //80年西安坐标系参数
ProjNo = (int) (longitude / ZoneWide);
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI;
longitude1 = longitude * iPI; // 经度转换为弧度
latitude1 = latitude * iPI; // 纬度转换为弧度
e2 = 2 * f - f * f;
ee = e2 * (1.0 - e2);
NN = a
/ Math.sqrt(1.0 - e2 * Math.sin(latitude1)
* Math.sin(latitude1));
T = Math.tan(latitude1) * Math.tan(latitude1);
C = ee * Math.cos(latitude1) * Math.cos(latitude1);
A = (longitude1 - longitude0) * Math.cos(latitude1);
M = a
* ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256)
* latitude1
- (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2
/ 1024) * Math.sin(2 * latitude1)
+ (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)
* Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072)
* Math.sin(6 * latitude1));
// 因为是以赤道为Y轴的,与我们南北为Y轴是相反的,所以xy与高斯投影的标准xy正好相反;
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
* Math.tan(latitude1)
* (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 = 1000000L * (ProjNo + 1) + 500000L;
Y0 = 0;
xval = xval + X0;
yval = yval + Y0;
output[0] = xval;
output[1] = yval;
return output;
}
java代码是网上找的, 原作者对该代码很有自信, 我根据与正算公式的比较, 发现了几个不同点, 对代码做了修改.
不同点1:ee不同
代码中的"ee = e2 * (1.0 - e2)",这对应了正算公式中的e'的平方, 代码中的f就是正算公式中的扁率α. 计算ee是否与e'的平方一致, e的平方=(a²-b²)/a², e'的平方=(a²-b²)/b², 过程如下:
一. e2=2*f-f*f=2*(a-b)/a-(a-b)/a * (a-b)/a=(a²-b²)/a²;// e2=公式中e的平方,正确
二. ee=e2*(1.0-e2)=(a²-b²)b²/(a² * a²);//ee!=e'的平方
不同点2: 代码与公式某常量值不同
公式中的X=.....(...+270C-330TC)...; Y=....(...+14C-58TC)....;
对应代码分别为...(...+600*C-330*ee)..; ...(...+72*C-58*ee);
修改后的代码:
/**
* 由经纬度反算成高斯投影坐标
*
* @param longitude
* @param latitude
* @return
*/
public static double[] GaussToBLToGauss(double longitude, double latitude) {
int ProjNo = 0;
int ZoneWide; // //带宽
double[] output = new double[2];
double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
double a, f, e2, ee, NN, T, C, A, M, iPI;
iPI = 0.0174532925199433; // //3.1415926535898/180.0;
ZoneWide = 6; // //6度带宽
a = 6378245.0;
f = 1.0 / 298.3; // 54年北京坐标系参数
// //a=6378140.0; f=1/298.257; //80年西安坐标系参数
ProjNo = (int) (longitude / ZoneWide);
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI;
longitude1 = longitude * iPI; // 经度转换为弧度
latitude1 = latitude * iPI; // 纬度转换为弧度
e2 = 2 * f - f * f;
ee = e2 / (1.0 - e2);
NN = a
/ Math.sqrt(1.0 - e2 * Math.sin(latitude1)
* Math.sin(latitude1));
T = Math.tan(latitude1) * Math.tan(latitude1);
C = ee * Math.cos(latitude1) * Math.cos(latitude1);
A = (longitude1 - longitude0) * Math.cos(latitude1);
M = a
* ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256)
* latitude1
- (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2
/ 1024) * Math.sin(2 * latitude1)
+ (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)
* Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072)
* Math.sin(6 * latitude1));
// 因为是以赤道为Y轴的,与我们南北为Y轴是相反的,所以xy与高斯投影的标准xy正好相反;
xval = NN
* (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 14
* C - 58 * ee)
* A * A * A * A * A / 120);
yval = M
+ NN
* Math.tan(latitude1)
* (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61
- 58 * T + T * T + 270 * C - 330 * ee)
* A * A * A * A * A * A / 720);
X0 = 1000000L * (ProjNo + 1) + 500000L;
Y0 = 0;
xval = xval + X0;
yval = yval + Y0;
output[0] = xval;
output[1] = yval;
return output;
}
![](http://blog.51cto.com/attachment/201112/110706360.png)
![](http://blog.51cto.com/attachment/201112/110751678.png)
![](http://blog.51cto.com/attachment/201112/110835917.png)
java代码(附有源代码和修改后的代码):
源代码:
/**
* 由经纬度反算成高斯投影坐标
*
* @param longitude
* @param latitude
* @return
*/
public static double[] GaussToBLToGauss(double longitude, double latitude) {
int ProjNo = 0;
int ZoneWide; // //带宽
double[] output = new double[2];
double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
double a, f, e2, ee, NN, T, C, A, M, iPI;
iPI = 0.0174532925199433; // //3.1415926535898/180.0;
ZoneWide = 6; // //6度带宽
a = 6378245.0;
f = 1.0 / 298.3; // 54年北京坐标系参数
// //a=6378140.0; f=1/298.257; //80年西安坐标系参数
ProjNo = (int) (longitude / ZoneWide);
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI;
longitude1 = longitude * iPI; // 经度转换为弧度
latitude1 = latitude * iPI; // 纬度转换为弧度
e2 = 2 * f - f * f;
ee = e2 * (1.0 - e2);
NN = a
/ Math.sqrt(1.0 - e2 * Math.sin(latitude1)
* Math.sin(latitude1));
T = Math.tan(latitude1) * Math.tan(latitude1);
C = ee * Math.cos(latitude1) * Math.cos(latitude1);
A = (longitude1 - longitude0) * Math.cos(latitude1);
M = a
* ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256)
* latitude1
- (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2
/ 1024) * Math.sin(2 * latitude1)
+ (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)
* Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072)
* Math.sin(6 * latitude1));
// 因为是以赤道为Y轴的,与我们南北为Y轴是相反的,所以xy与高斯投影的标准xy正好相反;
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
* Math.tan(latitude1)
* (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 = 1000000L * (ProjNo + 1) + 500000L;
Y0 = 0;
xval = xval + X0;
yval = yval + Y0;
output[0] = xval;
output[1] = yval;
return output;
}
java代码是网上找的, 原作者对该代码很有自信, 我根据与正算公式的比较, 发现了几个不同点, 对代码做了修改.
不同点1:ee不同
代码中的"ee = e2 * (1.0 - e2)",这对应了正算公式中的e'的平方, 代码中的f就是正算公式中的扁率α. 计算ee是否与e'的平方一致, e的平方=(a²-b²)/a², e'的平方=(a²-b²)/b², 过程如下:
一. e2=2*f-f*f=2*(a-b)/a-(a-b)/a * (a-b)/a=(a²-b²)/a²;// e2=公式中e的平方,正确
二. ee=e2*(1.0-e2)=(a²-b²)b²/(a² * a²);//ee!=e'的平方
不同点2: 代码与公式某常量值不同
公式中的X=.....(...+270C-330TC)...; Y=....(...+14C-58TC)....;
对应代码分别为...(...+600*C-330*ee)..; ...(...+72*C-58*ee);
修改后的代码:
/**
* 由经纬度反算成高斯投影坐标
*
* @param longitude
* @param latitude
* @return
*/
public static double[] GaussToBLToGauss(double longitude, double latitude) {
int ProjNo = 0;
int ZoneWide; // //带宽
double[] output = new double[2];
double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
double a, f, e2, ee, NN, T, C, A, M, iPI;
iPI = 0.0174532925199433; // //3.1415926535898/180.0;
ZoneWide = 6; // //6度带宽
a = 6378245.0;
f = 1.0 / 298.3; // 54年北京坐标系参数
// //a=6378140.0; f=1/298.257; //80年西安坐标系参数
ProjNo = (int) (longitude / ZoneWide);
longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
longitude0 = longitude0 * iPI;
longitude1 = longitude * iPI; // 经度转换为弧度
latitude1 = latitude * iPI; // 纬度转换为弧度
e2 = 2 * f - f * f;
ee = e2 / (1.0 - e2);
NN = a
/ Math.sqrt(1.0 - e2 * Math.sin(latitude1)
* Math.sin(latitude1));
T = Math.tan(latitude1) * Math.tan(latitude1);
C = ee * Math.cos(latitude1) * Math.cos(latitude1);
A = (longitude1 - longitude0) * Math.cos(latitude1);
M = a
* ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256)
* latitude1
- (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2
/ 1024) * Math.sin(2 * latitude1)
+ (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)
* Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072)
* Math.sin(6 * latitude1));
// 因为是以赤道为Y轴的,与我们南北为Y轴是相反的,所以xy与高斯投影的标准xy正好相反;
xval = NN
* (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 14
* C - 58 * ee)
* A * A * A * A * A / 120);
yval = M
+ NN
* Math.tan(latitude1)
* (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61
- 58 * T + T * T + 270 * C - 330 * ee)
* A * A * A * A * A * A / 720);
X0 = 1000000L * (ProjNo + 1) + 500000L;
Y0 = 0;
xval = xval + X0;
yval = yval + Y0;
output[0] = xval;
output[1] = yval;
return output;
}
相关文章推荐
- 手把手教你如何把java代码,打包成jar文件以及转换为exe可执行文件
- Android使用Java代码设置selector或drawable,以及使用自定义控件方式使用它
- BDF2-JAVA代码以及view页面获取当前登录用户的信息
- java S2SH项目框架Struts2.3+Spring3.1+Hibernate3.6框架搭建整合以及实例demo代码分享
- 大型JavaWeb分布式系统中关于maven多模块构建以及代码依赖管理
- Spring学习笔记1——IOC: 尽量使用注解以及java代码(转)
- Dijkstra算法以及java实现_02(代码部分)
- 二叉树定义及相关术语、节点数计算公式、代码实现(遍历,Java版)
- CharacterEncodingFilter.java以及修改web.xml所用代码
- 大型JavaWeb分布式系统中关于maven多模块构建以及代码依赖管理
- java 中的构造代码块和构造函数以及静态代码块
- 关于java代码中静态代码块、非静态构造代码块以及构造方法的执行顺序
- Java笔记3 多线程<2>线程间通信-代码分析以及多线程常见方法的运用
- 统计某个路径下所有的java文件,以及统计代码数量
- {@link}与@see的简单使用以及区别--如何写好Java代码注释。
- JNI详解以及使用Java调c代码
- java POI导出excel(带有图片的数据)以及ftp上传图片到另一台服务器的方法参考代码
- 手把手教你如何把java代码,打包成jar文件以及转换为exe可执行文件
- Eclipse中的快捷键快速生成常用代码(例如无参、带参构造,set、get方法),以及Java中重要的内存分析(栈、堆、方法区、常量池)
- 手把手教你如何把java代码,打包成jar文件以及转换为exe可执行文件