ezprojections.pas

很管用的GIS控件

  t := abs( lp.phi ) - HALFPI;
  If ( t > EPS ) Or ( abs( lp.lam ) > 10 ) Then
  Begin
    xy.y := HUGE_VAL;
    xy.x := xy.y;
    pj_errno := -14;
  End
  Else
  Begin { proceed with projection }
    pj_errno := 0;
    If abs( t ) <= EPS Then
    Begin
      If lp.phi < 0 Then
        lp.phi := -HALFPI
      Else
        lp.phi := HALFPI;
    End
    Else If geoc Then
      lp.phi := arctan( rone_es * tan( lp.phi ) );
    lp.lam := lp.lam - lam0; { compute del lp.lam }
    If Not over Then
      lp.lam := adjlon( lp.lam ); { adjust del longitude }
    xy := TForward( fwd )( lp, self ); { project }
    If pj_errno <> 0 Then
    Begin
      xy.y := HUGE_VAL;
      xy.x := xy.y;
    End
    Else
    Begin
      { adjust for major axis and easting/northings }
      xy.x := fr_meter * ( a * xy.x + x0 );
      xy.y := fr_meter * ( a * xy.y + y0 );
    End;
  End;
  result := xy;
End;

// ------------------------------------------------------------------ //
{ inverse projection entry }

resourcestring
  serr0 = 'Inverse conversion not implemented';
  
Function TEzGeoConvert.pj_inv( Var xy: TXY ): TLP;
Var
  lp: TLP;
