VERSION 5.00
Begin {C62A69F0-16DC-11CE-9E98-00AA00574A4F} frmDatumTran 
   Caption         =   "Datum Transformation Parameters"
   ClientHeight    =   5520
   ClientLeft      =   45
   ClientTop       =   330
   ClientWidth     =   4050
   OleObjectBlob   =   "frmDatumTran.frx":0000
   StartUpPosition =   1  'CenterOwner
End
Attribute VB_Name = "frmDatumTran"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Option Explicit
Private pMxDoc As IMxDocument
Private pMaps As IMaps
Private pMap As IMap
Private m_pEnumLayers As IEnumLayer
Private m_pFLayer As IFeatureLayer

Dim pFLayer1 As IFeatureLayer
Dim pFLayer2 As IFeatureLayer
Dim pSpRef_In As ISpatialReference
Dim pSpRFc As SpatialReferenceEnvironment
Dim pGC_In As IGeographicCoordinateSystem
Dim pGC_Out As IGeographicCoordinateSystem
Dim pSpRef_Out As ISpatialReference
Dim pTrans As ITransformation
Dim pDisplayTrans As IDisplayTransformation
Dim pUserTrans As IMolodenskyTransformation
Dim pGTOperationSet As IGeoTransformationOperationSet
Dim pCompGTSet As ICompositeGeoTransformation
Dim pMapGeoTrans As IMapGeographicTransformations
Dim pSpatialReferenceOld As ISpatialReference
Dim pFeatureClass As IFeatureClass
Dim pGeoDataset As IGeoDataset
Dim dA As Double
Dim dF As Double
Dim dE2 As Double 'spheroid params
Dim DX As Double
Dim DY As Double
Dim DZ As Double ' molodensky shifts

Private Sub cboMaps_Change()
    cmdFrameInfo.Enabled = cboMaps.Text <> ""
    lboMapLayers.Clear
    
    Dim i As Integer
    For i = 0 To pMaps.Count - 1
        Set pMap = pMaps.Item(i)
        If pMap.Name = cboMaps.Text Then
            Exit For
        End If
    Next i
    Set m_pEnumLayers = pMap.Layers
    Set m_pFLayer = m_pEnumLayers.Next
    Do Until m_pFLayer Is Nothing
        lboMapLayers.AddItem m_pFLayer.Name
        Set m_pFLayer = m_pEnumLayers.Next
    Loop
End Sub

Private Sub cmdClose_Click()
    Unload Me
End Sub

Private Sub CmdFindMol_Click()
Dim iFactCode As Long
  iFactCode = pSpatialReferenceOld.FactoryCode
  Set pGC_In = pSpRFc.CreateGeographicCoordinateSystem(iFactCode)
   Set pSpRef_In = pGC_In
   Set pGC_Out = pSpRFc.CreateGeographicCoordinateSystem(esriSRGeoCS_WGS1984)
  Set pSpRef_Out = pGC_Out
  Call FindDatumMolodParms
End Sub

Private Sub cmdFrameInfo_Click()
    MsgBox "Information for data frame " & pMap.Name & vbCrLf & _
         "*    The spatial reference is " & pMap.SpatialReference.Name & vbCrLf & _
         "*    It has " & pMap.LayerCount & " layers" & vbCrLf & _
         "*    It has " & pMap.SelectionCount & " features selected", vbInformation
End Sub

Private Sub cmdLayerInfo_Click()
    Set pFeatureClass = m_pFLayer.FeatureClass
    Dim pGeometry As IGeometry
    'Set pGeometry = m_pFLayer.ShapeType
    ' QI for the Geodataset from the FeatureClass
    Set pGeoDataset = pFeatureClass
    Set pSpatialReferenceOld = pGeoDataset.SpatialReference
    'MsgBox "Geometry type of sel Layer " & pGeometry.Name & vbCrLf & _

      MsgBox "Selected layer's coordsys is " & pSpatialReferenceOld.Name, vbInformation
End Sub

Private Sub cmdWriteToGeogTranFile_Click()
    'If DX * DY * DZ <> 0 Then
        Open "D:\Holland\geogtran" For Output As #1
        Print #1, "GEOGTRAN,208209,";
        Print #1, Chr(34); "DutchRD_To_WGS_FromCtrlPoints2"; Chr(34); ",";
        Print #1, "PE_GCS_AMERSFOORT,"; "PE_GCS_WGS_1984,";
        Print #1, "PE_MTH_MOLODENSKY,100040,";
        Print #1, FormatNumber(DX, 3);
        Print #1, ",100041,";
        Print #1, FormatNumber(DY, 3);
        Print #1, ",100042,";
        Print #1, FormatNumber(DZ, 3)
        Close #1
    'End If
