ezprojimpl.pas

很管用的GIS控件

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

  is_phits := P.pj_param.defined( 'lat_ts' );
  phits := 0;
  If is_phits Then
  Begin
    phits := abs( P.pj_param.asradians( 'lat_ts' ) );
    If phits >= HALFPI Then
    Begin
      P.pj_errno := -24;
      exit;
    End;
  End;
  If P.es <> 0 Then
  Begin (* ellipsoid *)
    If is_phits Then
      P.k0 := pj_msfn( sin( phits ), cos( phits ), P.es );
    P.inv := @e_inverse_merc;
    P.fwd := @e_forward_merc;
  End
  Else
  Begin (* sphere *)
    If is_phits Then
      P.k0 := cos( phits );
    P.inv := @s_inverse_merc;
    P.fwd := @s_forward_merc;
  End;
End;

{ ***************** tmerc - Traverse Mercator ***************** }

Const
  FC1 = 1.0;
  FC2 = 0.5;
  FC3 = 0.16666666666666666666;
  FC4 = 0.08333333333333333333;
  FC5 = 0.05;
  FC6 = 0.03333333333333333333;
  FC7 = 0.02380952380952380952;
  FC8 = 0.01785714285714285714;

Function e_forward_tmerc( const lp: TLP; P: TEzGeoConvert ): TXY; // ellipse
Var
  xy: TXY;
  al, als, n, cosphi, sinphi, t: double;
Begin
  sinphi := sin( lp.phi );
  cosphi := cos( lp.phi );
  If abs( cosphi ) > 1E-10 Then
    t := sinphi / cosphi
  Else
    t := 0;
  t := t * t;
  al := cosphi * lp.lam;
  als := al * al;
  al := al / sqrt( 1.0 - P.es * sinphi * sinphi );
  n := P.esp * cosphi * cosphi;
  xy.x := P.k0 * al * ( FC1 +
    FC3 * als * ( 1.0 - t + n +
    FC5 * als * ( 5.0 + t * ( t - 18 ) + n * ( 14.0 - 58.0 * t )
    + FC7 * als * ( 61.0 + t * ( t * ( 179.0 - t ) - 479.0 ) ) ) ) );
  xy.y := P.k0 * ( pj_mlfn( lp.phi, sinphi, cosphi, P ) - P.ml0 +
    sinphi * al * lp.lam * FC2 * ( 1.0 +
    FC4 * als * ( 5.0 - t + n * ( 9.0 + 4.0 * n ) +
    FC6 * als * ( 61.0 + t * ( t - 58.0 ) + n * ( 270.0 - 330 * t ) + FC8 * als * ( 1385.0 + t * ( t * ( 543.0 - t ) - 3111.0 )
      ) ) ) ) );
  result := xy;
End;

