Home > Others

GIS Coordinate Transformation and programming implementation

Source: Internet
Author: User
Tags gety
Developer on Alibaba Coud: Build your first app with APIs, SDKs, and tutorials on the Alibaba Cloud. Read more ＞

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.

This article is an English version of an article which is originally in the Chinese language on aliyun.com and is provided for information purposes only. This website makes no representation or warranty of any kind, either expressed or implied, as to the accuracy, completeness ownership or reliability of the article or any translations thereof. If you have any concerns or complaints relating to the article, please send an email, providing a detailed description of the concern or complaint, to info-contact@alibabacloud.com. A staff member will contact you within 5 working days. Once verified, infringing content will be removed immediately.

Related Keywords:
discovering gis and arcgis interface and implementation in java difference between interface and implementation security policies and implementation issues pdf download programming quotes and sayings visual and game programming coding and programming salary

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.

What's Trending
Top 10 Tags
datastax versions naming convention zookeeper client class definition md5 microsoft sql server 2005 data structures exception handling error handling
Top 10 Keywords
microsoft download center down wordpress address url site address url wordpress address url windows installer 4 0 download 302 not found web address url definition site address url wordpress db2 integer mac os installation step by step pdf abbreviation for return

A Free Trial That Lets You Build Big!

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

Get Started for Free

  • Sales Support

    1 on 1 presale consultation

    Chat Contact Sales

  • After-Sales Support

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

    Open a Ticket

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

    Learn More