📄 ezprojimpl.pas

📁 很管用的GIS控件
💻 PAS
📖 第 1 页 / 共 5 页
字号:
Function e_inverse_aea( xy: TXY; P: TEzGeoConvert ): TLP; (* ellipsoid & spheroid *)
Var
  lp: TLP;
Begin
  xy.y := P.rho0 - xy.y;
  P.rho := hypot( xy.x, xy.y );
  If P.rho <> 0 Then
  Begin
    If P.n < 0 Then
    Begin
      P.rho := -P.rho;
      xy.x := -xy.x;
      xy.y := -xy.y;
    End;
    lp.phi := P.rho / P.dd;
    If P.ellips Then
    Begin
      lp.phi := ( P.c - lp.phi * lp.phi ) / P.n;
      If abs( P.ec - abs( lp.phi ) ) > TOL7 Then
      Begin
        lp.phi := phi1_( lp.phi, P.e, P.one_es );
        If lp.phi = HUGE_VAL Then
        Begin
          P.pj_errno := -20;
          result := lp;
          exit;
        End;
      End
      Else
      Begin
        If lp.phi < 0 Then
          lp.phi := -HALFPI
        Else
          lp.phi := HALFPI;
      End;
    End;
    lp.phi := ( P.c - lp.phi * lp.phi ) / P.n2;
    If abs( lp.phi ) <= 1 Then
    Begin
      lp.phi := arcsin( lp.phi );
    End
    Else
    Begin
      If lp.phi < 0 Then
        lp.phi := -HALFPI
      Else
        lp.phi := HALFPI;
    End;
    lp.lam := arctan2( xy.x, xy.y ) / P.n;
  End
  Else
  Begin
    lp.lam := 0;
    If P.n > 0 Then
      lp.phi := HALFPI
    Else
      lp.phi := -HALFPI;
  End;
  result := lp;
End;

Procedure setup_aea( P: TEzGeoConvert );
Var
  cosphi, sinphi: double;
  secant: boolean;
  ml1, m1: double;
  ml2, m2: double;
Begin
  If abs( P.phi1 + P.phi2 ) < EPS10 Then
  Begin
    P.pj_errno := -21;
    exit;
  End;
  sinphi := sin( P.phi1 );
  P.n := sinphi;
  cosphi := cos( P.phi1 );
  secant := abs( P.phi1 - P.phi2 ) >= EPS10;
  P.ellips := P.es > 0;
  If P.ellips Then
  Begin
    pj_enfn( P );
    m1 := pj_msfn( sinphi, cosphi, P.es );
    ml1 := pj_qsfn( sinphi, P.e, P.one_es );
    If secant Then
    Begin (* secant cone *)
      sinphi := sin( P.phi2 );
      cosphi := cos( P.phi2 );
      m2 := pj_msfn( sinphi, cosphi, P.es );
      ml2 := pj_qsfn( sinphi, P.e, P.one_es );
      P.n := ( m1 * m1 - m2 * m2 ) / ( ml2 - ml1 );
    End;
    P.ec := 1 - 0.5 * P.one_es * ln( ( 1 - P.e ) /
      ( 1 + P.e ) ) / P.e;
    P.c := m1 * m1 + P.n * ml1;
    P.dd := 1 / P.n;
    P.rho0 := P.dd * sqrt( P.c - P.n * pj_qsfn( sin( P.phi0 ),
      P.e, P.one_es ) );
  End
  Else
  Begin
    If secant Then
      P.n := 0.5 * ( P.n + sin( P.phi2 ) );
    P.n2 := P.n + P.n;
    P.c := cosphi * cosphi + P.n2 * sinphi;
    P.dd := 1 / P.n;
    P.rho0 := P.dd * sqrt( P.c - P.n2 * sin( P.phi0 ) );
  End;
  P.inv := @e_inverse_aea;
  P.fwd := @e_forward_aea;
End;

Procedure aea( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;

  P.phi1 := P.pj_param.AsRadians( 'lat_1' );
  P.phi2 := P.pj_param.AsRadians( 'lat_2' );

  setup_aea( P );
End;

Procedure leac( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;

  P.phi2 := P.pj_param.asradians( 'lat_1' );
  If P.pj_param.AsBoolean( 'south' ) Then
    P.phi1 := -HALFPI
  Else
    P.phi1 := HALFPI;

  setup_aea( P );
End;

{ lcc - Lambert Conformal Conic }

Function e_forward_lcc( lp: TLP; P: TEzGeoConvert ): TXY; (* ellipsoid & spheroid *)
Var
  xy: TXY;
  tmp: double;
Begin
  If abs( abs( lp.phi ) - HALFPI ) < EPS10 Then
  Begin
    If ( lp.phi * P.n ) <= 0 Then
    Begin
      P.pj_errno := -20;
      result := xy;
      exit;
    End;
    P.rho := 0;
  End
  Else
  Begin
    If P.ellips Then
      tmp := power( pj_tsfn( lp.phi, sin( lp.phi ), P.e ), P.n )
    Else
      tmp := power( tan( FORTPI + 0.5 * lp.phi ), -P.n );
    P.rho := P.c * tmp;
  End;
  lp.lam := lp.lam * P.n;
  xy.x := P.k0 * ( P.rho * sin( lp.lam ) );
  xy.y := P.k0 * ( P.rho0 - P.rho * cos( lp.lam ) );
  result := xy;
End;

Function e_inverse_lcc( xy: TXY; P: TEzGeoConvert ): TLP; (* ellipsoid & spheroid *)
Var
  lp: TLP;
Begin
  xy.x := xy.x / P.k0;
  xy.y := xy.y / P.k0;
  xy.y := P.rho0 - xy.y;
  P.rho := hypot( xy.x, xy.y );
  If P.rho <> 0 Then
  Begin
    If P.n < 0 Then
    Begin
      P.rho := -P.rho;
      xy.x := -xy.x;
      xy.y := -xy.y;
    End;
    If P.ellips Then
    Begin
      lp.phi := pj_phi2( power( P.rho / P.c, 1 / P.n ), P.e );
      If lp.phi = HUGE_VAL Then
      Begin
        P.pj_errno := -20;
        result := lp;
        exit;
      End;
    End
    Else
      lp.phi := 2 * arctan( power( P.c / P.rho, 1 / P.n ) ) - HALFPI;
    lp.lam := arctan2( xy.x, xy.y ) / P.n;
  End
  Else
  Begin
    lp.lam := 0;
    If P.n > 0 Then
      lp.phi := HALFPI
    Else
      lp.phi := -HALFPI;
  End;
  result := lp;
End;

Procedure factors( const lp: TLP; P: TEzGeoConvert; Var fac: TFACTORS );
Var
  tmp: double;
Begin
  If abs( abs( lp.phi ) - HALFPI ) < EPS10 Then
  Begin
    If ( lp.phi * P.n ) <= 0 Then
      exit;
    P.rho := 0;
  End
  Else
  Begin
    If P.ellips Then
      tmp := power( pj_tsfn( lp.phi, sin( lp.phi ), P.e ), P.n )
    Else
      tmp := power( tan( FORTPI + 0.5 * lp.phi ), -P.n );
    P.rho := P.c * tmp;
  End;
  fac.code := fac.code Or ( IS_ANAL_HK + IS_ANAL_CONV );
  fac.h := P.k0 * P.n * P.rho / pj_msfn( sin( lp.phi ), cos( lp.phi ), P.es );
  fac.k := fac.h;
  fac.conv := -P.n * lp.lam;
End;

Procedure lcc( P: TEzGeoConvert; init: boolean );
Var
  cosphi, sinphi: double;
  secant: boolean;
  ml1, m1, tmp: double;
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;

  P.phi1 := P.pj_param.asradians( 'lat_1' );
  If P.pj_param.defined( 'lat_2' ) Then
    P.phi2 := P.pj_param.asradians( 'lat_2' )
  Else
  Begin
    P.phi2 := P.phi1;
    If Not P.pj_param.defined( 'lat_0' ) Then
      P.phi0 := P.phi1;
  End;
  If abs( P.phi1 + P.phi2 ) < EPS10 Then
  Begin
    P.pj_errno := -21;
    exit;
  End;
  sinphi := sin( P.phi1 );
  P.n := sinphi;
  cosphi := cos( P.phi1 );
  secant := abs( P.phi1 - P.phi2 ) >= EPS10;
  P.ellips := P.es <> 0;
  If P.ellips Then
  Begin
    P.e := sqrt( P.es );
    m1 := pj_msfn( sinphi, cosphi, P.es );
    ml1 := pj_tsfn( P.phi1, sinphi, P.e );
    If secant Then
    Begin (* secant cone *)
      sinphi := sin( P.phi2 );
      P.n := ln( m1 / pj_msfn( sinphi, cos( P.phi2 ), P.es ) );
      P.n := P.n / ln( ml1 / pj_tsfn( P.phi2, sinphi, P.e ) );
    End;
    P.rho0 := m1 * power( ml1, -P.n ) / P.n;
    P.c := P.rho0;
    If abs( abs( P.phi0 ) - HALFPI ) < EPS10 Then
      tmp := 0
    Else
      tmp := power( pj_tsfn( P.phi0, sin( P.phi0 ), P.e ), P.n );
    P.rho0 := P.rho0 * tmp;
  End
  Else
  Begin
    If secant Then
      P.n := ln( cosphi / cos( P.phi2 ) ) /
        ln( tan( FORTPI + 0.5 * P.phi2 ) /
        tan( FORTPI + 0.5 * P.phi1 ) );
    P.c := cosphi * power( tan( FORTPI + 0.5 * P.phi1 ), P.n ) / P.n;
    If abs( abs( P.phi0 ) - HALFPI ) < EPS10 Then
      tmp := 0
    Else
      tmp := P.c * power( tan( FORTPI + 0.5 * P.phi0 ), -P.n );
    P.rho0 := tmp;
  End;
  P.inv := @e_inverse_lcc;
  P.fwd := @e_forward_lcc;
  P.spc := @factors;

End;

{ bonne - Bonne (Werner lat_1=90)}

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

Function s_forward_bonne( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  E, rh: double;
  xy: TXY;
Begin
  rh := P.cphi1 + P.phi1 - lp.phi;
  If ( abs( rh ) > EPS10 ) Then
  Begin
    E := lp.lam * cos( lp.phi ) / rh;
    xy.x := rh * sin( E );
    xy.y := P.cphi1 - rh * cos( E );
  End
  Else
  Begin
    xy.y := 0.0;
    xy.x := 0.0;
  End;
  result := xy;
End;

Function s_inverse_bonne( xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  rh: double;
  lp: TLP;
Begin
  xy.y := P.cphi1 - xy.y;
  rh := hypot( xy.x, xy.y );
  lp.phi := P.cphi1 + P.phi1 - rh;
  If abs( lp.phi ) > HALFPI Then
    Raise Exception.Create( 'Error' );
  If abs( abs( lp.phi ) - HALFPI ) <= EPS10 Then
    lp.lam := 0.0
  Else
    lp.lam := rh * aatan2( xy.x, xy.y ) / cos( lp.phi );
  result := lp;
End;

Function e_inverse_bonne( xy: TXY; P: TEzGeoConvert ): TLP; // ellipsoid
Var
  s, rh: double;
  lp: TLP;
Begin
  xy.y := P.am1 - xy.y;
  rh := hypot( xy.x, xy.y );
  lp.phi := pj_inv_mlfn( P.am1 + P.m1 - rh, P );
  s := abs( lp.phi );
  If s < HALFPI Then
  Begin
    s := sin( lp.phi );
    lp.lam := rh * aatan2( xy.x, xy.y ) * sqrt( 1.0 - P.es * s * s ) / cos( lp.phi );
  End
  Else If abs( s - HALFPI ) <= EPS10 Then
    lp.lam := 0.0
  Else
    Raise Exception.Create( 'Error' );
  result := lp;
End;

Procedure bonne( P: TEzGeoConvert; init: boolean );
Var
  c: double;
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.phi1 := P.pj_param.asradians( 'lat_1' );
  If abs( P.phi1 ) < EPS10 Then
    Raise Exception.Create( SErr23 );
  If P.es <> 0 Then
  Begin
    pj_enfn( P );
    P.am1 := sin( P.phi1 );
    c := cos( P.phi1 );
    P.m1 := pj_mlfn( P.phi1, P.am1, c, P );
    P.am1 := c / ( sqrt( 1.0 - P.es * P.am1 * P.am1 ) * P.am1 );
    P.inv := @e_inverse_bonne;
    P.fwd := @e_forward_bonne;
  End
  Else
  Begin
    If abs( P.phi1 ) + EPS10 >= HALFPI Then
      P.cphi1 := 0.0
    Else
      P.cphi1 := 1.0 / tan( P.phi1 );
    P.inv := @s_inverse_bonne;
    P.fwd := @s_forward_bonne;
  End
End;

{ cassini }
Const
  C1 = 0.16666666666666666666;
  C2 = 0.00833333333333333333;
  C3 = 0.04166666666666666666;
  C4 = 0.33333333333333333333;
  C5 = 0.06666666666666666666;

Function e_forward_cass( const lp: TLP; P: TEzGeoConvert ): TXY; // ellipsoid
Var
  xy: TXY;
Begin
  P.n := sin( lp.phi );
  P.c := cos( lp.phi );
  xy.y := pj_mlfn( lp.phi, P.n, P.c, P );
  P.n := 1.0 / sqrt( 1.0 - P.es * P.n * P.n );
  P.tn := tan( lp.phi );
  P.t := P.tn * P.tn;
  P.a1 := lp.lam * P.c;
  P.c := P.c * ( P.es * P.c / ( 1.0 - P.es ) );
  P.a2 := P.a1 * P.a1;
  xy.x := P.n * P.a1 * ( 1.0 - P.a2 * P.t * ( C1 - ( 8.0 - P.t + 8.0 * P.c ) * P.a2 * C2 ) );
  xy.y := xy.y - ( P.m0 - P.n * P.tn * P.a2 * ( 0.5 + ( 5.0 - P.t + 6.0 * P.c ) * P.a2 * C3 ) );
  result := xy;
End;

Function s_forward_cass( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid */
Var
  xy: TXY;
Begin
  xy.x := aasin( cos( lp.phi ) * sin( lp.lam ) );
  xy.y := aatan2( tan( lp.phi ), cos( lp.lam ) ) - P.phi0;
  result := xy;
End;

Function e_inverse_cass( const xy: TXY; P: TEzGeoConvert ): TLP; // ellipsoid
Var
  ph1: double;
  lp: TLP;
Begin
  ph1 := pj_inv_mlfn( P.m0 + xy.y, P );
  P.tn := tan( ph1 );
  P.t := P.tn * P.tn;
  P.n := sin( ph1 );
  P.r := 1.0 / ( 1.0 - P.es * P.n * P.n );
  P.n := sqrt( P.r );
  P.r := P.r * ( ( 1.0 - P.es ) * P.n );
  P.dd := xy.x / P.n;
  P.d2 := P.dd * P.dd;
  lp.phi := ph1 - ( P.n * P.tn / P.r ) * P.d2 * ( 0.5 - ( 1.0 + 3.0 * P.t ) * P.d2 * C3 );
  lp.lam := P.dd * ( 1.0 + P.t * P.d2 * ( -C4 + ( 1.0 + 3.0 * P.t ) * P.d2 * C5 ) ) / cos( ph1 );
  result := lp;
End;

Function s_inverse_cass( const xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  lp: TLP;
Begin
  P.dd := xy.y + P.phi0;
  lp.phi := aasin( sin( P.dd ) * cos( xy.x ) );
  lp.lam := aatan2( tan( xy.x ), cos( P.dd ) );
  result := lp;
End;

Procedure cass( P: TEzGeoConvert; init: boolean );
Begin
  If P.es <> 0 Then
  Begin
    pj_enfn( P );
    //if (!(P.en = pj_enfn(P.es))) then raise exception.create('Error 0');
    P.m0 := pj_mlfn( P.phi0, sin( P.phi0 ), cos( P.phi0 ), P );
    P.inv := @e_inverse_cass;
    P.fwd := @e_forward_cass;
  End
  Else
  Begin
    P.inv := @s_inverse_cass;
    P.fwd := @s_forward_cass;
  End;
End;

{ cc - central cylindrical }

Function s_forward_cc( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
Begin
  If abs( abs( lp.phi ) - HALFPI ) <= EPS10 Then
    Raise exception.create( 'error' );
  xy.x := lp.lam;
  xy.y := tan( lp.phi );
  result := xy;
End;

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

Procedure cc( 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_cc;
  P.fwd := @s_forward_cc;
💿 文件大小 1389 K
👤 上传用户 dasdasdc
📂 所属分类人工智能/神经网络
🏷️ 相关标签

#GIS #控件
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -