ezprojimpl.pas

很管用的GIS控件

Unit EzProjImpl;

{***********************************************************}
{     EzGIS/CAD Components                                  }
{   (c) 2003 EzSoft Engineering                             }
{         All Rights Reserved                               }
{***********************************************************}

{$I EZ_FLAG.PAS}
Interface

Uses
  SysUtils, EzProjections, Math;

{ implemented projections  }

Procedure airy( P: TEzGeoConvert; init: boolean ); // Airy
Procedure merc( P: TEzGeoConvert; init: boolean ); // Mercator
Procedure tmerc( P: TEzGeoConvert; init: boolean ); // Transverse Mercator
Procedure utm( P: TEzGeoConvert; init: boolean ); // Universal Transverse Mercator (UTM)
Procedure omerc( P: TEzGeoConvert; init: boolean ); // Oblique mercator
Procedure aea( P: TEzGeoConvert; init: boolean ); // Albers Equal Area
Procedure leac( P: TEzGeoConvert; init: boolean ); // Lambert Equal Area Conic
Procedure lcc( P: TEzGeoConvert; init: boolean ); // Lambert Conformal Conic
Procedure bonne( P: TEzGeoConvert; init: boolean ); // Bonne (Werner lat_1=90)
Procedure cass( P: TEzGeoConvert; init: boolean ); // cassini
Procedure cc( P: TEzGeoConvert; init: boolean ); //  central cylindrical
Procedure cea( P: TEzGeoConvert; init: boolean ); // Equal Area Cylindrical
Procedure mill( P: TEzGeoConvert; init: boolean ); // Miller Cylindrical
Procedure moll( P: TEzGeoConvert; init: boolean ); // Mollweide
Procedure wag4( P: TEzGeoConvert; init: boolean ); // Wagner IV
Procedure wag5( P: TEzGeoConvert; init: boolean ); // Wagner V
Procedure eck4( P: TEzGeoConvert; init: boolean ); // Eckert IV
Procedure eck5( P: TEzGeoConvert; init: boolean ); // Eckert IV
Procedure eck6( P: TEzGeoConvert; init: boolean ); // Eckert IV
Procedure sinu( P: TEzGeoConvert; init: boolean );
Procedure gn_sinu( P: TEzGeoConvert; init: boolean );
Procedure mbtfps( P: TEzGeoConvert; init: boolean );
Procedure wag7( P: TEzGeoConvert; init: boolean );
Procedure tpeqd( P: TEzGeoConvert; init: boolean ); // Two Point Equidistant
Procedure tcea( P: TEzGeoConvert; init: boolean ); // Transverse Cylindrical Equal Area
Procedure vandg4( P: TEzGeoConvert; init: boolean ); // van der Grinten IV
Procedure laea( P: TEzGeoConvert; init: boolean ); // Lambert Azimuthal Equal Area
Procedure stere( P: TEzGeoConvert; init: boolean ); //Stereographic
Procedure ups( P: TEzGeoConvert; init: boolean ); //Universal Polar Stereographic
Procedure ortho( P: TEzGeoConvert; init: boolean ); // Orthographic
Procedure imw_p( P: TEzGeoConvert; init: boolean ); // International Map of the World Polyconic
Procedure poly( P: TEzGeoConvert; init: boolean ); // Polyconic (American)

{ general procedures }

//Procedure pj_enfn( P: TEzGeoConvert );
//Function pj_mlfn( phi, sphi, cphi: double; P: TEzGeoConvert ): double;
//Function pj_inv_mlfn( arg: double; P: TEzGeoConvert ): double;
//Function pj_tsfn( phi, sinphi, e: double ): double;
//Function pj_qsfn( sinphi, e, one_es: double ): double;
//Function pj_msfn( sinphi, cosphi, es: double ): double;
Function adjlon( Const lon: double ): double;
//Function pj_phi2( Const ts, e: double ): double;
//Function asqrt( Const v: double ): double;
// determine latitude from authalic latitude
//Procedure pj_authset( P: TEzGeoConvert );
//Function pj_authlat( Const beta: double; P: TEzGeoConvert ): double;

Implementation

{ general procedures }

Const
  ONE_TOL = 1.00000000000001;
  TOL =     0.000000001;
  ATOL =    1E-50;
  SPI =     3.14159265359;
  TWOPI =   6.2831853071795864769;
  HALFPI =  1.5707963267948966;
  TOL2 =    1.0E-10;
  N_ITER =  15;
  EPSILON = 1.0E-7;
  C00 =     1.0;
  C02 =     0.25;
  C04 =     0.046875;
  C06 =     0.01953125;
  C08 =     0.01068115234375;
  C22 =     0.75;
  C44 =     0.46875;
  C46 =     0.01302083333333333333;
  C48 =     0.00712076822916666666;
  C66 =     0.36458333333333333333;
  C68 =     0.00569661458333333333;
  C88 =     0.3076171875;
  EPS11 =   1E-11;
  MAX_ITER = 10;