End Sub

Private Sub lboMapLayers_Click()
    'lboFields.Clear
    cmdLayerInfo.Enabled = lboMapLayers.Text <> ""
    m_pEnumLayers.Reset
    Set m_pFLayer = m_pEnumLayers.Next
    Do Until m_pFLayer Is Nothing
        If m_pFLayer.Name = lboMapLayers.Text Then
            Exit Do
        End If
        Set m_pFLayer = m_pEnumLayers.Next
    Loop
    'Set m_pFLayer = m_pFLayer
    Set pFeatureClass = m_pFLayer.FeatureClass
    ' QI for the Geodataset from the FeatureClass
    Set pGeoDataset = pFeatureClass
    Set pSpatialReferenceOld = pGeoDataset.SpatialReference
      MsgBox m_pFLayer.Name + " SpatialRef of selected Layer is " + pSpatialReferenceOld.Name
End Sub

Private Sub obSpatRefChoice1_Click()
  Dim iFactCode As Long
  iFactCode = pSpatialReferenceOld.FactoryCode
  Set pGC_In = pSpRFc.CreateGeographicCoordinateSystem(iFactCode)
  Call TransformCoordsToWGS
  MsgBox "Points will be transformed to WGS "
End Sub

Private Sub obSpatRefChoice2_Click()
  Dim iFactCode As Long
  iFactCode = pSpatialReferenceOld.FactoryCode
  Set pGC_In = pSpRFc.CreateGeographicCoordinateSystem(iFactCode)
  Call UndoTransformCoordsToWGS
  MsgBox " Points will be set back to initial positions"
End Sub

