PureBasic Forums - English

Code: Select all

type
  TGeoidGrid = Array[0..719, 0..1439] of SmallInt;

var
  GeoidData: TGeoidGrid;
  dH,H,H0: Double;

// Reading the WW15MGH.DAC file
Procedure LoadGeoidFile(const FileName: string);
var
  fs: TFileStream;
  i, j: Integer;
  tmp: SmallInt;
begin
  //create a stream To Read the file
  fs:= TFileStream.Create('WW15MGH.DAC', fmOpenRead Or fmShareDenyWrite);
  try
    //go through all the rows (latitude) And columns (longitude)
    For i := 0 To 719 do
      For j := 0 To 1439 do
      begin
        //Reading 2 bytes (SmallInt) from the stream
        fs.Read(tmp, SizeOf(SmallInt));
        //Convert byte order (Big Endian into Little Endian)
        tmp:=((tmp And $FF) shl 8) Or ((tmp And $FF00) shr 8);
        //save the value into an Array
        GeoidData[i,j]:=tmp;
      End;
  finally
    //Release of the stream
    fs.Free;
  End;
End;

Code: Select all

EnableExplicit

Define.i File, i, j, result = #False
Define *File, *Ptr.Word, *FileEnd
Define Filename$

Dim GeoidData.w(719, 1439)

;Filename$ = OpenFileRequester("Choose a file", "", "All|*.*", 0)
Filename$ = "c:\WW15MGH.DAC"
If Filename$
  File = ReadFile(#PB_Any, Filename$)
  If File
    *File = AllocateMemory(Lof(File), #PB_Memory_NoClear)
    If *File
      If ReadData(File, *File, MemorySize(*File)) = MemorySize(*File)
        result = #True
      EndIf
    EndIf
    CloseFile(File)
  EndIf
EndIf

If result = #True
  *Ptr = *File
  ;Temporary limited. The file size is 2 Mb.
  *FileEnd = *File + 1000 ;+ MemorySize(*File) - 1
  While *Ptr <= *FileEnd
    ;Convert byte order (Big Endian into Little Endian)
    GeoidData(i, j) = ((*Ptr\w & $FF) << 8) | ((*Ptr\w & $FF00) >> 8)
    *Ptr + 2  ; SizeOf(Word) = 2 bytes
    i+1 : j+1
  Wend
  FreeMemory(*File)
EndIf

Code: Select all

*FileEnd = *File + MemorySize(*File) - 1

Code: Select all

*Ptr.Unicode

Code: Select all

EnableExplicit


Define.i File, i, j, result, FileSize
Define *File, *Ptr.Word
Define Filename$

Dim GeoidData.w(719, 1439)

;Filename$ = OpenFileRequester("Choose a file", "", "All|*.*", 0)
Filename$ = "WW15MGH.DAC"
If Filename$
  
  *File = #Null
  
  FileSize = FileSize(Filename$)
  If FileSize >= 720 * 1440 * 2  ; 2.073.600
    File = ReadFile(#PB_Any, Filename$)
    If File
      *File = AllocateMemory(FileSize, #PB_Memory_NoClear)
      If *File
        If ReadData(File, *File, FileSize) = FileSize
          result = #True
        Else
          Debug "Was not able to read complete file"
        EndIf
      Else
        Debug "Was not able to allocate memory"
      EndIf
      CloseFile(File)
    Else
      Debug "Was not able to read file"
    EndIf
    
    If result = #True
      *Ptr = *File
      ;Temporary limited. The file size is 2 Mb.
      
      For i = 0 To 719
        For j = 0 To 1439
          ;Convert byte order (Big Endian into Little Endian)
          GeoidData(i, j) = ((*Ptr\w & $FF) << 8) | ((*Ptr\w & $FF00) >> 8)
          *Ptr + 2  ; SizeOf(Word) = 2 bytes
        Next j
      Next i
    EndIf
    
    If *File
      FreeMemory(*File)
    EndIf
    
  Else
    Debug "File too small (not valid): " + Str(FileSize)
  EndIf
  
EndIf

Code: Select all

GeoidData(i, j) = ((*Ptr\w) << 8) | ((*Ptr\w) >> 8)

; or use Macro

Macro BSwap16(_WordValue)
   ((_WordValue<<8 | _WordValue>>8) & $FFFF)
EndMacro

Code: Select all

Macro BSwap16(_WordValue)
   ((_WordValue<<8 | _WordValue>>8) & $FFFF)
EndMacro

a.w = $F1F2

b.u = ((a) << 8) | ((a) >> 8)
Debug Hex(b)  ; FFF1

b = BSwap16(a)
Debug Hex(b)  ; FFF1

b = ((a) << 8) | (((a) >> 8) & $FF)
Debug Hex(b)  ; F2F1



c.u = $F1F2

d.u = ((c) << 8) | ((c) >> 8)
Debug Hex(d)  ; F2F1

d = BSwap16(c)
Debug Hex(d)  ; F2F1

Code: Select all

  ; Swaps the 2 Bytes of a 16Bit value
  ; use it : result = BSwap16(x); x=[0..$FFFF])
  ; because of Arithmetic right Shift we need to Mask out irrelevant Bits
  Macro BSwap16(_WordValue)
    ((_WordValue & $FF00)<<8 | (_WordValue & $FF )>>8))
  EndMacro

Code: Select all

Procedure bswap16(value.u)
  CompilerIf #PB_Compiler_Backend = #PB_Backend_C
    !return __builtin_bswap16(v_value);
  CompilerElse
    !xor eax,eax
    !mov ax, word [p.v_value]
    !rol ax, 8
    ProcedureReturn
  CompilerEndIf
EndProcedure

Code: Select all

Procedure Show()
  Protected f.d, a.d, b.d, e.d

  f = 1.0 / ValD(GetGadgetText(#Edit8)) ;polar compression
  a = ValD(GetGadgetText(#Edit7)) ;Big axle shaft
  b = a * (1 - f) ;Small axle shaft
  e = 2 * f - f * f ;square of eccentricity

  ; loading a geoid file
  ; ...

  SetGadgetText(#Edit9,  StrD(b))
  SetGadgetText(#Edit10, StrD(e))

  CheckBox1Click()
EndProcedure


; Converting WGS-84 to ECEF coordinates

Procedure Button1Click()
  Protected a.d, e.d
  Protected B.d, L.d
  Protected x.d, y.d, z.d, N.d ;abscissa, ordinate, applicate (in meters)

  ;latitude and longitude (in degrees)
  B = ValD(GetGadgetText(#Edit2))
  L = ValD(GetGadgetText(#Edit1))

  ; height correction
  If GetGadgetState(#CheckBox1)
    H0 = ValD(GetGadgetText(#Edit3)) ;orthometric height
    H  = ValD(GetGadgetText(#Edit11)) ;geodetic height
  Else
    H  = ValD(GetGadgetText(#Edit3))
    H0 = ValD(GetGadgetText(#Edit11))
  EndIf

  ; degrees to radians
  B = Radian(B)
  L = Radian(L)

  e = ValD(GetGadgetText(#Edit10));square of eccentricity
  a = ValD(GetGadgetText(#Edit7)) ;Big axle shaft

  ; radius of curvature of the first vertical
  N = a / Sqr(1 - e * Pow(Sin(B), 2))

  ; ECEF coordinates
  x = (N + H) * Cos(B) * Cos(L)
  y = (N + H) * Cos(B) * Sin(L)
  z = (N * (1 - e) + H) * Sin(B)

  ; result
  SetGadgetText(#Edit4, StrD(x)) ;abscissa
  SetGadgetText(#Edit5, StrD(y)) ;ordinate
  SetGadgetText(#Edit6, StrD(z)) ;applicate
EndProcedure

Procedure Button2Click()
  Protected a.d, e.d
  Protected B.d, B0.d, L.d
  Protected x.d, y.d, z.d
  Protected N.d, p.d

  e = ValD(GetGadgetText(#Edit10))
  a = ValD(GetGadgetText(#Edit7))

  x = ValD(GetGadgetText(#Edit4))
  y = ValD(GetGadgetText(#Edit5))
  z = ValD(GetGadgetText(#Edit6))

  ; longitude
  L = ATan2(y, x)

  ; auxiliary parameter
  p = Sqr(x * x + y * y) ; p:=sqrt(Sqr(x)+Sqr(y))  in the original

  ; initial latitude value
  B = ATan2(z, p * (1 - e))

  ; iterations for calculating latitude
  Repeat
    B0 = B
    N  = a / Sqr(1 - e * Sin(B) * Sin(B))
    B  = ATan2(z + e * N * Sin(B), p)
  Until Abs(B0 - B) < ValD(GetGadgetText(#ComboBox1)) ; continue until we get the required accuracy

  ; geodetic height
  H = p * Cos(B) + z * Sin(B) - a * Sqr(1 - e * Sin(B) * Sin(B))

  ; radians to degrees
  L = Degree(L)
  B = Degree(B)

  ; geoid height
  dH = GetGeoidHeight(B, L)

  ; orthometric height
  H0 = H - dH

  ; output of coordinates
  SetGadgetText(#Edit2, StrD(B))
  SetGadgetText(#Edit1, StrD(L))

  ; height correction
  If GetGadgetState(#CheckBox1)
    SetGadgetText(#Edit3,  StrD(H0))
    SetGadgetText(#Edit11, StrD(H))
  Else
    SetGadgetText(#Edit3,  StrD(H))
    SetGadgetText(#Edit11, StrD(H0))
  EndIf

  SetGadgetText(#Edit12, StrD(dH))
EndProcedure