ezprojimpl.pas

很管用的GIS控件

End;

{ cea - Equal Area Cylindrical }

Function e_forward_cea( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
Begin
  xy.x := P.k0 * lp.lam;
  xy.y := 0.5 * pj_qsfn( sin( lp.phi ), P.e, P.one_es ) / P.k0;
  result := xy;
End;

Function s_forward_cea( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
Begin
  xy.x := P.k0 * lp.lam;
  xy.y := sin( lp.phi ) / P.k0;
  result := xy;
End;

Function e_inverse_cea( const xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  lp: TLP;
Begin
  lp.phi := pj_authlat( arcsin( 2.0 * xy.y * P.k0 / P.qp ), P );
  lp.lam := xy.x / P.k0;
  result := lp;
End;

Function s_inverse_cea( xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  t: double;
  lp: TLP;
Begin
  xy.y := xy.y * P.k0;
  t := abs( xy.y );
  If t - EPS <= 1.0 Then
  Begin
    If ( t >= 1.0 ) Then
    Begin
      If xy.y < 0 Then
        lp.phi := -HALFPI
      Else
        lp.phi := HALFPI;
    End
    Else
      lp.phi := arcsin( xy.y );
    lp.lam := xy.x / P.k0;
  End
  Else
    Raise exception.create( 'error' );
  result := lp;
End;

Procedure cea( P: TEzGeoConvert; init: boolean );
Var
  t: double;
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  If P.pj_param.defined( 'lat_ts' ) Then
  Begin
    t := P.pj_param.asradians( 'lat_ts' );
    P.k0 := cos( t );
    If P.k0 < 0 Then
      Raise exception.create( serr24 )
  End
  Else
    t := 0;
  If P.es <> 0 Then
  Begin
    t := sin( t );
    P.k0 := P.k0 / sqrt( 1.0 - P.es * t * t );
    P.e := sqrt( P.es );
    pj_authset( P );
    //if (!(P->apa = pj_authset(P->es))) E_ERROR_0;
    P.qp := pj_qsfn( 1., P.e, P.one_es );
    P.inv := @e_inverse_cea;
    P.fwd := @e_forward_cea;
  End
  Else
  Begin
    P.inv := @s_inverse_cea;
    P.fwd := @s_forward_cea;
  End;
End;

{ mill - miller cylindrical }

Function s_forward_mill( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
Begin
  xy.x := lp.lam;
  xy.y := ln( tan( FORTPI + lp.phi * 0.4 ) ) * 1.25;
  result := xy;
End;

Function s_inverse_mill( const xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  lp: TLP;
Begin
  lp.lam := xy.x;
  lp.phi := 2.5 * ( arctan( exp( 0.8 * xy.y ) ) - FORTPI );
  result := lp;
End;

Procedure mill( P: TEzGeoConvert; init: boolean ); //Miller Cylindrical
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.es := 0.0;
  P.inv := @s_inverse_mill;
  P.fwd := @s_forward_mill;
End;

{ Mollweide, Wagner IV, Wagner V}
Const
  LOOP_TOL = 1E-7;

Function s_forward_moll( lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
  k, V: double;
  i: integer;
Begin
  k := P.C_p * sin( lp.phi );
  i := MAX_ITER;
  While i > 0 Do
  Begin
    V := ( lp.phi + sin( lp.phi ) - k ) / ( 1.0 + cos( lp.phi ) );
    lp.phi := lp.phi - V;
    If abs( V ) < LOOP_TOL Then
      break;
    dec( i );
  End;
  If i = 0 Then
  Begin
    If lp.phi < 0 Then
      lp.phi := -HALFPI
    Else
      lp.phi := HALFPI;
  End
  Else
    lp.phi := lp.phi * 0.5;
  xy.x := P.C_x * lp.lam * cos( lp.phi );
  xy.y := P.C_y * sin( lp.phi );
  result := xy;
End;

Function s_inverse_moll( const xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  lp: TLP;
  //th, s:double;
Begin
  lp.phi := aasin( xy.y / P.C_y );
  lp.lam := xy.x / ( P.C_x * cos( lp.phi ) );
  lp.phi := lp.phi + lp.phi;
  lp.phi := aasin( ( lp.phi + sin( lp.phi ) ) / P.C_p );
  result := lp;
End;

Procedure setup_moll( P: TEzGeoConvert; Const dp: double );
Var
  r, sp, p2: double;
Begin
  p2 := dp + dp;
  P.es := 0;
  sp := sin( dp );
  r := sqrt( TWOPI * sp / ( p2 + sin( p2 ) ) );
  P.C_x := 2.0 * r / PI;
  P.C_y := r / sp;
  P.C_p := p2 + sin( p2 );
  P.inv := @s_inverse_moll;
  P.fwd := @s_forward_moll;
End;

Procedure moll( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  setup_moll( P, HALFPI );
End;

Procedure wag4( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  setup_moll( P, PI / 3 );
End;

Procedure wag5( P: TEzGeoConvert; init: boolean ); //Miller Cylindrical
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.es := 0;
  P.C_x := 0.90977;
  P.C_y := 1.65014;
  P.C_p := 3.00896;
  P.inv := @s_inverse_moll;
  P.fwd := @s_forward_moll;
End;

{ eck4 - Eckert IV }
Const
  C_x = 0.42223820031577120149;
  C_y = 1.32650042817700232218;
  RC_y = 0.75386330736002178205;
  C_p = 3.57079632679489661922;
  RC_p = 0.28004957675577868795;
  EPS7 = 1E-7;
  NITER = 6;

Function s_forward_eck4( lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
  dp, V, s, c: double;
  i: integer;
Begin
  dp := C_p * sin( lp.phi );
  V := lp.phi * lp.phi;
  lp.phi := lp.phi * ( 0.895168 + V * ( 0.0218849 + V * 0.00826809 ) );
  i := NITER;
  While i > 0 Do
  Begin
    c := cos( lp.phi );
    s := sin( lp.phi );
    V := ( lp.phi + s * ( c + 2.0 ) - dp ) / ( 1.0 + c * ( c + 2.0 ) - s * s );
    lp.phi := lp.phi - V;
    If abs( V ) < EPS7 Then
      break;
    dec( i );
  End;
  If i = 0 Then
  Begin
    xy.x := C_x * lp.lam;
    If lp.phi < 0 Then
      xy.y := -C_y
    Else
      xy.y := C_y;
  End
  Else
  Begin
    xy.x := C_x * lp.lam * ( 1.0 + cos( lp.phi ) );
    xy.y := C_y * sin( lp.phi );
  End;
  result := xy;
End;

Function s_inverse_eck4( const xy: TXY; P: TEzGeoConvert ): TLP; // spheroid */
Var
  lp: TLP;
  c: double;
Begin
  lp.phi := aasin( xy.y / C_y );
  c := cos( lp.phi );
  lp.lam := xy.x / ( C_x * ( 1.0 + c ) );
  lp.phi := aasin( ( lp.phi + sin( lp.phi ) * ( c + 2.0 ) ) / C_p );
  result := lp;
End;

Procedure eck4( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.es := 0.0;
  P.inv := @s_inverse_eck4;
  P.fwd := @s_forward_eck4;
End;

{ eck5 - Eckert V}
Const
  XF = 0.44101277172455148219;
  RXF = 2.26750802723822639137;
  YF = 0.88202554344910296438;
  RYF = 1.13375401361911319568;

Function s_forward( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
Begin
  xy.x := XF * ( 1.0 + cos( lp.phi ) ) * lp.lam;
  xy.y := YF * lp.phi;
  result := xy;
End;

Function s_inverse( const xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  lP: TLP;
Begin
  lp.phi := RYF * xy.y;
  lp.lam := RXF * xy.x / ( 1.0 + cos( lp.phi ) );
  result := lp;
End;

Procedure eck5( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.es := 0.0;
  P.inv := @s_inverse;
  P.fwd := @s_forward;
End;

{ eck6 - Eckert VI }
Const
  MAX_ITER8 = 8;

  // Ellipsoidal Sinusoidal only

Function e_forward_eck6( const lp: TLP; P: TEzGeoConvert ): TXY; // ellipsoid
Var
  s, c: double;
  xy: TXY;
Begin
  s := sin( lp.phi );
  c := cos( lp.phi );
  xy.y := pj_mlfn( lp.phi, s, c, P );
  xy.x := lp.lam * c / sqrt( 1.0 - P.es * s * s );
  result := xy;
End;

Function e_inverse_eck6( const xy: TXY; P: TEzGeoConvert ): TLP; // ellipsoid
Var
  lp: TLP;
  s: double;
Begin
  lp.phi := pj_inv_mlfn( xy.y, P );
  s := abs( lp.phi );
  If s < HALFPI Then
  Begin
    s := sin( lp.phi );
    lp.lam := xy.x * sqrt( 1.0 - P.es * s * s ) / cos( lp.phi );
  End
  Else If ( ( s - EPS10 ) < HALFPI ) Then
    lp.lam := 0.0
  Else
    Raise Exception.create( 'error' );
  result := lp;
End;

// General spherical sinusoidals
Const
  MAXITER8 = 8;

Function s_forward_eck6( lp: TLP; P: TEzGeoConvert ): TXY; // sphere
Var
  xy: TXY;
  k, V: double;
  i: integer;
Begin
  If P.m = 0 Then
  Begin
    If P.n <> 1.0 Then
      lp.phi := aasin( P.n * sin( lp.phi ) );
  End
  Else
  Begin
    k := P.n * sin( lp.phi );
    i := MAXITER8;
    While i > 0 Do
    Begin
      V := ( P.m * lp.phi + sin( lp.phi ) - k ) / ( P.m + cos( lp.phi ) );
      lp.phi := lp.phi - V;
      If abs( V ) < LOOP_TOL Then
        break;
      dec( i );
    End;
    If i = 0 Then
      Raise Exception.create( 'error' );
  End;
  xy.x := P.C_x * lp.lam * ( P.m + cos( lp.phi ) );
  xy.y := P.C_y * lp.phi;
  result := xy;
End;

Function s_inverse_eck6( xy: TXY; P: TEzGeoConvert ): TLP; // ellipsoid
Var
  //s:double;
  lp: TLP;
Begin
  xy.y := xy.y / P.C_y;
  If P.m <> 0 Then
    lp.phi := aasin( ( P.m * xy.y + sin( xy.y ) ) / P.n )
  Else
  Begin
    If P.n <> 1.0 Then
      lp.phi := aasin( sin( xy.y ) / P.n )
    Else
      lp.phi := xy.y;
  End;
  lp.lam := xy.x / ( P.C_x * ( P.m + cos( xy.y ) ) );
  result := lp;
End;

// for spheres, only

Procedure setup_eck6( P: TEzGeoConvert );
Begin
  P.es := 0;
  P.C_y := sqrt( ( P.m + 1.0 ) / P.n );
  P.C_x := P.C_y / ( P.m + 1.0 );
  P.inv := @s_inverse_eck6;
  P.fwd := @s_forward_eck6;
End;

Procedure sinu( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  pj_enfn( P );
  If P.es <> 0 Then
  Begin
    pj_enfn( P );
    P.inv := @e_inverse_eck6;
    P.fwd := @e_forward_eck6;
  End
  Else
  Begin
    P.n := 1.0;
    P.m := 0.0;
    setup_eck6( P );
  End;
End;

Procedure eck6( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.m := 1.0;
  P.n := 2.570796326794896619231321691;
  setup_eck6( P );
End;

Procedure mbtfps( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.m := 0.5;
  P.n := 1.785398163397448309615660845;
  setup_eck6( P );
End;

Procedure gn_sinu( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  If P.pj_param.defined( 'n' ) And P.pj_param.defined( 'm' ) Then
  Begin
    P.n := P.pj_param.asfloat( 'n' );
    P.m := P.pj_param.asfloat( 'm' );
  End
  Else
    exception.create( 'error 99' );
  setup_eck6( P );
End;

{ wag7 - Wagner VII }

Function s_forward_wag7( lp: TLP; P: TEzGeoConvert ): TXY; // sphere
Var
  xy: TXY;
  theta, ct, D: double;
Begin
  xy.y := 0.90630778703664996 * sin( lp.phi );
  theta := arcsin( xy.y );
  ct := cos( theta );
  lp.lam := lp.lam / 3.0;
  xy.x := 2.66723 * ct * sin( lp.lam );
  D := 1 / ( sqrt( 0.5 * ( 1 + ct * cos( lp.lam ) ) ) );
  xy.y := xy.y * ( 1.24104 * D );
  xy.x := xy.x * D;
  result := xy;
End;

Procedure wag7( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.fwd := @s_forward_wag7;
  P.inv := Nil;
  P.es := 0.0;