ezprojimpl.pas

很管用的GIS控件

End;

{ tpeqd - two point equidistant }

Function s_forward_tpeqd( const lp: TLP; P: TEzGeoConvert ): TXY; // sphere
Var
  t, z1, z2, dl1, dl2, sp, cp: double;
  xy: TXY;
Begin
  sp := sin( lp.phi );
  cp := cos( lp.phi );
  dl1 := lp.lam + P.dlam2;
  z1 := aacos( P.sp1 * sp + P.cp1 * cp * cos( dl1 ) );
  dl2 := lp.lam - P.dlam2;
  z2 := aacos( P.sp2 * sp + P.cp2 * cp * cos( dl2 ) );
  z1 := z1 * z1;
  z2 := z2 * z2;
  t := z1 - z2;
  xy.x := P.r2z0 * t;
  t := P.z02 - t;
  xy.y := P.r2z0 * asqrt( 4.0 * P.z02 * z2 - t * t );
  If ( P.ccs * sp - cp * ( P.cs * sin( dl1 ) - P.sc * sin( dl2 ) ) ) < 0.0 Then
    xy.y := -xy.y;
  result := xy;
End;

Function s_inverse_tpeqd( const xy: TXY; P: TEzGeoConvert ): TLP; // sphere
Var
  cz1, cz2, s, d, cp, sp: double;
  lp: TLP;
Begin
  cz1 := cos( hypot( xy.y, xy.x + P.hz0 ) );
  cz2 := cos( hypot( xy.y, xy.x - P.hz0 ) );
  s := cz1 + cz2;
  d := cz1 - cz2;
  lp.lam := -aatan2( d, ( s * P.thz0 ) );
  lp.phi := aacos( hypot( P.thz0 * s, d ) * P.rhshz0 );
  If xy.y < 0.0 Then
    lp.phi := -lp.phi;
  // lam--phi now in system relative to P1--P2 base equator
  sp := sin( lp.phi );
  cp := cos( lp.phi );
  lp.lam := lp.lam - P.lp;
  s := cos( lp.lam );
  lp.phi := aasin( P.sa * sp + P.ca * cp * s );
  lp.lam := arctan2( cp * sin( lp.lam ), P.sa * cp * s - P.ca * sp ) + P.lamc;
  result := lp;
End;

Procedure tpeqd( P: TEzGeoConvert; init: boolean );
Var
  lam_1, lam_2, phi_1, phi_2, A12, pp: double;
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;

  // get control point locations
  phi_1 := P.pj_param.asradians( 'lat_1' );
  lam_1 := P.pj_param.asradians( 'lon_1' );
  phi_2 := P.pj_param.asradians( 'lat_2' );
  lam_2 := P.pj_param.asradians( 'lon_2' );
  If ( phi_1 = phi_2 ) And ( lam_1 = lam_2 ) Then
    Raise exception.create( serr25 );
  P.lam0 := adjlon( 0.5 * ( lam_1 + lam_2 ) );
  P.dlam2 := adjlon( lam_2 - lam_1 );
  P.cp1 := cos( phi_1 );
  P.cp2 := cos( phi_2 );
  P.sp1 := sin( phi_1 );
  P.sp2 := sin( phi_2 );
  P.cs := P.cp1 * P.sp2;
  P.sc := P.sp1 * P.cp2;
  P.ccs := P.cp1 * P.cp2 * sin( P.dlam2 );
  P.z02 := aacos( P.sp1 * P.sp2 + P.cp1 * P.cp2 * cos( P.dlam2 ) );
  P.hz0 := 0.5 * P.z02;
  A12 := arctan2( P.cp2 * sin( P.dlam2 ), P.cp1 * P.sp2 - P.sp1 * P.cp2 * cos( P.dlam2 ) );
  pp := aasin( P.cp1 * sin( A12 ) );
  P.ca := cos( pp );
  P.sa := sin( pp );
  P.lp := adjlon( arctan2( P.cp1 * cos( A12 ), P.sp1 ) - P.hz0 );
  P.dlam2 := p.dlam2 * 0.5;
  P.lamc := HALFPI - arctan2( sin( A12 ) * P.sp1, cos( A12 ) ) - P.dlam2;
  P.thz0 := tan( P.hz0 );
  P.rhshz0 := 0.5 / sin( P.hz0 );
  P.r2z0 := 0.5 / P.z02;
  P.z02 := P.z02 * P.z02;
  P.inv := @s_inverse_tpeqd;
  P.fwd := @s_forward_tpeqd;
  P.es := 0.0;
End;