Resourcestring
  serr0 = 'Inverse conversion not implemented';
  //serr1=	'no arguments in initialization list'; future use
  //serr2=	'no options found in ''init'' file"';
  //serr3=	'no colon in init= string';
  //serr4=	'projection not named';
  //serr5=	'unknown projection id';
  //serr6=	'effective eccentricity = 1.0';
  //serr7=	'unknown unit conversion id';
  //serr8=	'invalid boolean param argument';
  //serr9=	'unknown elliptical parameter name';
  //serr10=	'reciprocal flattening (1/f) = 0';
  //serr11=	'|radius reference latitude| > 90';
  //serr12=	'squared eccentricity < 0';
  //serr13=	'major axis or radius = 0 or not given';
  //serr14=	'latitude or longitude exceeded limits';
  //serr15=	'invalid x or y';
  //serr16=	'improperly formed DMS value';
  //serr17=	'non-convergent inverse meridinal dist';
  serr18 = 'non-convergent inverse phi2';
  serr19 = 'acos/asin: |arg| >1.+1e-14';
  //serr20=	'tolerance condition error';
  //serr21=	'conic lat_1 = -lat_2';
  //serr22=	'lat_1 >= 90';
  serr23 = 'lat_1 = 0';
  serr24 = 'lat_ts >= 90';
  serr25 = 'no distance between control points';
  //serr26=	'projection not selected to be rotated';
  //serr27=	'W <= 0 or M <= 0';
  //serr28=	'lsat not in 1-5 range';
  //serr29=	'path not in range';
  //serr30=	'h <= 0';
  //serr31=	'k <= 0';
  serr32 = 'lat_0 = 0 or 90 or alpha = 90';
  serr33 = 'lat_1=lat_2 or lat_1=0 or lat_2=90';
  serr34 = 'elliptical usage required';
  serr35 = 'invalid UTM zone number';
  serr36 = 'arg(s) out of range for Tcheby eval';
  //serr37=	'failed to find projection to be rotated';
  //serr38=	'failed to load NAD27-83 correction file';
  //serr39=	'both n & m must be spec''d and > 0';
  //serr40=	'n <= 0, n > 1 or not specified';
  //serr41=	'lat_1 or lat_2 not specified';
  //serr42=	'|lat_1| = |lat_2|';
  //serr43=	'lat_0 is pi/2 from mean lat';

  //projection dialog messages
  //SToXY = 'X = %.4f, Y = %.4f';
  //SToLatLong= 'Long = %.5n, Lat = %.5n';
  //SToLatLongCaption= 'Test X,Y to Lon,Lat';


Function adjlon( Const lon: double ): double;
Begin
  result := lon;
  While abs( result ) > SPI Do
    If result < 0 Then
      result := result + TWOPI
    Else
      result := result - TWOPI;
End;

// arc sin, cosine, tan2 and sqrt that will NOT fail
Function aasin( Const v: double ): double;
Var
  av: double;
Begin
  av := abs( v );
  If av >= 1.0 Then
  Begin
    If av > ONE_TOL Then
    Begin
      Raise Exception.Create( serr19 );
      //pj_errno := -19;
      exit;
    End;
    If v < 0 Then
      result := -HALFPI
    Else
      result := HALFPI;
    exit;
  End;
  result := arcsin( v );
End;

Function aacos( Const v: double ): double;
Var
  av: double;
Begin
  av := abs( v );
  If av >= 1.0 Then
  Begin
    If av > ONE_TOL Then
    Begin
      Raise Exception.Create( serr19 );
      //pj_errno := -19;
      exit;
    End;
    If v < 0 Then
      result := PI
    Else
      result := 0;
    exit;
  End;
  result := arccos( v );
End;

Function asqrt( Const v: double ): double;
Begin
  If v <= 0 Then
    result := 0
  Else
    result := sqrt( v );
End;

Function aatan2( Const n, d: double ): double;
Begin
  If ( abs( n ) < ATOL ) And ( abs( d ) < ATOL ) Then
    result := 0
  Else
    result := arctan2( n, d );
End;

Function pj_phi2( Const ts, e: double ): double;
Var
  eccnth, Phi, con, dphi: double;
  i: integer;
Begin
  eccnth := 0.5 * e;
  Phi := HALFPI - 2 * arctan( ts );
  i := N_ITER;
  Repeat
    con := e * sin( Phi );
    dphi := HALFPI - 2 * arctan( ts * power( ( 1 - con ) /
      ( 1 + con ), eccnth ) ) - Phi;
    Phi := Phi + dphi;
    Dec( i );
  Until ( abs( dphi ) <= TOL2 ) Or ( i <= 0 );
  If i <= 0 Then
    Raise Exception.Create( serr18 );
  //pj_errno := -18;
  result := Phi;
