📄 coord.java
字号:
x %= 1e6;
x -= 500000;
return c.convertFromTM( x, y, 0.0, merid, 1.0 );
}
/**
* The main program for the COORD class
*
*@param args The command line arguments
*/
public static void main( String[] args )
{
XY result = COORD.convertToGK( Double.valueOf( args[0] ).doubleValue(),
Double.valueOf( args[1] ).doubleValue() );
System.out.println( result.x + "," + result.y );
result = COORD.convertToLL( result.x, result.y );
System.out.println( result.x + "," + result.y );
System.exit( -1 );
double x = 49.4096;
double y = 8.7145;
double dist = 1e6;
double d;
int bk = -1;
int bi = -1;
COORD c = null;
for ( int k = 100; k < 101; k++ )
{
for ( int i = 28; i < 103; i++ )
{
result = COORD.convertToLatLong( Double.valueOf( args[0] ).doubleValue(),
Double.valueOf( args[1] ).doubleValue(),
k );
c = new COORD( k );
result = c.translate( true, result.x, result.y, i );
d = ( result.x - x ) * ( result.x - x ) + ( result.y - y ) * ( result.y - y );
if ( d < dist )
{
dist = d;
bi = i;
bk = k;
System.out.println( "X=" + result.x + " Y=" + result.y );
System.out.println( bi + "," + bk + "," + d );
}
}
}
System.out.println( bi + "," + bk + "," + dist );
System.out.println( c.gDatum[bi].name );
System.out.println( c.gDatum[bk].name );
System.exit( -1 );
}
/**
* Converts Gauss-Krueger to lat/long.
*
*@param x Parameter
*@param y Parameter
*@return Returned Value
*/
public static XY convertToLL( double x, double y )
{
XY gk = COORD.convertToLatLong( x, y );
return new COORD( 100 ).translate( false, gk.x, gk.y, 100 );
// return new COORD( 102 ).translate( false, gk.x, gk.y, 28 );
}
/**
* Converts lat/long to Gauss-Krueger.
*
*@param lat Parameter
*@param lon Parameter
*@return Returned Value
*/
public static XY convertToGK( double lat, double lon )
{
XY gk = new COORD( 102 ).translate( true, lat, lon, 100 );
// XY gk = COORD.translate( 100,28, lat, lon );
//System.out.println("PX="+gk.x+" PY="+gk.y);
return COORD.convertToGaussKrueger( gk.x, gk.y, 102 );
}
/**
* Retrieve earth datum parameters for given ID.
*
*@param datum Parameter
*@return Returned Value
*/
protected XY datumParams( int datum )
{
double f = 1.0 / gEllipsoid[gDatum[datum].ellipsoid].invf;
// flattening
double a;
// flattening
double es;
es = ( 2 * f - f * f );
// eccentricity^2
a = ( gEllipsoid[gDatum[datum].ellipsoid].a );
// semimajor axis
return new XY( a, es );
}
/**
* Performs a transmercator (TM) projection.
*
*@param lat Parameter
*@param lon Parameter
*@param lat0 Parameter
*@param lon0 Parameter
*@param k0 Parameter
*@return Returned Value
*/
protected XY convertToTM( double lat, double lon, double lat0, double lon0, double k0 )
{
double m;
double et2;
double n;
double t;
double c;
double A;
double a;
double m0;
double es;
double lambda;
double phi;
double lambda0;
double phi0;
XY result = new XY();
XY val = datumParams( datum );
a = val.x;
es = val.y;
lambda = lon * Degree;
phi = lat * Degree;
phi0 = lat0 * Degree;
lambda0 = lon0 * Degree;
m0 = M( phi0, a, es );
m = M( phi, a, es );
et2 = es / ( 1 - es );
n = a / Math.sqrt( 1 - es * Math.pow( Math.sin( phi ), 2.0 ) );
t = Math.pow( Math.tan( phi ), 2.0 );
c = et2 * Math.pow( Math.cos( phi ), 2.0 );
A = ( lambda - lambda0 ) * Math.cos( phi );
result.x = k0 * n * ( A + ( 1.0 - t + c ) * A * A * A / 6.0
+ ( 5.0 - 18.0 * t + t * t + 72.0 * c - 58.0 * et2 ) * Math.pow( A, 5.0 ) / 120.0 );
result.y = k0 * ( m - m0 + n * Math.tan( phi ) * ( A * A / 2.0
+ ( 5.0 - t + 9.0 * c + 4 * c * c ) * Math.pow( A, 4.0 ) / 24.0
+ ( 61.0 - 58.0 * t + t * t + 600.0 * c - 330.0 * et2 ) *
Math.pow( A, 6.0 ) / 720.0 ) );
return result;
}
/**
* Performs re-projection from TM.
*
*@param x Parameter
*@param y Parameter
*@param lat0 Parameter
*@param lon0 Parameter
*@param k0 Parameter
*@return Returned Value
*/
protected XY convertFromTM( double x, double y, double lat0, double lon0, double k0 )
{
double a;
double m0;
double es;
double et2;
double m;
double e1;
double mu;
double phi1;
double c1;
double t1;
double n1;
double r1;
double d;
double phi0;
double lambda0;
double lat;
double lon;
phi0 = lat0 * Degree;
lambda0 = lon0 * Degree;
//datumParams(gFilePrefs.datum, &a, &es);
XY val = datumParams( datum );
a = val.x;
es = val.y;
m0 = M( phi0, a, es );
et2 = es / ( 1.0 - es );
m = m0 + y / k0;
e1 = ( 1.0 - Math.sqrt( 1.0 - es ) ) / ( 1.0 + Math.sqrt( 1.0 - es ) );
mu = m / ( a * ( 1.0 - es / 4.0 - 3.0 * es * es / 64.0 - 5.0 * es * es * es / 256.0 ) );
phi1 = mu + ( 3.0 * e1 / 2.0 - 27.0 * Math.pow( e1, 3.0 ) / 32.0 ) * Math.sin( 2.0 * mu )
+ ( 21.0 * e1 * e1 / 16.0 - 55.0 * Math.pow( e1, 4.0 ) / 32.0 )
* Math.sin( 4.0 * mu ) + 151.0 * Math.pow( e1, 3.0 ) / 96.0 * Math.sin( 6.0 * mu )
+ 1097.0 * Math.pow( e1, 4.0 ) / 512.0 * Math.sin( 8.0 * mu );
c1 = et2 * Math.pow( Math.cos( phi1 ), 2.0 );
t1 = Math.pow( Math.tan( phi1 ), 2.0 );
n1 = a / Math.sqrt( 1 - es * Math.pow( Math.sin( phi1 ), 2.0 ) );
r1 = a * ( 1.0 - es ) / Math.pow( 1.0 - es * Math.pow( Math.sin( phi1 ), 2.0 ), 1.5 );
d = x / ( n1 * k0 );
lat = ( phi1 - n1 * Math.tan( phi1 ) / r1
* ( d * d / 2.0 - ( 5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1 * c1 - 9.0 * et2 )
* Math.pow( d, 4.0 ) / 24.0 + ( 61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1 * t1
- 252.0 * et2 - 3.0 * c1 * c1 ) * Math.pow( d, 6.0 ) / 720.0 ) ) / Degree;
lon = ( lambda0 + ( d - ( 1.0 + 2.0 * t1 + c1 ) * Math.pow( d, 3.0 ) / 6.0
+ ( 5.0 - 2.0 * c1 + 28.0 * t1 - 3.0 * c1 * c1 + 8.0 * et2 + 24.0 * t1 * t1 )
* Math.pow( d, 5.0 ) / 120.0 ) / Math.cos( phi1 ) ) / Degree;
return new XY( lat, lon );
}
/**
* Calculate M.
*
*@param phi Parameter
*@param a Parameter
*@param es Parameter
*@return Returned Value
*/
protected double M( double phi, double a, double es )
{
double fix;
if ( phi == 0.0 )
{
return 0.0;
}
else
{
fix = a * (
( 1.0 - es / 4.0 - 3.0 * es * es / 64.0 - 5.0 * es * es * es / 256.0 ) * phi -
( 3.0 * es / 8.0 + 3.0 * es * es / 32.0 + 45.0 * es * es * es / 1024.0 ) *
Math.sin( 2.0 * phi ) +
( 15.0 * es * es / 256.0 + 45.0 * es * es * es / 1024.0 ) * Math.sin( 4.0 * phi ) -
( 35.0 * es * es * es / 3072.0 ) * Math.sin( 6.0 * phi ) );
return fix;
}
}
static class XYZ
{
double X,Y,Z;
XYZ(double _x, double _y, double _z)
{
X=_x; Y=_y; Z=_z;
}
}
final static double awgs = 6378137.0; // WGS84 Semi-Major Axis = Equatorial Radius in meters
final static double bwgs = 6356752.314; // WGS84 Semi-Minor Axis = Polar Radius in meters
final static double abes = 6377397.155; // Bessel Semi-Major Axis = Equatorial Radius in meters
final static double bbes = 6356078.962; // Bessel Semi-Minor Axis = Polar Radius in meters
final static double cbes = 111120.6196; // Bessel latitude to Gauss-Krueger meters
final static double dx = -585.7; // Translation Parameter 1
final static double dy = -87.0; // Translation Parameter 2
final static double dz = -409.2; // Translation Parameter 3
final static double rotx = 2.540423689E-6; // Rotation Parameter 1
final static double roty = 7.514612057E-7; // Rotation Parameter 2
final static double rotz = -1.368144208E-5; // Rotation Parameter 3
final static double sc = 1/0.99999122; // Scaling Factor
static double l1;
static double b1;
static double l2;
static double b2;
static double h1;
static double h2;
static double R;
static double H;
static double a;
static double b;
static double eq;
static double N;
/*static double Xq;
static double Yq;
static double Zq;
static double X;
static double Y;
static double Z;
*/
public static XY convertToGK_NEW( double lat, double lon )
{
double l1=Math.PI*lon/180;
double b1=Math.PI*lat/180;
a=awgs;
b=bwgs;
eq=(a*a-b*b)/(a*a);
N=a/Math.sqrt(1-eq*Math.sin(b1)*Math.sin(b1));
double Xq=(N+h1)*Math.cos(b1)*Math.cos(l1);
double Yq=(N+h1)*Math.cos(b1)*Math.sin(l1);
double Zq=((1-eq)*N+h1)*Math.sin(b1);
XYZ xyz = HelmertTransformation(Xq,Yq,Zq);
//System.out.println("HX="+xyz.X+" HY="+xyz.Y+" HZ="+xyz.Z);
a=abes;
b=bbes;
eq=(a*a-b*b)/(a*a);
XYZ blh2 = BLRauenberg(xyz.X,xyz.Y,xyz.Z);
//System.out.println("RX="+blh2.X+" RY="+blh2.Y+" RZ="+blh2.Z);
XY RH = BesselBLnachGaussKrueger(blh2.X,blh2.Y);
b2=blh2.X*180/Math.PI;
l2=blh2.Y*180/Math.PI;
//System.out.println("PX="+b2+" PY="+l2);
return RH;
}
static XYZ HelmertTransformation(double x,double y,double z)
{
double xo,yo,zo;
xo=dx+(sc*(1*x+rotz*y-roty*z));
yo=dy+(sc*(-rotz*x+1*y+rotx*z));
zo=dz+(sc*(roty*x-rotx*y+1*z));
return new XYZ(xo,yo,zo);
}
static XY BesselBLnachGaussKrueger(double b,double ll)
{
double l;
double l0;
double bg;
double lng;
double Ng;
double k;
double t;
double eq;
double Vq;
double v;
double nk;
double X;
double gg;
double SS;
double Y;
double kk;
double pii;
double RVV;
double Re,Ho;
bg=180*b/Math.PI;
lng=180*ll/Math.PI;
l0=3*Math.round((180*ll/Math.PI)/3);
l0=Math.PI*l0/180;
l=ll-l0;
k=Math.cos(b);
t=Math.sin(b)/k;
eq=(abes*abes-bbes*bbes)/(bbes*bbes);
Vq=1+eq*k*k;
v=Math.sqrt(Vq);
Ng=abes*abes/(bbes*v);
nk=(abes-bbes)/(abes+bbes);
X=((Ng*t*k*k*l*l)/2)+((Ng*t*(9*Vq-t*t-4)*k*k*k*k*l*l*l*l)/24);
gg=b+(((-3*nk/2)+(9*nk*nk*nk/16))*Math.sin(2*b)+15*nk*nk*Math.sin(4*b)/16-35*nk*nk*nk*Math.sin(6*b)/48);
SS=gg*180*cbes/Math.PI;
Ho=(SS+X);
Y=Ng*k*l+Ng*(Vq-t*t)*k*k*k*l*l*l/6+Ng*(5-18*t*t+t*t*t*t)*k*k*k*k*k*l*l*l*l*l/120;
kk=500000;
pii=Math.PI;
RVV=Math.round((180*ll/pii)/3);
Re=RVV*1000000+kk+Y;
return new XY(Re,Ho);
}
static XYZ BLRauenberg (double x,double y,double z)
{
double f;
double f1;
double f2;
double ft;
double p;
double b,l,h;
f=Math.PI*50/180;
p=z/Math.sqrt(x*x+y*y);
do
{
f1=neuF(f,x,y,p);
f2=f;
f=f1;
ft=180*f1/Math.PI;
}
while(!(Math.abs(f2-f1)<10E-10));
b=f;
l=Math.atan(y/x);
h=Math.sqrt(x*x+y*y)/Math.cos(f1)-a/Math.sqrt(1-eq*Math.sin(f1)*Math.sin(f1));
return new XYZ(b,l,h);
}
static double neuF(double f,double x,double y,double p)
{
double zw;
double nnq;
zw=a/Math.sqrt(1-eq*Math.sin(f)*Math.sin(f));
nnq=1-eq*zw/(Math.sqrt(x*x+y*y)/Math.cos(f));
return(Math.atan(p/nnq));
}
double round(double src)
{
return Math.round(src);
/* double theInteger;
double theFraction;
double criterion = 0.5;
theFraction = modf(src,&theInteger);
if (!(theFraction < criterion))
{
theInteger += 1;
}
return theInteger;
*/
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -