正算公式(將經緯度轉化為座標):
java代碼(附有源代碼和修改後的代碼):
源代碼:
<code>/**
* 由經緯度反算成高斯投影座標
*
* @param longitude
* @param latitude
* @return
*/
public static double[] GaussToBLToGauss(
\t\t\t\t\t\t\t\tdouble longitude, double latitude)
{
\tint ProjNo = 0; int ZoneWide; // //帶寬
\tdouble[] output = new double[2];
\tdouble longitude1, latitude1, longitude0, X0, Y0, xval, yval;
\tdouble a, f, e2, ee, NN, T, C, A, M, iPI;
\tiPI = 0.0174532925199433; // //3.1415926535898/180.0;
\tZoneWide = 6; // //6度帶寬
\ta = 6378245.0;
\tf = 1.0 / 298.3; // 54年北京座標系參數
\t// //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;
}/<code>
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);
修改後的代碼:
<code>/**
\t * 由經緯度反算成高斯投影座標
*
* @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;
}/<code>
評論：而且xval, yval和公式中的X, Y正好反過來了
分享到:
閱讀更多 愛碼廝小妖 的文章