End;

Function pj_msfn( sinphi, cosphi, es: double ): double;
Begin
  result := ( cosphi / sqrt( 1 - es * sinphi * sinphi ) );
End;

//---------------------------

Function pj_qsfn( const sinphi, e, one_es: double ): double;
Var
  con: double;
Begin
  If e >= EPSILON Then
  Begin
    con := e * sinphi;
    result := ( one_es * ( sinphi / ( 1 - con * con ) -
      ( 0.5 / e ) * ln( ( 1 - con ) / ( 1 + con ) ) ) );
  End
  Else
    result := sinphi + sinphi;
End;

//---------------------------------------

Function pj_tsfn( const phi, sinphi, e: double ): double;
var
  sinphie: Double;
Begin
  sinphie := sinphi * e;
  result := ( tan( 0.5 * ( HALFPI - phi ) ) /
    power( ( 1.0 - sinphie ) / ( 1.0 + sinphie ), 0.5 * e ) );
End;

//-------------------------------------------------------------------//
(* meridiOnal distance for ellipsoid and inverse
**	8th degree - accurate to < 1e-5 meters when used in conjunction
**		with typical major axis values.
**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
*)

Procedure pj_enfn( P: TEzGeoConvert );
Var
  t, es: double;
Begin
  es := P.es;
  P.en[0] := C00 - es * ( C02 + es * ( C04 + es * ( C06 + es * C08 ) ) );
  P.en[1] := es * ( C22 - es * ( C04 + es * ( C06 + es * C08 ) ) );
  t := es * es;
  P.en[2] := t * ( C44 - es * ( C46 + es * C48 ) );
  t := t * es;
  P.en[3] := t * ( C66 - es * C68 );
  P.en[4] := t * es * C88;
End;

Function pj_mlfn( phi, sphi, cphi: double; P: TEzGeoConvert ): double;
Begin
  cphi := cphi * sphi;
  sphi := sphi * sphi;
  result := ( P.en[0] * phi - cphi * ( P.en[1] + sphi * ( P.en[2]
    + sphi * ( P.en[3] + sphi * P.en[4] ) ) ) );
End;

Function pj_inv_mlfn( arg: double; P: TEzGeoConvert ): double;
Var
  s, t, phi, k, es: double;
  i: integer;
Begin
  es := P.es;
  k := 1.0 / ( 1.0 - es );

  phi := arg;
  For i := 1 To MAX_ITER Do
  Begin (* rarely goes over 2 iterations *)
    s := sin( phi );
    t := 1.0 - es * s * s;
    t := ( pj_mlfn( phi, s, cos( phi ), P ) - arg ) * ( t * sqrt( t ) ) * k;
    phi := phi - t;
    If abs( t ) < EPS11 Then
    Begin
      result := phi;
      exit;
    End;
  End;
  //pj_errno := -17;
  result := phi;
End;

Const
  P00 = 0.33333333333333333333;
  P01 = 0.17222222222222222222;
  P02 = 0.10257936507936507936;
  P10 = 0.06388888888888888888;
  P11 = 0.06640211640211640211;
  P20 = 0.01641501294219154443;

Procedure pj_authset( P: TEzGeoConvert );
Var
  t: double;
  es: double;
Begin
  es := P.es;
  P.APA[0] := es * P00;
  t := es * es;
  P.APA[0] := P.APA[0] + ( t * P01 );
  P.APA[1] := t * P10;
  t := t * es;
  P.APA[0] := P.APA[0] + ( t * P02 );
  P.APA[1] := P.APA[1] + ( t * P11 );
  P.APA[2] := P.APA[2] + ( t * P20 );
End;

Function pj_authlat( Const beta: double; P: TEzGeoConvert ): double;
Var
  t: double;
Begin
  t := beta + beta;
  result := beta + P.APA[0] * sin( t ) + P.APA[1] * sin( t + t ) + P.APA[2] * sin( t + t + t );
End;

{ ***************** airy - Airy ***************** }

Const
  EPS = 1.E-10;
  N_POLE = 0;
  S_POLE = 1;
  EQUIT = 2;
  OBLIQ = 3;

Function s_forward_airy( lp: TLP; P: TEzGeoConvert ): TXY; (* spheroid *)
Var
  xy: TXY;
  sinlam, coslam, cosphi, sinphi, t, s, Krho, cosz: double;