Begin
  { can't do as much preliminary checking as with forward }
  If inv = Nil Then
    Raise Exception.Create( serr0 );
  If ( xy.x = HUGE_VAL ) Or ( xy.y = HUGE_VAL ) Then
  Begin
    lp.phi := HUGE_VAL;
    lp.lam := lp.phi;
    pj_errno := -15;
    exit;
  End;
  pj_errno := 0;
  xy.x := ( xy.x * to_meter - x0 ) * ra; { de-scale and de-offset }
  xy.y := ( xy.y * to_meter - y0 ) * ra;
  lp := TInverse( inv )( xy, Self ); { inverse project }
  If pj_errno <> 0 Then
  Begin
    lp.phi := HUGE_VAL;
    lp.lam := lp.phi;
  End
  Else
  Begin
    lp.lam := lp.lam + lam0; { reduce from del lp.lam }
    If Not over Then
      lp.lam := adjlon( lp.lam ); { adjust longitude to CM }
    If geoc And ( abs( abs( lp.phi ) - HALFPI ) > EPS ) Then
      lp.phi := arctan( one_es * tan( lp.phi ) );
  End;
  result := lp;
End;

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

Procedure TEzGeoConvert.Geo_CoordSysInit( Params: TStringList );
Var
  name, s, temp: String;
  proj: TProj;
  pc: TEzProjectionCode;
  found: boolean;
  pl: TStringList;

  Procedure SetUnits;
  Var
    i: TEzCoordsUnits;
    //en: String;
  Begin
    { set units }
    s := '';
    name := pj_param.AsString( 'units' );
    If Length( name ) > 0 Then
    Begin
      units := name;
      found := false;
      For i := low( pj_units ) To high( pj_units ) Do
      Begin
        //en := GetEnumName( System.TypeInfo( TEzCoordsUnits ), Ord( i ) );
        If AnsiCompareText( pj_units[i].ID{copy( en, 3, length( en ) )}, Name ) = 0 Then 
        Begin
          found := true;
          break;
        End;
      End;
      If Not found Then
      Begin
        pj_errno := -7;
        exit;
      End;
      s := FloatToStr( pj_units[i].to_meter );
    End;
    temp := pj_param.AsString( 'to_meter' );
    If Length( temp ) > 0 Then
      s := temp;
    If Length( s ) > 0 Then
    Begin
      to_meter := StrToFloat( s );
      fr_meter := 1 / to_meter;
    End
    Else
    Begin
      fr_meter := 1;
      to_meter := fr_meter;
    End;
  End;

Begin
  { this signals as not having a projection }
  fwd := Nil;
  inv := Nil;
  spc := Nil;

  fParaList.Assign( Params );

  pl := fParaList;
  pj_errno := 0;

  If pl.count = 0 Then
  Begin
    pj_errno := -1;
    exit;
  End;

  { find projection selection }
  name := pj_param.AsString( 'proj' );
  If Length( name ) = 0 Then
  Begin
    pj_errno := -4;
    exit;
  End;
  found := false;
  pc := ProjCodeFromID( name, found );
  {for pc := low(pj_list) to high(pj_list) do
  begin
     if AnsiCompareText(pj_list[pc].id,name)=0 then
     begin
        found := true;
        break;
     end;
  end;}
  If Not found Then
  Begin
    pj_errno := -5;
    exit;
  End;
  proj := pj_list[pc].proj;
  { initialize projection structure }
  proj( self, true );
  { set ellipsoid/sphere parameters }
  If pj_ell_set( a, es ) <> 0 Then
  Begin
    exit;
  End;
  e := sqrt( es );
  ra := 1 / a;
  one_es := 1 - es;
  If one_es = 0 Then
  Begin
    pj_errno := -6;
    exit;
  End;
  rone_es := 1 / one_es;
  { set PIN.geoc coordinate system }
  self.geoc := ( es <> 0 ) And pj_param.AsBoolean( 'geoc' );
  { over-ranging flag }
  self.over := pj_param.AsBoolean( 'over' );
  { central meridian }
  self.lam0 := pj_param.AsRadians( 'lon_0' );
  { central latitude }
  self.phi0 := pj_param.AsRadians( 'lat_0' );
  { false easting and northing }
  self.x0 := pj_param.asfloat( 'x_0' );
  self.y0 := pj_param.asfloat( 'y_0' );
  { general scaling factor }
  If pj_param.Defined( 'k_0' ) Then
    self.k0 := pj_param.asfloat( 'k_0' )
  Else If pj_param.Defined( 'k' ) Then
    self.k0 := pj_param.asfloat( 'k' )
  Else
    self.k0 := 1;
  If self.k0 <= 0 Then
  Begin
    pj_errno := -31;
    exit;
  End;

  { set units of projection }
  SetUnits;
  If pj_errno <> 0 Then
  Begin
    exit;
  End;

  { projection specific initialization }
  proj( self, false );
  //if pj_errno <> 0 then
  //begin
     // you can raise error here and clean up some memory if needed (not now)
  //end;
End;

Procedure TEzGeoConvert.Geo_CoordSysToLatLong( Const x, y: Double; Var Long, Lat: Double );
Var
  lp: TLP;
  xy: TXY;
Begin
  { returns Long, Lat in degrees }
  xy.x := x;
  xy.y := y;
  lp := pj_inv( xy );
  long := RadToDeg( lp.lam );
  lat := RadToDeg( lp.phi );
End;

Procedure TEzGeoConvert.Geo_CoordSysFromLatLong( Const Long, Lat: Double; Var x, y: Double );
Var
  lp: TLP;
  xy: TXY;
Begin
  { receives long,lat in degrees }
  lp.lam := DegToRad( long );
  lp.phi := DegToRad( lat );
  xy := pj_fwd( lp );
  x := xy.x;
  y := xy.y;
End;

{calculate distance from (long1,lat1) to (long2,lat2)
 geodetic inverse problem is used}

Function TEzGeoConvert.Geo_Distance( Const Long1, Lat1, Long2, Lat2: Double ): Double;
Var
  faz, baz, b: double; //Forward and Backward azimuth
Begin
  b := a * sqrt( 1 - es );
  invert1( DegToRad( lat1 ),
    DegToRad( long1 ),
    DegToRad( lat2 ),
    DegToRad( long2 ),
    a,
    ( a - b ) / a,
    faz,
    baz,
    Result );
  Result := Result * fr_meter; // return in units configured, example feet
End;
// ----------------------------------------------------------------------

Initialization
  pj_list[aea].proj := @EzProjImpl.aea;
  //pj_list[aeqd     ].proj := nil;
  pj_list[airy].proj := @EzProjImpl.airy;
  //pj_list[aitoff   ].proj := nil;
  //pj_list[alsk     ].proj := nil;
  //pj_list[apian    ].proj := nil;
  //pj_list[august   ].proj := nil;
  //pj_list[bacon    ].proj := nil;
  //pj_list[bipc     ].proj := nil;
  //pj_list[boggs    ].proj := nil;
  pj_list[bonne].proj := @EzProjImpl.bonne;
  pj_list[cass].proj := @EzProjImpl.cass;
  pj_list[cc].proj := @EzProjImpl.cc;
  pj_list[cea].proj := @EzProjImpl.cea;
  //pj_list[chamb    ].proj := nil;
  //pj_list[collg    ].proj := nil;
  //pj_list[crast    ].proj := nil;
  //pj_list[denoy    ].proj := nil;
  //pj_list[eck1     ].proj := nil;
  //pj_list[eck2     ].proj := nil;
  //pj_list[eck3     ].proj := nil;
  pj_list[eck4].proj := @EzProjImpl.eck4;
  pj_list[eck5].proj := @EzProjImpl.eck5;
  pj_list[eck6].proj := @EzProjImpl.eck6;
  //pj_list[eqc      ].proj := nil;
  //pj_list[eqdc     ].proj := nil;
  //pj_list[euler    ].proj := nil;
  //pj_list[fahey    ].proj := nil;
  //pj_list[fouc     ].proj := nil;
  //pj_list[fouc_s   ].proj := @EzProjImpl.fouc_s;
  //pj_list[gall     ].proj := nil;
  //pj_list[gins8    ].proj := nil;
  pj_list[gn_sinu].proj := @EzProjImpl.gn_sinu;
  //pj_list[gnom     ].proj := nil;
  //pj_list[goode    ].proj := nil;
  //pj_list[gs48     ].proj := nil;
  //pj_list[gs50     ].proj := nil;
  //pj_list[hammer   ].proj := nil;
  //pj_list[hatano   ].proj := nil;
  pj_list[imw_p].proj := @EzProjImpl.imw_p;
  //pj_list[kav5     ].proj := nil;
  //pj_list[kav7     ].proj := nil;
  //pj_list[labrd    ].proj := nil;
  pj_list[laea].proj := @EzProjImpl.laea;
  //pj_list[lagrng   ].proj := nil;
  //pj_list[larr     ].proj := nil;
  //pj_list[lask     ].proj := nil;
  pj_list[lcc].proj := @EzProjImpl.lcc;
  pj_list[leac].proj := @EzProjImpl.leac;
  //pj_list[lee_os   ].proj := nil;
  //pj_list[loxim    ].proj := nil;
  //pj_list[lsat     ].proj := nil;
  //pj_list[mbt_s    ].proj := nil;
  //pj_list[mbt_fps  ].proj := nil;
  //pj_list[mbtfpp   ].proj := nil;
  //pj_list[mbtfpq   ].proj := nil;
  pj_list[mbtfps].proj := @EzProjImpl.mbtfps;
  pj_list[merc].proj := @EzProjImpl.merc;
  //pj_list[mil_os   ].proj := nil;
  pj_list[mill].proj := @EzProjImpl.mill;
  //pj_list[mpoly    ].proj := nil;
  pj_list[moll].proj := @EzProjImpl.moll;
  //pj_list[murd1    ].proj := nil;
  //pj_list[murd2    ].proj := nil;
  //pj_list[murd3    ].proj := nil;
  //pj_list[nell     ].proj := nil;
  //pj_list[nell_h   ].proj := nil;
  //pj_list[nicol    ].proj := nil;
  //pj_list[nsper    ].proj := nil;
  //pj_list[nzmg     ].proj := nil;
  //pj_list[ob_tran  ].proj := nil;
  //pj_list[ocea     ].proj := nil;
  //pj_list[oea      ].proj := nil;
  pj_list[omerc].proj := @EzProjImpl.omerc;
  //pj_list[ortel    ].proj := nil;
  pj_list[ortho].proj := @EzProjImpl.ortho;
  //pj_list[pconic   ].proj := nil;
  pj_list[poly].proj := @EzProjImpl.poly;
  //pj_list[putp1    ].proj := nil;
  //pj_list[putp2    ].proj := nil;
  //pj_list[putp3    ].proj := nil;
  //pj_list[putp3p   ].proj := nil;
  //pj_list[putp4p   ].proj := nil;
  //pj_list[putp5    ].proj := nil;
  //pj_list[putp5p   ].proj := nil;
  //pj_list[putp6    ].proj := nil;
  //pj_list[putp6p   ].proj := nil;
  //pj_list[qua_aut  ].proj := nil;
  //pj_list[robin    ].proj := nil;
  //pj_list[rpoly    ].proj := nil;
  pj_list[sinu].proj := @EzProjImpl.sinu;
  //pj_list[somerc   ].proj := nil;
  pj_list[stere].proj := @EzProjImpl.stere;
  //pj_list[tcc      ].proj := nil;
  pj_list[tcea].proj := @EzProjImpl.tcea;
  //pj_list[tissot   ].proj := nil;
  pj_list[tmerc].proj := @EzProjImpl.tmerc;
  pj_list[tpeqd].proj := @EzProjImpl.tpeqd;
  //pj_list[tpers    ].proj := nil;
  pj_list[ups].proj := @EzProjImpl.ups;
  //pj_list[urm5     ].proj := nil;
  //pj_list[urmfps   ].proj := nil;
  pj_list[utm].proj := @EzProjImpl.utm;
  //pj_list[vandg    ].proj := nil;
  //pj_list[vandg2   ].proj := nil;
  //pj_list[vandg3   ].proj := nil;
  pj_list[vandg4].proj := @EzProjImpl.vandg4;
  //pj_list[vitk1    ].proj := nil;
  //pj_list[wag1     ].proj := nil;
  //pj_list[wag2     ].proj := nil;
  //pj_list[wag3     ].proj := nil;
  pj_list[wag4].proj := @EzProjImpl.wag4;
  pj_list[wag5].proj := @EzProjImpl.wag5;
  //pj_list[wag6     ].proj := nil;
  pj_list[wag7].proj := @EzProjImpl.wag7;
  //pj_list[weren    ].proj := nil;
  //pj_list[wink1    ].proj := nil;
  //pj_list[wink2    ].proj := nil;
  //pj_list[wintri   ].proj := nil;

End.