Function s_forward_tmerc( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
  b, cosphi: double;
Begin
  cosphi := cos( lp.phi );
  b := cosphi * sin( lp.lam );
  If abs( abs( b ) - 1.0 ) <= EPS10 Then
  Begin
    P.pj_errno := -20;
    result := xy;
    exit;
  End;
  xy.x := P.ml0 * ln( ( 1.0 + b ) / ( 1.0 - b ) ); //falta checar funcion log
  xy.y := cosphi * cos( lp.lam ) / sqrt( 1.0 - b * b );
  b := abs( xy.y );
  If b >= 1.0 Then
  Begin
    If ( b - 1.0 ) > EPS10 Then
    Begin
      P.pj_errno := -20;
      result := xy;
      exit;
    End
    Else
      xy.y := 0.0;
  End
  Else
    xy.y := arccos( xy.y );
  If lp.phi < 0 Then
    xy.y := -xy.y;
  xy.y := P.esp * ( xy.y - P.phi0 );
  result := xy;
End;

Function e_inverse_tmerc( const xy: TXY; P: TEzGeoConvert ): TLP; // ellipsoid
Var
  lp: TLP;
  n, con, cosphi, d, ds, sinphi, t: double;
Begin
  lp.phi := pj_inv_mlfn( P.ml0 + xy.y / P.k0, P );
  If abs( lp.phi ) >= HALFPI Then
  Begin
    If xy.y < 0.0 Then
      lp.phi := -HALFPI
    Else
      lp.phi := HALFPI;
    lp.lam := 0.0;
  End
  Else
  Begin
    sinphi := sin( lp.phi );
    cosphi := cos( lp.phi );
    If abs( cosphi ) > 1E-10 Then
      t := sinphi / cosphi
    Else
      t := 0;
    n := P.esp * cosphi * cosphi;
    con := 1.0 - P.es * sinphi * sinphi;
    d := xy.x * sqrt( con ) / P.k0;
    con := con * t;
    t := t * t;
    ds := d * d;
    lp.phi := lp.phi - ( con * ds / ( 1.0 - P.es ) ) * FC2 * ( 1.0 -
      ds * FC4 * ( 5.0 + t * ( 3.0 - 9.0 * n ) + n * ( 1.0 - 4 * n ) -
      ds * FC6 * ( 61.0 + t * ( 90.0 - 252.0 * n +
      45.0 * t ) + 46.0 * n - ds * FC8 * ( 1385.0 + t * ( 3633.0 + t * ( 4095.0 + 1574.0 * t ) ) ) ) ) );
    lp.lam := d * ( FC1 - ds * FC3 * ( 1.0 + 2. * t + n - ds * FC5 * ( 5.0 + t * ( 28.0 + 24.0 * t + 8. * n ) + 6.0 * n
      - ds * FC7 * ( 61.0 + t * ( 662.0 + t * ( 1320.0 + 720.0 * t ) ) ) ) ) ) / cosphi;
  End;
  result := lp;
End;

Function s_inverse_tmerc( const xy: TXY; P: TEzGeoConvert ): TLP; // sphere
Var
  lp: TLP;
  h, g: double;
Begin
  h := exp( xy.x / P.esp );
  g := 0.5 * ( h - 1.0 / h );
  h := cos( P.phi0 + xy.y / P.esp );
  lp.phi := arcsin( sqrt( ( 1.0 - h * h ) / ( 1.0 + g * g ) ) );
  If xy.y < 0 Then
    lp.phi := -lp.phi;
  If ( g <> 0 ) And ( h <> 0 ) Then
    lp.lam := arctan2( g, h )
  Else
    lp.lam := 0;
  result := lp;
End;

Procedure setup_tmerc( P: TEzGeoConvert );
Begin
  If P.es <> 0 Then
  Begin
    pj_enfn( P );
    P.ml0 := pj_mlfn( P.phi0, sin( P.phi0 ), cos( P.phi0 ), P );
    P.esp := P.es / ( 1.0 - P.es );
    P.inv := @e_inverse_tmerc;
    P.fwd := @e_forward_tmerc;
  End
  Else
  Begin
    P.esp := P.k0;
    P.ml0 := 0.5 * P.esp;
    P.inv := @s_inverse_tmerc;
    P.fwd := @s_forward_tmerc;
  End;
End;

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

  setup_tmerc( P );
End;

Procedure utm( P: TEzGeoConvert; init: boolean );
Var
  zone: integer;
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;

  If p.es = 0 Then
    Raise Exception.Create( SErr34 );
  If P.pj_param.asboolean( 'south' ) Then
    p.y0 := 10000000.0
  Else
    p.y0 := 0;
  p.x0 := 500000.;
  If P.pj_param.defined( 'zone' ) Then
  Begin //* zone input ? */
    zone := P.pj_param.asinteger( 'zone' );
    If ( zone > 0 ) And ( zone <= 60 ) Then
      Dec( zone )
    Else
      Raise Exception.Create( SErr35 );
  End
  Else
  Begin //* nearest central meridian input */
    zone := floor( ( adjlon( P.lam0 ) + PI ) * 30. / PI );
    If zone < 0 Then
      zone := 0
    Else If zone >= 60 Then
      zone := 59;
  End;
  p.lam0 := ( zone + 0.5 ) * PI / 30.0 - PI;
  p.k0 := 0.9996;
  p.phi0 := 0.0;

  setup_tmerc( P );
End;

{ ***************** omerc - oblique mercator ***************** }

Const
  TOL7 = 1.0E-7;

Function TSFN0( Const x: double ): double;
Begin
  result := tan( 0.5 * ( HALFPI - x ) );
End;

Function omerc_forward( const lp: TLP; P: TEzGeoConvert ): TXY; // ellipsoid & spheroid */
Var
  con, q, s, ul, us, vl, vs, temp: double;
  xy: TXY;
Begin
  vl := sin( P.bl * lp.lam );
  If ( abs( abs( lp.phi ) - HALFPI ) <= EPS ) Then
  Begin
    If lp.phi < 0 Then
      ul := -P.singam
    Else
      ul := P.singam;
    us := P.al * lp.phi / P.bl;
  End
  Else
  Begin
    If P.ellips Then
      temp := power( pj_tsfn( lp.phi, sin( lp.phi ), P.e ), P.bl )
    Else
      temp := TSFN0( lp.phi );
    q := P.el / temp;
    s := 0.5 * ( q - 1 / q );
    ul := 2 * ( s * P.singam - vl * P.cosgam ) / ( q + 1 / q );
    con := cos( P.bl * lp.lam );
    If ( abs( con ) >= TOL7 ) Then
    Begin
      us := P.al * arctan( ( s * P.cosgam + vl * P.singam ) / con ) / P.bl;
      If ( con < 0 ) Then
        us := us + PI * P.al / P.bl;
    End
    Else
      us := P.al * P.bl * lp.lam;
  End;
  If ( abs( abs( ul ) - 1 ) <= EPS ) Then
    exit; //F_ERROR;
  vs := 0.5 * P.al * ln( ( 1 - ul ) / ( 1 + ul ) ) / P.bl; // falta checar si log = ln
  us := us - P.u_0;
  If P.rot Then
  Begin
    xy.x := us;
    xy.y := vs;
  End
  Else
  Begin
    xy.x := vs * P.cosrot + us * P.sinrot;
    xy.y := us * P.cosrot - vs * P.sinrot;
  End;
  result := xy;
End;

Function omerc_inverse( const xy: TXY; P: TEzGeoConvert ): TLP; // ellipsoid & spheroid */
Var
  q, s, ul, us, vl, vs: double;
  lp: TLP;
Begin
  If P.rot Then
  Begin
    us := xy.x;
    vs := xy.y;
  End
  Else
  Begin
    vs := xy.x * P.cosrot - xy.y * P.sinrot;
    us := xy.y * P.cosrot + xy.x * P.sinrot;
  End;
  us := us + P.u_0;
  q := exp( -P.bl * vs / P.al );
  s := 0.5 * ( q - 1 / q );
  vl := sin( P.bl * us / P.al );
  ul := 2 * ( vl * P.cosgam + s * P.singam ) / ( q + 1 / q );
  If ( abs( abs( ul ) - 1 ) < EPS ) Then
  Begin
    lp.lam := 0;
    If ul < 0 Then
      lp.phi := -HALFPI
    Else
      lp.phi := HALFPI;
  End
  Else
  Begin
    lp.phi := P.el / sqrt( ( 1 + ul ) / ( 1 - ul ) );
    If P.ellips Then
    Begin
      lp.phi := pj_phi2( power( lp.phi, 1 / P.bl ), P.e );
      If ( lp.phi = HUGE_VAL ) Then
        exit; // I_ERROR;
    End
    Else
      lp.phi := HALFPI - 2 * abs( lp.phi );
    lp.lam := -aatan2( ( s * P.cosgam - vl * P.singam ), cos( P.bl * us / P.al ) ) / P.bl;
  End;
  result := lp;
End;

Procedure omerc( P: TEzGeoConvert; init: boolean );
Var
  con, com, cosph0, d, f, h, l, sinph0, pd, j: double;
  azi: boolean;
Begin
  P.rot := P.pj_param.asboolean( 'no_rot' ) = false;
  azi := P.pj_param.defined( 'alpha' );
  If azi Then
  Begin
    P.lamc := P.pj_param.asradians( 'lonc' );
    P.alpha := P.pj_param.asradians( 'alpha' );
    If ( ( abs( P.alpha ) <= TOL7 ) Or
      ( abs( abs( P.phi0 ) - HALFPI ) <= TOL7 ) Or
      ( abs( abs( P.alpha ) - HALFPI ) <= TOL7 ) ) Then
      Raise Exception.Create( SErr32 );
    //E_ERROR(-32);
  End
  Else
  Begin
    P.lam1 := P.pj_param.asradians( 'lon_1' );
    P.phi1 := P.pj_param.asradians( 'lat_1' );
    P.lam2 := P.pj_param.asradians( 'lon_2' );
    P.phi2 := P.pj_param.asradians( 'lat_2' );
    con := abs( P.phi1 );
    If ( ( abs( P.phi1 - P.phi2 ) <= TOL7 ) Or
      ( con <= TOL7 ) Or
      ( abs( con - HALFPI ) <= TOL7 ) Or
      ( abs( abs( P.phi0 ) - HALFPI ) <= TOL7 ) Or
      ( abs( abs( P.phi2 ) - HALFPI ) <= TOL7 ) ) Then
      Raise Exception.Create( SErr33 );
    //  E_ERROR(-33);
  End;
  P.ellips := P.es > 0;
  If P.ellips Then
    com := sqrt( P.one_es )
  Else
    com := 1;
  If abs( P.phi0 ) > EPS Then
  Begin
    sinph0 := sin( P.phi0 );
    cosph0 := cos( P.phi0 );
    If ( P.ellips ) Then
    Begin
      con := 1 - P.es * sinph0 * sinph0;
      P.bl := cosph0 * cosph0;
      P.bl := sqrt( 1 + P.es * P.bl * P.bl / P.one_es );
      P.al := P.bl * P.k0 * com / con;
      d := P.bl * com / ( cosph0 * sqrt( con ) );
    End
    Else
    Begin
      P.bl := 1;
      P.al := P.k0;
      d := 1 / cosph0;
    End;
    f := d * d - 1;
    If ( f <= 0 ) Then
      f := 0
    Else
    Begin
      f := sqrt( f );
      If ( P.phi0 < 0 ) Then
        f := -f;
    End;
    f := f + d;
    P.el := f;
    If ( P.ellips ) Then
      P.el := P.el * power( pj_tsfn( P.phi0, sinph0, P.e ), P.bl )
    Else
      P.el := P.el * TSFN0( P.phi0 );
  End
  Else
  Begin
    P.bl := 1 / com;
    P.al := P.k0;
    f := 1;
    d := 1;
    P.el := 1;
  End;
  If azi Then
  Begin
    P.Gamma := arcsin( sin( P.alpha ) / d );
    P.lam0 := P.lamc - arcsin( ( 0.5 * ( f - 1 / f ) ) * tan( P.Gamma ) ) / P.bl;
  End
  Else
  Begin
    If P.ellips Then
    Begin
      h := power( pj_tsfn( P.phi1, sin( P.phi1 ), P.e ), P.bl );
      l := power( pj_tsfn( P.phi2, sin( P.phi2 ), P.e ), P.bl );
    End
    Else
    Begin
      h := TSFN0( P.phi1 );
      l := TSFN0( P.phi2 );
    End;
    f := P.el / h;
    pd := ( l - h ) / ( l + h );
    j := P.el * P.el;
    j := ( j - l * h ) / ( j + l * h );
    con := P.lam1 - P.lam2;
    If ( con < -PI ) Then
      P.lam2 := P.lam2 - TWOPI
    Else If ( con > PI ) Then
      P.lam2 := P.lam2 + TWOPI;
    P.lam0 := adjlon( 0.5 * ( P.lam1 + P.lam2 ) - arctan(
      j * tan( 0.5 * P.bl * ( P.lam1 - P.lam2 ) ) / pd ) / P.bl );
    P.Gamma := arctan( 2 * sin( P.bl * adjlon( P.lam1 - P.lam0 ) ) /
      ( f - 1 / f ) );
    P.alpha := arcsin( d * sin( P.Gamma ) );
  End;
  P.singam := sin( P.Gamma );
  P.cosgam := cos( P.Gamma );
  If P.pj_param.asboolean( 'rot_conv' ) Then
    f := P.Gamma
  Else
    f := P.alpha;
  P.sinrot := sin( f );
  P.cosrot := cos( f );
  If P.pj_param.asboolean( 'no_uoff' ) Then
    P.u_0 := 0
  Else
    P.u_0 := abs( P.al * arctan( sqrt( d * d - 1 ) / P.cosrot ) / P.bl );
  If P.phi0 < 0 Then
    P.u_0 := -P.u_0;
  P.inv := @omerc_inverse;
  P.fwd := @omerc_forward;
End;

{ ***************** aea - Albers Equal Area *****************}

// determine latitude angle phi-1
Const
  TOL10 = 1.0E-10;

Function phi1_( Const qs, Te, Tone_es: double ): double;
Var
  i: integer;
  Phi, sinpi, cospi, con, com, dphi: double;
Begin
  Phi := arcsin( 0.5 * qs );
  If Te < EPSILON Then
  Begin
    result := Phi;
    exit;
  End;
  i := N_ITER;
  Repeat
    sinpi := sin( Phi );
    cospi := cos( Phi );
    con := Te * sinpi;
    com := 1.0 - con * con;
    dphi := 0.5 * com * com / cospi * ( qs / Tone_es -
      sinpi / com + 0.5 / Te * ln( ( 1.0 - con ) / ( 1.0 + con ) ) );
    Phi := Phi + dphi;
  Until ( abs( dphi ) <= TOL10 ) Or ( i <= 0 );
  If i > 0 Then
    result := Phi
  Else
    result := HUGE_VAL;
End;

Function e_forward_aea( lp: TLP; P: TEzGeoConvert ): TXY; (* ellipsoid & spheroid *)
Var
  tmp: double;
  xy: TXY;
Begin
  If P.ellips Then
    tmp := P.n * pj_qsfn( sin( lp.phi ), P.e, P.one_es )
  Else
    tmp := P.n2 * sin( lp.phi );
  P.rho := P.c - tmp;
  If P.rho < 0 Then
  Begin
    P.pj_errno := -20;
    result := xy;
    exit;
  End;
  P.rho := P.dd * sqrt( P.rho );
  lp.lam := lp.lam * P.n;
  xy.x := P.rho * sin( lp.lam );
  xy.y := P.rho0 - P.rho * cos( lp.lam );
  result := xy;
End;