Begin
  sinlam := sin( lp.lam );
  coslam := cos( lp.lam );
  Case P.mode Of
    EQUIT, OBLIQ:
      Begin
        sinphi := sin( lp.phi );
        cosphi := cos( lp.phi );
        cosz := cosphi * coslam;
        If P.mode = OBLIQ Then
          cosz := P.sinph0 * sinphi + P.cosph0 * cosz;
        If ( P.no_cut = false ) And ( cosz < -EPS ) Then
        Begin
          P.pj_errno := -20;
          result := xy;
          exit;
        End;
        s := 1 - cosz;
        If abs( s ) > EPS Then
        Begin
          t := 0.5 * ( 1 + cosz );
          Krho := -ln( t ) / s - P.Cb / t;
        End
        Else
          Krho := 0.5 - P.Cb;
        xy.x := Krho * cosphi * sinlam;
        If P.mode = OBLIQ Then
          xy.y := Krho * ( P.cosph0 * sinphi - P.sinph0 * cosphi * coslam )
        Else
          xy.y := Krho * sinphi;
      End;
    S_POLE, N_POLE:
      Begin
        lp.phi := abs( P.p_halfpi - lp.phi );
        If ( P.no_cut = false ) And ( ( lp.phi - EPS ) > HALFPI ) Then
        Begin
          P.pj_errno := -20;
          result := xy;
          exit;
        End;
        lp.phi := lp.phi * 0.5;
        If lp.phi > EPS Then
        Begin
          t := tan( lp.phi );
          Krho := -2 * ( ln( cos( lp.phi ) ) / t + t * P.Cb );
          xy.x := Krho * sinlam;
          xy.y := Krho * coslam;
          If P.mode = N_POLE Then
            xy.y := -xy.y
        End
        Else
        Begin
          xy.y := 0;
          xy.x := 0;
        End;
      End;
  End;
  result := xy;
End;

Procedure airy( P: TEzGeoConvert; init: boolean );
Var
  beta: double;
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;

  P.no_cut := P.pj_param.asboolean( 'no_cut' );
  beta := 0.5 * ( HALFPI - P.pj_param.asradians( 'lat_b' ) );
  If abs( beta ) < EPS Then
    P.Cb := -0.5
  Else
  Begin
    P.Cb := 1 / tan( beta );
    P.Cb := P.Cb * P.Cb * ln( cos( beta ) );
  End;
  If abs( abs( P.phi0 ) - HALFPI ) < EPS Then
  Begin
    If P.phi0 < 0 Then
    Begin
      P.p_halfpi := -HALFPI;
      P.mode := S_POLE;
    End
    Else
    Begin
      P.p_halfpi := HALFPI;
      P.mode := N_POLE;
    End;
  End
  Else
  Begin
    If abs( P.phi0 ) < EPS Then
      P.mode := EQUIT
    Else
    Begin
      P.mode := OBLIQ;
      P.sinph0 := sin( P.phi0 );
      P.cosph0 := cos( P.phi0 );
    End;
  End;
  P.fwd := @s_forward_airy;
  P.es := 0;

End;

{ ***************** merc - Mercator ***************** }

Const
  EPS10 = 1.E-10;

Function e_forward_merc( const lp: TLP; P: TEzGeoConvert ): TXY; (* ellipsoid *)
Var
  xy: TXY;
Begin
  If abs( abs( lp.phi ) - HALFPI ) <= EPS10 Then
  Begin
    P.pj_errno := -20;
    result := xy;
    exit;
  End;
  xy.x := P.k0 * lp.lam;
  xy.y := -P.k0 * ln( pj_tsfn( lp.phi, sin( lp.phi ), P.e ) );
  result := xy;
End;

Function s_forward_merc( const lp: TLP; P: TEzGeoConvert ): TXY; (* spheroid *)
Var
  xy: TXY;
Begin
  If abs( abs( lp.phi ) - HALFPI ) <= EPS10 Then
  Begin
    P.pj_errno := -20;
    result := xy;
    exit;
  End;
  xy.x := P.k0 * lp.lam;
  xy.y := P.k0 * ln( tan( FORTPI + 0.5 * lp.phi ) );
  result := xy;
End;

Function e_inverse_merc( const xy: TXY; P: TEzGeoConvert ): TLP; (* ellipsoid *)
Var
  lp: TLP;
Begin
  lp.phi := pj_phi2( exp( -xy.y / P.k0 ), P.e );
  If lp.phi = HUGE_VAL Then
  Begin
    P.pj_errno := -20;
    result := lp;
    exit;
  End;
  lp.lam := xy.x / P.k0;
  result := lp;
End;

Function s_inverse_merc( const xy: TXY; P: TEzGeoConvert ): TLP; (* spheroid *)
Begin
  result.phi := HALFPI - 2 * arctan( exp( -xy.y / P.k0 ) );
  result.lam := xy.x / P.k0;
End;