GIS Coordinate Transformation and programming implementation

Source: Internet
Author: User
Tags gety

Recently, we have been working on this coordinate transformation task, which involves many profound things such as surveyors. Here I will not explain these difficult theories. Here, I will post the code to share with you. In fact, if you want to write this code, you must be familiar with the knowledge related to the land survey. Otherwise, you will not be able to understand the code. Okay.

Source code

/**
* Space earth Cartesian coordinates-> Earth coordinates
*/
Public geopoint xyz_blh (INT ellipse, point3d point)
{
Geopoint = new geopoint ();
If (geopoint = NULL | point = NULL)
{
System. Out. println ("the object is blank! ");
}
Double A = 0, B = 0, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the square of the first eccentric Heart Rate
Double DETA = 0.0000001;
Switch (ellipse) // select an elliptical body
{
Case 0: // WGS84 elliptical body
{
A = 6378137.0; // long Half Axis
B = 6356752.3142; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = 0.006694384999588;
Break;
}
Case 1: // Beijing 54 elliptical
{
A = 6378245.0; // long Half Axis
B = 6356863.0187730473; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = 0.006693421622966;
Break;
}
Case 2: // Xi'an 80 elliptical
{
A = 6378140.0; // long Half Axis
B = 6356755.2881575287; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = 0.006694384999588;
Break;
}
Default: break;
}
// Double W = math. SQRT (1.0-E1 * E1 * Math. Pow (math. Sin (point. getx (), 2 ));
// Double N = A/W;
Double X = point. getx ();
Double Y = point. Gety ();
Double Z = point. Getz ();
Double M = math. SQRT (math. Pow (x, 2) + math. Pow (Y, 2 ));
Double L = geopoint. getl (); // geopoint longitude
Double B = geopoint. getb (); //
Double H = geopoint. geth (); // large height
L = math. atan (y/X );
If (L <0)
{
L + = math. Pi;
}
Double E2 _ = e2/(1-e2 );
Double C = A * Math. SQRT (1 + E2 _);
Double CE2 = C * E2;
Double K = 1 + E2 _;
Double front = z/m;
Double temp = front;
Int COUNT = 0; // number of iterations
Do
{
Front = temp;
M = math. SQRT (math. Pow (x, 2) + math. Pow (Y, 2 ));
Temp = z/m + CE2 * front/(m * Math. SQRT (K + math. Pow (front, 2 )));
}
While (math. Abs (front-Temp) <DETA & count <100000); // whether it is within the error range
B = math. atan (temp); // obtain the latitude.
If (B <0)
{
B + = math. Pi;
}
Double W = math. SQRT (1-e1 * E1 * Math. Sin (B) * Math. Sin (B ));
Double N = A/W;
// N = (A * m-C * C)/(2 * B * z );
System. Out. println ("n =" + n );
H = m/Math. Cos (B)-N; // height
// H = z/Math. Sin (B)-N * (1.0-E1 * E1 );
// H = x/(math. Cos (B) * Math. Cos (L)-N;
// Convert to the angle value
L = math. todegrees (L );
B = math. todegrees (B );
// H = math. todegrees (h );
Geopoint. SETl (L );
Geopoint. SETB (B );
Geopoint. Seth (h );
Return geopoint;
}
 
/**
* Coordinate-> space ry Cartesian coordinates
* @ Param ellipse elliptical body
* @ Param geopoint coordinate before conversion
* @ Return returns the Cartesian coordinates of the space.
*/
Public point3d blh_xyz (INT ellipse, geopoint)
{
Point3d point1 = new point3d (1, 1 );
If (geopoint = NULL | point1 = NULL)
{
System. Out. println ("the object is blank! ");
}
Double A = 0, B, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the square of the first eccentric Heart Rate
Switch (ellipse) // select an elliptical body
{
Case 0: // WGS84 elliptical body
{
A = 6378137.0; // long Half Axis
B = 6356752.3142; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = 0.006694384999588;
Break;
}
Case 1: // Beijing 54 elliptical
{
A = 6378245.0; // long Half Axis
B = 6356863.0187730473; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = 0.006693421622966;
Break;
}
Case 2: // Xi'an 80 elliptical
{
A = 6378140.0; // long Half Axis
B = 6356755.2881575287; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = 0.006694384999588;
Break;
}
Default: break;
}
Double B = geopoint. getb () * Math. PI/180; // converts to radians.
Double L = geopoint. getl () * Math. PI/180;
Double H = geopoint. geth ();
// Calculate the radius of curvature of a uniform Ring
Double N = A/Math. SQRT (1.0-E1 * E1 * Math. Pow (math. Sin (B), 2 ));
// Calculate space Cartesian coordinates
Double X = (n + H) * Math. Cos (B) * Math. Cos (L );
Double Y = (n + H) * Math. Cos (B) * Math. Sin (L );
Double Z = (N * (1-e1 * E1) + H) * Math. Sin (B );
Point1.x = X;
Point1.y = y;
Point1.z = z;
Return point1;
}
 
/**
* Gaussian-kerluge coordinate normal calculation (from the earth coordinate to the plane Cartesian coordinate)
* @ Param ellipse elliptical body
* @ Param zonewide bandwidth (3 or 6 degrees)
* @ Param geopoint ry
* @ Param point Cartesian coordinates
*/
Public void bl_xy (INT ellipse, int zonewide, geopoint, point)
{
If (geopoint = NULL | point = NULL)
{
System. Out. println ("the object is blank! ");
}
Double A = 0, B, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the square of the first eccentric Heart Rate
Switch (ellipse) // select an elliptical body
{
Case 0: // WGS84 elliptical body
{
A = 6378137.0; // long Half Axis
B = 6356752.3142; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B;
// E2e = 0.006694384999588;
Break;
}
Case 1: // Beijing 54 elliptical
{
A = 6378245.0; // long Half Axis
B = 6356863.0187730473; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B;
// E2e = 0.006693421622966;
Break;
}
Case 2: // Xi'an 80 elliptical
{
A = 6378140.0; // long Half Axis
B = 6356755.2881575287; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B;
// E2e = 0.006694384999588;
Break;
}
Default: break;
}

Double B = geopoint. getb ();
Double L = geopoint. getl ();
Double centerl = 0.0; // degree of projection with central longitude
Int n = 0; // The belt number of the projection band (6 ° band: 1-60, 3 ° band: 1-120)
Double L = 0; // The Difference Between the longitude of the region and the central longitude
If (6 = zonewide) // if it is a 6-degree band
{
N = (INT) L/6;
If (L % 6> 0)
{
N = n + 1;
}
Centerl = 6 * n-3;
L = L-centerl;
}
If (3 = zonewide) // if it is a 3-degree band
{
N = (INT) (L-1.5)/3;
If (L-1.5) % 3> 0)
{
N = n + 1;
}
Centerl = 3 * N;
L = L-centerl;
}
// Convert to radians
B = B * Math. PI/180;
L = L * Math. PI/180;
L = L * Math. PI/180;
Double X; // the meridian arc length calculated from the equator
// Calculate the coefficient of meridian Arc Length
Double a0 = 1.0 + 3.0/4.0 * Math. Pow (E1, 2) + 45.0/64.0 * Math. Pow (E1, 4)
+ 350.0/512.0 * Math. Pow (E1, 6) + 11025.0/16384.0 * Math. Pow (E1, 8 );
Double a2 =-3.0/4 * Math. Pow (E1, 2) + 60.0/64 * Math. Pow (E1, 4) +
525.0/512 * Math. Pow (E1, 6) + 17640.0/16384.0 * Math. Pow (E1, 8)/2.0;
Double A4 = 15.0/64.0 * Math. Pow (E1, 4) + 210.0/512.0 * Math. Pow (E1, 6) +
8820.0/16384.0 * Math. Pow (E1, 8)/4.0;
Double A6 =-35.0/512.0 * Math. Pow (E1, 6) + 2520.0/16384.0 * Math. Pow (E1, 8)/6.0;
Double A8 = 315.0/16384.0 * Math. Pow (E1, 8)/8.0;
// Calculate the meridian arc length x
X = A * (1.0-math. pow (E1, 2) * a0 * B + A2 * Math. sin (2 * B) + A4 * Math. sin (4 * B)
+ A6 * Math. Sin (6 * B) + A8 * Math. Sin (8 * B );
Double T = math. Tan (B );
Double Anke = e2 * Math. Cos (B );
Double N = A/Math. SQRT (1.0-math. Pow (E1, 2) * Math. Pow (math. Sin (B), 2 ));
// Coordinate Calculation
Double X = point. getx ();
Double Y = point. Gety ();
X = x + 1.0/2.0 * n * T * Math. Pow (math. Cos (B), 2) * Math. Pow (L, 2) +
1.0/24.0 * n * T * (5.0-math. pow (T, 2) + 9 * Math. pow (Anke, 2) + 4 * Math. pow (Anke, 4) * Math. pow (math. cos (B), 4) * Math. pow (L, 4)
+ 1.0/720 * n * T * (61.0-58 * Math. pow (T, 2) + math. pow (T, 4) + 270.0 * Math. pow (Anke, 2)-330.0 * Math. pow (Anke, 2) * Math. pow (T, 2) * Math. pow (math. cos (B), 6 );

Y = N * Math. cos (B) * L + 1.0/6 * n * (1.0-math. pow (T, 2) + math. pow (Anke, 2) * Math. pow (math. cos (B), 3) * Math. pow (L, 3) +
1.0/120 * n * (5.0-18 * T * t + math. pow (T, 4) + 14 * Math. pow (Anke, 2)-58 * Anke * T * t )*
Math. Pow (math. Cos (B), 3) * Math. Pow (L, 5 );
Y + = 500000.0; // Add 500 km
Point. setx (X );
Point. sety (y );
}
 
/**
* The Cartesian coordinate transformation of the Gaussian plane is the earth coordinate.
* @ Param ellipse elliptical body
* @ Param zonewide bandwidth
* @ Param point plane coordinate
* @ Param geopoint ry
*/
Public void xy_bl (INT ellipse, double centerl, point, geopoint)
{
If (geopoint = NULL | point = NULL)
{
System. Out. println ("the object is blank! ");
}
Double A = 0, B, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the second eccentric Heart Rate
Switch (ellipse) // select an elliptical body
{
Case 0: // WGS84 elliptical body
{
A = 6378137.0; // long Half Axis
B = 6356752.3142; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B;
// E2e = 0.006694384999588;
Break;
}
Case 1: // Beijing 54 elliptical
{
A = 6378245.0; // long Half Axis
B = 6356863.0187730473; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B;
// E2e = 0.006693421622966;
Break;
}
Case 2: // Xi'an 80 elliptical
{
A = 6378140.0; // long Half Axis
B = 6356755.2881575287; // short half axis
E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/;
E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B;
// E2e = 0.006694384999588;
Break;
}
Default: break;
}
Double X = point. getx ();
Double Y = point. Gety ()-500000;
Double BF = 0.0; // latitude of the base point
Double b0 = 1.0 + 3.0/4.0 * E1 * E1 + 45.0/64.0 * Math. Pow (E1, 4) + 350.0/312.0 * Math. Pow (E1, 6) +
11025.0/16384.0 * Math. Pow (E1, 8) + 43659.0/65536.0 * Math. Pow (E1, 10 );
BF = x/(A * (1.0-math. Pow (E1, 2) * B0 );
Double TF = math. Tan (BF );
Double mf = A * (1.0-E1 * E1)/Math. SQRT (math. Pow (1.0-math. Pow (E1, 2), 3 ));
Double NF = (A/math. SQRT (1.0-math. pow (E1, 2)/math. SQRT (1.0 + math. pow (e2 * Math. cos (BF), 2 ));
Double Anke = e2 * Math. Cos (BF );
Double B = geopoint. getb ();
Double L = geopoint. getl ();
// Start Coordinate Calculation
B = BF-TF/(2 * MF * NF * Math. cos (BF) * Math. pow (Y, 2) + TF/(24 * MF * Math. pow (NF, 3 ))*
(5.0 + 3.0 * tF * TF + Anke * anke-9.0 * Anke * tF) * Math. Pow (Y, 4 )-
1.0/(720.0 * Math. pow (NF, 5) * Math. cos (BF) * (61.0 + 90.0 * tF * TF + 45 * Math. pow (TF, 4) * Math. pow (Y, 6 );
L = y/(NF * Math. cos (BF)-(1.0 + 2 * tF * TF + Anke * Anke) * Math. pow (Y, 3)/(6.0 * Math. pow (NF, 3) * Math. cos (BF ))
+ (5.0 + 28 * tF * TF + 24 * Math. pow (TF, 4) + 6 * Anke + 8 * Anke * tF) * Math. pow (Y, 5)/(120.0 * Math. pow (NF, 5) * Math. cos (BF ));

B = B * 180/Math. Pi;
L = L * 180/Math. Pi;
L = L + centerl; // central longitude + longitude difference
Geopoint. SETB (B );
Geopoint. SETl (L );
}

The possible reason for these codes is that, when the space Cartesian coordinates are converted to the Earth coordinates, the calculation of the Earth height is very different.

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.