{ tcea - Transverse Cylindrical Equal Area}

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

Function s_inverse_tcea( xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  t: double;
  lp: TLP;
Begin
  xy.y := xy.y * P.rk0 + P.phi0;
  xy.x := xy.x * P.k0;
  t := sqrt( 1.0 - xy.x * xy.x );
  lp.phi := arcsin( t * sin( xy.y ) );
  lp.lam := arctan2( xy.x, t * cos( xy.y ) );
  result := lp;
End;

Procedure tcea( P: TEzGeoConvert; init: boolean );
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  P.rk0 := 1 / P.k0;
  P.inv := @s_inverse_tcea;
  P.fwd := @s_forward_tcea;
  P.es := 0.0;
End;

{ vandg4 - van der Grinten IV}
Const

  TWORPI = 0.63661977236758134308;

Function s_forward_vandg4( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Var
  xy: TXY;
  x1, t, bt, ct, ft, bt2, ct2, dt, dt2: double;
Begin
  If abs( lp.phi ) < TOL10 Then
  Begin
    xy.x := lp.lam;
    xy.y := 0.0;
  End
  Else If ( abs( lp.lam ) < TOL10 ) Or ( abs( abs( lp.phi ) - HALFPI ) < TOL10 ) Then
  Begin
    xy.x := 0.0;
    xy.y := lp.phi;
  End
  Else
  Begin
    bt := abs( TWORPI * lp.phi );
    bt2 := bt * bt;
    ct := 0.5 * ( bt * ( 8.0 - bt * ( 2.0 + bt2 ) ) - 5.0 ) / ( bt2 * ( bt - 1.0 ) );
    ct2 := ct * ct;
    dt := TWORPI * lp.lam;
    dt := dt + 1. / dt;
    dt := sqrt( dt * dt - 4.0 );
    If ( ( abs( lp.lam ) - HALFPI ) < 0.0 ) Then
      dt := -dt;
    dt2 := dt * dt;
    x1 := bt + ct;
    x1 := x1 * x1;
    t := bt + 3 * ct;
    ft := x1 * ( bt2 + ct2 * dt2 - 1.0 ) + ( 1.0 - bt2 ) *
      ( bt2 * ( t * t + 4.0 * ct2 ) + ct2 * ( 12.0 * bt * ct + 4.0 * ct2 ) );
    x1 := ( dt * ( x1 + ct2 - 1.0 ) + 2.0 * sqrt( ft ) ) / ( 4.0 * x1 + dt2 );
    xy.x := HALFPI * x1;
    xy.y := HALFPI * sqrt( 1.0 + dt * abs( x1 ) - x1 * x1 );
    If ( lp.lam < 0.0 ) Then
      xy.x := -xy.x;
    If ( lp.phi < 0.0 ) Then
      xy.y := -xy.y;
  End;
  result := xy;
End;

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

{ laea - lambert azimuthal equal area }
Var
  sinph0: double;
  cosph0: double;

Function e_forward_laea( const lp: TLP; P: TEzGeoConvert ): TXY; // ellipsoid
Label
  eqcon;
Var
  xy: TXY;
  coslam, sinlam, sinphi, q, sinb, cosb, b: double;
Begin
  sinb := 0.0;
  cosb := 0.0;
  b := 0.0;
  coslam := cos( lp.lam );
  sinlam := sin( lp.lam );
  sinphi := sin( lp.phi );
  q := pj_qsfn( sinphi, P.e, P.one_es );
  If ( P.mode = OBLIQ ) Or ( P.mode = EQUIT ) Then
  Begin
    sinb := q / P.qp;
    cosb := sqrt( 1.0 - sinb * sinb );
  End;
  Case P.mode Of
    OBLIQ:
      b := 1.0 + P.sinb1 * sinb + P.cosb1 * cosb * coslam;
    EQUIT:
      b := 1.0 + cosb * coslam;
    N_POLE:
      Begin
        b := HALFPI + lp.phi;
        q := P.qp - q;
      End;
    S_POLE:
      Begin
        b := lp.phi - HALFPI;
        q := P.qp + q;
      End;
  End;
  If abs( b ) < EPS10 Then
    Raise exception.create( 'error' );
  Case P.mode Of
    OBLIQ:
      Begin
        b := sqrt( 2.0 / b );
        xy.y := P.ymf * b * ( P.cosb1 * sinb - P.sinb1 * cosb * coslam );
        Goto eqcon;
      End;
    EQUIT:
      Begin
        b := sqrt( 2.0 / ( 1.0 + cosb * coslam ) );
        xy.y := b * sinb * P.ymf;
        eqcon:
        xy.x := P.xmf * b * cosb * sinlam;
      End;
    N_POLE, S_POLE:
      If q >= 0.0 Then
      Begin
        b := sqrt( q );
        xy.x := b * sinlam;
        If P.mode = S_POLE Then
          xy.y := coslam * b
        Else
          xy.y := coslam * ( -b );
      End
      Else
      Begin
        xy.y := 0.0;
        xy.x := 0.0;
      End;
  End;
  result := xy;
End;

Function s_forward_laea( const lp: TLP; P: TEzGeoConvert ): TXY; // spheroid
Label
  oblcon;
Var
  xy: TXY;
  coslam, cosphi, sinphi: double;
Begin
  sinphi := sin( lp.phi );
  cosphi := cos( lp.phi );
  coslam := cos( lp.lam );
  Case P.mode Of
    EQUIT:
      Begin
        xy.y := 1.0 + cosphi * coslam;
        Goto oblcon;
      End;
    OBLIQ:
      Begin
        xy.y := 1.0 + sinph0 * sinphi + cosph0 * cosphi * coslam;
        oblcon:
        If xy.y <= EPS10 Then
          Raise exception.create( 'error' );
        xy.y := sqrt( 2.0 / xy.y );
        xy.x := xy.y * cosphi * sin( lp.lam );
        If P.mode = EQUIT Then
          xy.y := xy.y * sinphi
        Else
          xy.y := xy.y * ( cosph0 * sinphi - sinph0 * cosphi * coslam );
      End;
    N_POLE, S_POLE:
      Begin
        If P.Mode = N_POLE Then
          coslam := -coslam;
        If abs( lp.phi + P.phi0 ) < EPS10 Then
          Raise exception.create( 'error' );
        xy.y := FORTPI - lp.phi * 0.5;
        If P.mode = S_POLE Then
          xy.y := 2.0 * cos( xy.y )
        Else
          xy.y := 2.0 * sin( xy.y );
        xy.x := xy.y * sin( lp.lam );
        xy.y := xy.y * coslam;
      End;
  End;
  result := xy;
End;

Function e_inverse_laea( xy: TXY; P: TEzGeoConvert ): TLP; // ellipsoid
Var
  lp: TLP;
  cCe, sCe, q, rho, ab: double;
Begin
  ab := 0.0;
  Case P.mode Of
    EQUIT,
      OBLIQ:
      Begin
        xy.x := xy.x / P.dd;
        xy.y := xy.y * P.dd;
        rho := hypot( xy.x, xy.y );
        If rho < EPS10 Then
        Begin
          lp.lam := 0.0;
          lp.phi := P.phi0;
          result := lp;
          exit;
        End;
        sCe := 2.0 * arcsin( 0.5 * rho / P.rq );
        cCe := cos( sCe );
        sCe := sin( sCe );
        xy.x := xy.x * sCe;
        If P.mode = OBLIQ Then
        Begin
          ab := cCe * P.sinb1 + xy.y * sCe * P.cosb1 / rho;
          //q := P.qp * ab;
          xy.y := rho * P.cosb1 * cCe - xy.y * P.sinb1 * sCe;
        End
        Else
        Begin
          ab := xy.y * sCe / rho;
          //q := P.qp * ab;
          xy.y := rho * cCe;
        End;
      End;
    N_POLE, S_POLE:
      Begin
        If P.Mode = N_POLE Then
          xy.y := -xy.y;
        q := xy.x * xy.x + xy.y * xy.y;
        If q = 0 Then
        Begin
          lp.lam := 0.0;
          lp.phi := P.phi0;
          result := lp;
        End;
        (*
        q := P.qp - q;
        *)
        ab := 1.0 - q / P.qp;
        If P.mode = S_POLE Then
          ab := -ab;
      End;
  End;
  lp.lam := arctan2( xy.x, xy.y );
  lp.phi := pj_authlat( arcsin( ab ), P );
  result := lp;
End;

Function s_inverse_laea( xy: TXY; P: TEzGeoConvert ): TLP; // spheroid
Var
  cosz, rh, sinz: double;
  lp: TLP;
Begin
  cosz := 0.0;
  sinz := 0.0;
  rh := hypot( xy.x, xy.y );
  lp.phi := rh * 0.5;
  If lp.phi > 1.0 Then
    Raise exception.create( 'error' );
  lp.phi := 2.0 * arcsin( lp.phi );
  If ( P.mode = OBLIQ ) Or ( P.mode = EQUIT ) Then
  Begin
    sinz := sin( lp.phi );
    cosz := cos( lp.phi );
  End;
  Case P.mode Of
    EQUIT:
      Begin
        If abs( rh ) <= EPS10 Then
          lp.phi := 0
        Else
          lp.phi := arcsin( xy.y * sinz / rh );
        xy.x := xy.x * sinz;
        xy.y := cosz * rh;
      End;
    OBLIQ:
      Begin
        If abs( rh ) <= EPS10 Then
          lp.phi := P.phi0
        Else
          lp.phi := arcsin( cosz * sinph0 + xy.y * sinz * cosph0 / rh );
        xy.x := xy.x * ( sinz * cosph0 );
        xy.y := ( cosz - sin( lp.phi ) * sinph0 ) * rh;
      End;
    N_POLE:
      Begin
        xy.y := -xy.y;
        lp.phi := HALFPI - lp.phi;
      End;
    S_POLE:
      lp.phi := lp.phi - HALFPI;
  End;
  If ( xy.y = 0.0 ) And ( ( P.mode = EQUIT ) Or ( P.mode = OBLIQ ) ) Then
    lp.lam := 0
  Else
    lp.lam := arctan2( xy.x, xy.y );
  result := lp;
End;

Procedure laea( P: TEzGeoConvert; init: boolean );
Var
  t: double;
  sinphi: double;
Begin
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  sinph0 := P.sinb1;
  cosph0 := P.cosb1;
  If init Then
  Begin
    P.fwd := Nil;
    P.inv := Nil;
    exit;
  End;
  t := abs( P.phi0 );
  If abs( t - HALFPI ) < EPS10 Then
  Begin
    If P.phi0 < 0.0 Then
      P.mode := S_POLE
    Else
      P.mode := N_POLE;
  End
  Else If abs( t ) < EPS10 Then
    P.mode := EQUIT
  Else
    P.mode := OBLIQ;
  If P.es <> 0 Then
  Begin
    P.e := sqrt( P.es );
    P.qp := pj_qsfn( 1.0, P.e, P.one_es );
    P.mmf := 0.5 / ( 1.0 - P.es );
    pj_authset( P );
    Case P.mode Of
      N_POLE,
        S_POLE:
        P.dd := 1.0;
      EQUIT:
        Begin
          P.rq := sqrt( 0.5 * P.qp );
          P.dd := 1.0 / P.rq;
          P.xmf := 1.0;
          P.ymf := 0.5 * P.qp;
        End;
      OBLIQ:
        Begin
          P.rq := sqrt( 0.5 * P.qp );
          sinphi := sin( P.phi0 );
          P.sinb1 := pj_qsfn( sinphi, P.e, P.one_es ) / P.qp;
          P.cosb1 := sqrt( 1.0 - P.sinb1 * P.sinb1 );
          P.dd := cos( P.phi0 ) / ( sqrt( 1.0 - P.es * sinphi * sinphi ) * P.rq * P.cosb1 );
          P.xmf := P.rq;
          P.ymf := P.xmf / P.dd;
          P.xmf := P.dd;
        End;
    End;
    P.inv := @e_inverse_laea;
    P.fwd := @e_forward_laea;
  End
  Else
  Begin
    If P.mode = OBLIQ Then
    Begin
      sinph0 := sin( P.phi0 );
      cosph0 := cos( P.phi0 );
    End;
    P.inv := @s_inverse_laea;
    P.fwd := @s_forward_laea;
  End;
End;

{ stere, ups - Universal Polar Stereographic }
Const
  TOL8 = 1.E-8;
  NITER8 = 8;
  S_POLEa = 0;
  N_POLEa = 1;
  OBLIQa = 2;
  EQUITa = 3;
  CONV = 1E-10;

Function ssfn_( Const phit: double; Var sinphi: double; Const eccen: double ): double;
Begin
  sinphi := sinphi * eccen;
  result := ( tan( 0.5 * ( HALFPI + phit ) ) * power( ( 1.0 - sinphi ) / ( 1.0 + sinphi ), 0.5 * eccen ) );
End;

Function e_forward_stere( lp: TLP; P: TEzGeoConvert ): TXY; // ellipsoide
Label
  xmul;
Var
  xy: TXY;
  coslam, sinlam, sinX, cosX, X, A, sinphi: double;
Begin
  sinX := 0.0;
  cosX := 0.0;
  coslam := cos( lp.lam );
  sinlam := sin( lp.lam );
  sinphi := sin( lp.phi );
  If ( P.mode = OBLIQa ) Or ( P.mode = EQUITa ) Then
  Begin
    X := 2 * arctan( ssfn_