Public Sub UserForm_Initialize()
    Set pMxDoc = ThisDocument
    Set pMaps = pMxDoc.Maps
    Dim i As Integer
    For i = 0 To pMaps.Count - 1
        cboMaps.AddItem pMaps.Item(i).Name
    Next i
    lblFocusMap.Caption = "Active Data Frame: " & pMxDoc.FocusMap.Name
    Set pMxDoc = Application.Document
    ' Start by getting a handle on the current FocusMap
    Set pMap = pMxDoc.FocusMap
    Set pSpRFc = New SpatialReferenceEnvironment
    Dim iLayerCount As Integer
    iLayerCount = pMap.LayerCount ''''moet nog verbeterd: alleen pointlayers tellen
    If iLayerCount < 2 Then
        MsgBox "load 2 layers"
        Exit Sub
    Else
        Set pFLayer1 = pMap.Layer(0)
        Set pFLayer2 = pMap.Layer(1)
  End If
  DX = 0
  DY = 0
  DZ = 0
End Sub

Private Sub TransformCoordsToWGS()
  '  Set the input spatial reference to CS_Choice, the cordsys of the selected layer
  'Set pGC_In = pSpRFc.CreateGeographicCoordinateSystem(CS_choice)
  'Set pGC_In = New GeographicCoordinateSystem
  Set pSpRef_In = pGC_In
  pSpRef_In.SetFalseOriginAndUnits -180, -90, 1000000
  '  Set the output spatial reference to WGS84
  Set pGC_Out = pSpRFc.CreateGeographicCoordinateSystem(esriSRGeoCS_WGS1984)
  Set pSpRef_Out = pGC_Out
  pSpRef_Out.SetFalseOriginAndUnits -180, -90, 1000000
  '  Define the usermade geographic transformation
  Set pMapGeoTrans = pMap
  Set pGTOperationSet = pMapGeoTrans.GeographicTransformations
  pGTOperationSet.RemoveAll
  Set pUserTrans = New MolodenskyTransformation
  pUserTrans.Name = "DutchRD_To_WGS_FromCtrlPoints"
  Call FindDatumMolodParms
  pUserTrans.PutParameters DX, DY, DZ
  'pUserTrans.PutParameters 565.036, 49.914, 465.839
  pUserTrans.PutSpatialReferences pGC_In, pGC_Out
  'Set pUserTrans = pSpRFc.CreateGeoTransformation(esriGeoT
  
  '  Add the transformation to the operation set
  Set pGTOperationSet = New GeoTransformationOperationSet
  '  The Map uses IMapGeographicTransformations to access the geogtrans
  Set pMapGeoTrans = pMap
  Set pGTOperationSet = pMapGeoTrans.GeographicTransformations
  
  pGTOperationSet.Set esriTransformForward, pUserTrans
  pGTOperationSet.Set esriTransformReverse, pUserTrans
  Dim DXtest, DYtest, DZtest
  'pUserTrans.GetParameters DXtest, DYtest, DZtest
  'Now refresh the map
  DoEvents
  pMxDoc.ActiveView.Refresh
End Sub
Private Sub UndoTransformCoordsToWGS()
  Dim pSpatRefFact2 As ISpatialReferenceFactory2
  Set pSpatRefFact2 = New SpatialReferenceEnvironment
  Set pMapGeoTrans = pMap
  Set pGTOperationSet = pMapGeoTrans.GeographicTransformations
  pGTOperationSet.RemoveAll
  Set pTrans = pSpatRefFact2.CreateGeoTransformation(esriSRGeoTransformation_NAD1927_To_WGS1984_1)
  pGTOperationSet.Set esriTransformForward, pTrans
  pGTOperationSet.Set esriTransformReverse, pTrans
  Dim pActiveview As IActiveView
  Set pActiveview = pMap
  pActiveview.Refresh
End Sub
Public Sub FindDatumMolodParms()
 'Open the 1st shape file
    Dim pFClass1 As IFeatureClass
    Set pFClass1 = pFLayer1.FeatureClass
    Dim pSpheroid As ISpheroid
    Set pSpheroid = pGC_In.Datum.Spheroid
    dA = pSpheroid.SemiMajorAxis
    dF = pSpheroid.Flattening
    dE2 = 2 * dF - dF * dF
    'Open the 2nd shape file
    Dim pFClass2 As IFeatureClass
    Set pFClass2 = pFLayer2.FeatureClass
    Dim pFeature1 As IFeature  ' Deze variabele is nodig om straks alle feature objecten uit de Search op te vangen en te behandelen
    'Set pFeature1 = New Feature   ' hier creeer je een nieuw punt. Dat is niet wat je wilt
    Dim pFeature2 As IFeature
    'Set pFeature2 = New Feature
    Dim pPoint1 As IPoint
    Dim pPoint2 As IPoint
    'Set pPoint2 = New Point
    Dim Lon1 As Double
    Dim Lat1 As Double
    Dim H1 As Double
    Dim Lon2 As Double
    Dim Lat2 As Double
    Dim H2 As Double
    Dim pField As IField
    Dim i, j As Long
    i = 0
    Dim LaLoHe1 As LatLonHei
    Dim LaLoHe2 As LatLonHei
    Set LaLoHe1 = New LatLonHei
    Set LaLoHe2 = New LatLonHei
    'LaLoHe1.Class_Initialize
    Dim pPointCursor1 As IFeatureCursor
    Dim pPointCursor2 As IFeatureCursor


    Dim numfeatures As Long
    numfeatures = pFClass1.FeatureCount(Nothing)
    
    Dim pFields As IFields
    Set pFields = pFClass1.Fields
    Dim pField2 As IField
    Dim numfields As Long
    numfields = pFields.FieldCount
    Dim i2 As Long
    Dim pGeometrydef As IGeometryDef
    Dim hasz As Boolean
    For i2 = 0 To numfields - 1
      Set pField2 = pFields.Field(i2)
      If pField2.Type = esriFieldTypeGeometry Then
        Set pGeometrydef = pField2.GeometryDef
        hasz = pGeometrydef.hasz
        Exit For
      End If
    Next i2
    Set pPointCursor1 = pFClass1.Search(Nothing, False)
    Set pFeature1 = pPointCursor1.NextFeature
    Set pPointCursor2 = pFClass2.Search(Nothing, False)
    Set pFeature2 = pPointCursor2.NextFeature
    Dim XYZ1(100, 2) As Double
    Dim XYZ2(100, 2) As Double
    Dim XYZShifts(100, 2) As Double
    Do Until pFeature1 Is Nothing
        If i > 18 Then Exit Do
        Set pPoint1 = pFeature1.Shape
        Set pPoint2 = pFeature2.Shape
        If pPoint1 Is Nothing Then End
        If pPoint2 Is Nothing Then End
        Lon1 = pPoint1.X
        Lat1 = pPoint1.Y
        H1 = pFeature1.Value(pPointCursor1.FindField("VALUE"))
        Lon2 = pPoint2.X
        Lat2 = pPoint2.Y
        H2 = pFeature2.Value(pPointCursor2.FindField("VALUE"))
        LaLoHe1.SetLat (Lat1)
         LaLoHe1.SetLon (Lon1)
         LaLoHe1.SetHei (H1)
         LaLoHe1.LLH_to_XYZ dA, dE2
         XYZ1(i, 0) = LaLoHe1.GetX
         XYZ1(i, 1) = LaLoHe1.GetY
         XYZ1(i, 2) = LaLoHe1.GetZ
         LaLoHe2.SetLat (Lat2)
         LaLoHe2.SetLon (Lon2)
         LaLoHe2.SetHei (H2)
         dA = 6378137 ' wgs84 spheroid params
         dE2 = 0.00669437999014 'eccentr. squared
         LaLoHe2.LLH_to_XYZ dA, dE2
         XYZ2(i, 0) = LaLoHe2.GetX
         XYZ2(i, 1) = LaLoHe2.GetY
         XYZ2(i, 2) = LaLoHe2.GetZ
         XYZShifts(i, 0) = LaLoHe2.GetX - LaLoHe1.GetX
         XYZShifts(i, 1) = LaLoHe2.GetY - LaLoHe1.GetY
         XYZShifts(i, 2) = LaLoHe2.GetZ - LaLoHe1.GetZ
         DX = (DX * i + XYZShifts(i, 0)) / (i + 1)
         DY = (DY * i + XYZShifts(i, 1)) / (i + 1)
         DZ = (DZ * i + XYZShifts(i, 2)) / (i + 1)
         Set pFeature1 = pPointCursor1.NextFeature
        Set pFeature2 = pPointCursor2.NextFeature
        i = i + 1
    Loop
    MsgBox "Molodensky shifts calculated: " & vbCrLf & _
         "*    DX = " & DX & vbCrLf & _
         "*    DY = " & DY & vbCrLf & _
         "*    DZ = " & DZ & vbCrLf, vbInformation
End Sub

VERSION 1.0 CLASS
BEGIN
  MultiUse = -1  'True
END
Attribute VB_Name = "LatLonHei"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = False
Attribute VB_Exposed = False
Option Explicit
Private Lat As Double
Private Lon As Double
Private Hei As Double
Private ell_a As Double
Private ell_f_inv As Double
Private ell_e2 As Double
Private X As Double
Private Y As Double
Private Z As Double
Private Const M_PI = 3.14159265358979
Private Sub Class_Initialize()
  Lat = 0
  Lon = 0
  Hei = 0
  ell_a = 6378137
  ell_f_inv = 298.257223563
  ell_e2 = 6694380 / 1000000000
End Sub

Public Function GetLat() As Double
  GetLat = Lat
End Function

Public Sub SetLat(LatIn As Double)
  Lat = LatIn
End Sub
Public Function GetLon() As Double
  GetLon = Lon
End Function

Public Sub SetLon(LonIn As Double)
  Lon = LonIn
End Sub
Public Function GetHei() As Double
  GetHei = Hei
End Function

Public Sub SetHei(HeiIn As Double)
  Hei = HeiIn
End Sub
Public Function GetX() As Double
  GetX = X
End Function

Public Sub SetX(XIn As Double)
  X = XIn
End Sub
Public Function GetY() As Double
  GetY = Y
End Function

Public Sub SetY(YIn As Double)
  Y = YIn
End Sub
Public Function GetZ() As Double
  GetZ = Z
End Function

Public Sub SetZ(ZIn As Double)
  Z = ZIn
End Sub
Public Sub SetEllipsoidA(a As Double)
  ell_a = a
End Sub
Public Sub SetEllipsoidE2(e2 As Double)
  ell_e2 = e2
End Sub

Public Sub LLH_to_XYZ(ByVal a As Double, ByVal e2 As Double)
  Dim sinLat As Double
  Dim cosLat As Double
  Dim sinLon As Double
  Dim cosLon As Double
  sinLat = Sin(Lat * M_PI / 180)
  cosLat = Cos(Lat * M_PI / 180)
  sinLon = Sin(Lon * M_PI / 180)
  cosLon = Cos(Lon * M_PI / 180)
  Dim N As Double
  Dim D As Double
  D = 1 - e2 * sinLat * sinLat
  D = Sqr(D)
  N = a / Sqr(1 - e2 * sinLat * sinLat)
  X = (N + Hei) * cosLat * cosLon
  Y = (N + Hei) * cosLat * sinLon
  Z = (N * (1 - e2) + Hei) * sinLat
End Sub