ezprojections.pas

很管用的GIS控件

  While i <= length( work ) Do
  Begin
    If work[i] In ['D', #39, #34] Then
    Begin
      Val( copy( work, last, ( i - last ) ), tv, Code );
      Case work[i] Of
        'D':
          Begin
            v := v + tv * 0.0174532925199433;
            nl := 1;
          End;
        #39: // '
          Begin
            v := v + tv * 0.0002908882086657216;
            nl := 2;
          End;
        #34: // "
          Begin
            v := v + tv * 0.0000048481368110953599;
            nl := 3;
          End;
      End;
      last := i + 1;
    End;
    inc( i );
  End;
  If last <= length( work ) Then
  Begin
    p := AnsiPos( work[length( work )], sym );
    If p > 0 Then
    Begin
      temp := copy( work, last, length( work ) - last );
      If p >= 4 Then
        sign := '-'
      Else
        sign := '+';
    End
    Else
    Begin
      temp := copy( work, last, length( work ) );
    End;
    If length( temp ) > 0 Then
    Begin
      val( temp, tv, code );
      Case nl Of
        0: v := v + tv * 0.0174532925199433;
        1: v := v + tv * 0.0002908882086657216;
        2: v := v + tv * 0.0000048481368110953599;
      End;
    End;
  End;
  If sign = '-' Then
    v := -v;
  result := v;

End;

{ SOLUTION OF THE GEODETIC INVERSE PROBLEM
  MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
  EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
  STANDPOINT / FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE

  A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
  F IS THE FLATTENING (NOT RECIPROCAL) OF THE REFERENCE ELLIPSOID
  LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST
  FORWARD AZIMUTHS AT BOTH POINTS RETURNED IN RADIANS FROM NORTH }

Procedure invert1( Const glat1, glon1, glat2, glon2, a, f: double;
  Var faz, baz, s: double );
Const
  EPS = 0.5E-13;
Var
  r, tu1, tu2, cu1, su1, cu2, x, sx, cx, sy, cy, y, sa, c2a, cz, e, c, d: double;
Begin
  r := 1.0 - f;
  tu1 := r * sin( glat1 ) / cos( glat1 );
  tu2 := r * sin( glat2 ) / cos( glat2 );
  cu1 := 1.0 / sqrt( tu1 * tu1 + 1.0 );
  su1 := cu1 * tu1;
  cu2 := 1.0 / sqrt( tu2 * tu2 + 1.0 );
  s := cu1 * cu2;
  baz := s * tu2;
  faz := baz * tu1;
  x := glon2 - glon1;
  Repeat
    sx := sin( x );
    cx := cos( x );
    tu1 := cu2 * sx;
    tu2 := baz - su1 * cu2 * cx;
    sy := sqrt( tu1 * tu1 + tu2 * tu2 );
    cy := s * cx + faz;
    y := Math.arctan2( sy, cy );
    sa := s * sx / sy;
    c2a := -sa * sa + 1.0;
    cz := faz + faz;
    If c2a > 0 Then
      cz := -cz / c2a + cy;
    e := cz * cz * 2.0 - 1.0;
    c := ( ( -3.0 * c2a + 4.0 ) * f + 4.0 ) * c2a * f / 16.0;
    d := x;
    x := ( ( e * cy * c + cz ) * sy * c + y ) * sa;
    x := ( 1.0 - c ) * x * f + glon2 - glon1;
  Until abs( d - x ) <= EPS;
  faz := Math.arctan2( tu1, tu2 );
  baz := Math.arctan2( cu1 * sx, baz * cx - su1 * cu2 ) + System.PI;
  x := sqrt( ( 1.0 / r / r - 1.0 ) * c2a + 1.0 ) + 1.0;
  x := ( x - 2.0 ) / x;
  c := 1.0 - x;
  c := ( x * x / 4.0 + 1.0 ) / c;
  d := ( 0.375 * x * x - 1.0 ) * x;
  x := e * cy;
  s := 1.0 - e - e;
  s := ( ( ( ( sy * sy * 4.0 - 3.0 ) * s * cz * d / 6.0 - x ) * d / 4.0 + cz ) * sy * d + y ) * c * a * r;
End;

Function ProjCodeFromID( Const ID: String; Var found: boolean ): TEzProjectionCode;
Var
  i: TEzProjectionCode;
Begin
  found := false;
  result := Low( TEzProjectionCode );
  For i := Low( pj_list ) To High( pj_list ) Do
    If AnsiCompareText( GetEnumName( System.TypeInfo( TEzProjectionCode ), Ord( i ) ), ID ) = 0 Then
    Begin
      result := i;
      found := true;
      exit;
    End;
End;

Function EllpsCodeFromID( Const ID: String ): TEzEllipsoidCode;
Var
  i: TEzEllipsoidCode;
  en: String;
Begin
  result := Low( TEzEllipsoidCode );
  For i := Low( TEzEllipsoidCode ) To High( TEzEllipsoidCode ) Do
  Begin
    en := GetEnumName( System.TypeInfo( TEzEllipsoidCode ), Ord( i ) );
    If AnsiCompareText( copy( en, 3, length( en ) ), ID ) = 0 Then //pj_ellps[i].ID
    Begin
      result := i;
      exit;
    End;
  End;
End;

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

Constructor TEzProjectParam.Create( GeoConvert: TEzGeoConvert );
Begin
  Inherited Create;
  fGeoConvert := GeoConvert;
End;

Function TEzProjectParam.Defined( Const opt: String ): boolean;
Begin
  result := Length( fGeoConvert.fParaList.Values[opt] ) > 0;
End;

Function TEzProjectParam.AsString( Const opt: String ): String;
Begin
  Result := AnsiLowerCase( fGeoConvert.fParaList.Values[opt] );
End;

Function TEzProjectParam.AsInteger( Const opt: String ): Integer;
Var
  Value: String;
Begin
  Value := AnsiLowerCase( fGeoConvert.fParaList.Values[opt] );
  If Length( Value ) > 0 Then
    result := StrToInt( Value )
  Else
    result := 0;
End;

Function TEzProjectParam.AsFloat( Const opt: String ): double;
Var
  Value: String;
  code: integer;
Begin
  Value := AnsiLowerCase( fGeoConvert.fParaList.Values[opt] );
  If Length( Value ) > 0 Then
    val( Value, result, code )
  Else
    result := 0;
End;

Function TEzProjectParam.AsBoolean( Const opt: String ): boolean;
Begin
  result := AnsiCompareText( fGeoConvert.fParaList.Values[opt], 't' ) = 0;
End;

Function TEzProjectParam.AsRadians( Const opt: String ): double;
Var
  Value: String;
Begin
  Value := AnsiLowerCase( fGeoConvert.fParaList.Values[opt] );
  If Length( Value ) > 0 Then
    result := TEzProjector.dmstor( Value )
  Else
    result := 0;
End;

//  -------------------------------------------------------------------- //
Const
  SIXTH = 0.1666666666666666667; // 1/6
  RA4 = 0.04722222222222222222; // 17/360
  RA6 = 0.02215608465608465608; // 67/3024
  RV4 = 0.06944444444444444444; // 5/72
  RV6 = 0.04243827160493827160; // 55/1296

Constructor TEzGeoConvert.Create;
Begin
  Inherited Create;
  fpj_param := TEzProjectParam.Create( Self );
  fParaList := TStringList.Create;
  { default parameters }
  fParaList.Add( 'units=m' );
  fParaList.Add( 'ellps=WGS84' );
  fParaList.Add( 'proj=utm' );
  fParaList.Add( 'zone=12' ); // our zone :-)
End;

Destructor TEzGeoConvert.Destroy;
Begin
  fpj_param.free;
  fParaList.free;
  Inherited Destroy;
End;

{ initialize geographic shape parameters }

Function TEzGeoConvert.pj_ell_set( Var a, es: double ): integer;
Label
  bomb;
Var
  i, last: Integer;
  b, e: double;
  name, s, en: String;
  found, defined: boolean;
  tmp, tmp2: double;
  pl: TStringList;
  c: TEzEllipsoidCode;
Begin
  Result := 0;
  pl := fParaList;

  last := pl.count;
  // check for varying forms of ellipsoid input
  a := 0;
  es := 0;
  // R takes precedence
  If pj_param.Defined( 'R' ) Then
    a := pj_param.AsFloat( 'R' )
  Else
  Begin { probable elliptical figure }
    name := pj_param.AsString( 'ellps' );
    { check if ellps present and temporarily append its values to pl }
    If Length( name ) > 0 Then
    Begin
      found := false;
      For c := low( TEzEllipsoidCode ) To high( TEzEllipsoidCode ) Do
      Begin
        en := GetEnumName( System.TypeInfo( TEzEllipsoidCode ), Ord( c ) );
        If AnsiCompareText( copy( en, 3, length( en ) ), Name ) = 0 Then //pj_ellps[i].ID
        Begin
          pl.add( pj_ellps[c].major );
          pl.add( pj_ellps[c].ell );
          found := true;
          break;
        End;
      End;
      If Not found Then
      Begin
        pj_errno := -9;
        exit;
      End;
    End;
    a := pj_param.AsFloat( 'a' );
    b := 0;
    If pj_param.Defined( 'es' ) Then { eccentricity squared }
      es := pj_param.AsFloat( 'es' )
    Else If pj_param.Defined( 'e' ) Then
    Begin { eccentricity }
      e := pj_param.AsFloat( 'e' );
      es := e * e;
    End
    Else If pj_param.Defined( 'rf' ) Then
    Begin { reciprocal flattening }
      es := pj_param.AsFloat( 'rf' );
      If es = 0 Then
      Begin
        pj_errno := -10;
        Goto bomb;
      End;
      es := 1 / es;
      es := es * ( 2 - es );
    End
    Else If pj_param.Defined( 'f' ) Then
    Begin { flattening }
      es := pj_param.AsFloat( 'f' );
      es := es * ( 2 - es );
    End
    Else If pj_param.Defined( 'b' ) Then
    Begin { minor axis }
      b := pj_param.AsFloat( 'b' );
      es := 1 - ( b * b ) / ( a * a );
    End; { else es = 0 and sphere of radius a }
    If b = 0 Then
      b := a * sqrt( 1 - es );
    { following options turn ellipsoid into equivalent sphere }
    If pj_param.AsBoolean( 'R_A' ) Then
    Begin { sphere -- area of ellipsoid }
      a := a * ( 1 - es * ( SIXTH + es * ( RA4 + es * RA6 ) ) );
      es := 0;
    End
    Else If pj_param.AsBoolean( 'R_V' ) Then
    Begin { sphere -- vol. of ellipsoid }
      a := a * ( 1 - es * ( SIXTH + es * ( RV4 + es * RV6 ) ) );
      es := 0.0;
    End
    Else If pj_param.AsBoolean( 'R_a' ) Then
    Begin { sphere -- arithmetic mean }
      a := 0.5 * ( a + b );
      es := 0;
    End
    Else If pj_param.AsBoolean( 'R_g' ) Then
    Begin { sphere -- geometric mean }
      a := sqrt( a * b );
      es := 0.0;
    End
    Else If pj_param.AsBoolean( 'R_h' ) Then
    Begin { sphere -- harmonic mean }
      a := 2 * a * b / ( a + b );
      es := 0;
    End
    Else
    Begin
      defined := pj_param.Defined( 'R_lat_a' );
      If defined Or { sphere -- arith. }
      pj_param.Defined( 'R_lat_g' ) Then
      Begin { or geom. mean at latitude }
        If defined Then
          s := 'R_lat_a'
        Else
          s := 'R_lat_g';
        tmp := Sin( pj_param.AsFloat( s ) );
        If abs( tmp ) > HALFPI Then
        Begin
          pj_errno := -11;
          Goto bomb;
        End;
        tmp := 1 - es * tmp * tmp;
        If defined Then
          tmp2 := 0.5 * ( 1 - es + tmp ) / ( tmp * sqrt( tmp ) )
        Else
          tmp2 := sqrt( 1 - es ) / tmp;
        a := a * tmp2;
        es := 0;
      End;
    End;
    bomb:
    For i := last To pl.count - 1 Do
      pl.delete( last );
    If pj_errno <> 0 Then
    Begin
      result := 1;
      exit;
    End;
  End;
  { some remaining checks }
  If es < 0 Then
  Begin
    pj_errno := -12;
    result := 1;
    exit;
  End;
  If a <= 0 Then
  Begin
    pj_errno := -13;
    result := 1;
    exit;
  End;
  result := 0;
End;

// ----------------------------------------------------------------------
Const
  EPS = 1.0E-12;

  { forward projection entry }

Function TEzGeoConvert.pj_fwd( Var lp: TLP ): TXY;
Var
  xy: TXY;
  t: double;
Begin
  { check for forward and latitude or longitude overange }