地理座標系変換による座標変換

'地理座標系変換
Public Sub Project()

Dim pMxDocument As IMxDocument
Set pMxDocument = ThisDocument

Dim pGCS_Tokyo As ISpatialReference
Set pGCS_Tokyo = fncSpatialReference(esriSRGeoCSType.esriSRGeoCS_Tokyo, True)

Dim pGCS_JGD_2000 As ISpatialReference
Set pGCS_JGD_2000 = fncSpatialReference(esriSRGeoCS3Type.esriSRGeoCS_JapanGeodeticDatum2000, True)


Dim pGCS_JGD_2011 As ISpatialReference
Set pGCS_JGD_2011 = fncSpatialReference(6668, True)

Dim pPoint As IPoint
Set pPoint = New Point
Dim pGeometry As IGeometry2
Set pGeometry = pPoint
pPoint.PutCoords 135, 35
Set pPoint.SpatialReference = pGCS_Tokyo

Dim dd As Long
Dim mm As Long
Dim ss As Double


'元座標
DD2DMS pPoint.X, dd, mm, ss
Debug.Print "Tokyo E", dd, mm, ss

DD2DMS pPoint.Y, dd, mm, ss
Debug.Print "Tokyo N", dd, mm, ss


Dim pTKY2JGD2000_NTv2 As IGeoTransformation
Set pTKY2JGD2000_NTv2 = CreateTKY2JGD2000_NTv2Transformation


Dim pTKY2JGD2000_1 As IGeoTransformation
Set pTKY2JGD2000_1 = CreateGeoTransformationTKY2JGD2000

Dim pJGD20002JGD2011_NTv2 As IGeoTransformation
Set pJGD20002JGD2011_NTv2 = CreateTKY2JGD2011_NTv2Transformation

'計測
Dim time As Double
time = Timer

'地理座標系変換
'pGeometry.ProjectEx pGCS_JGD_2000, esriTransformForward, pTKY2JGD2000_NTv2, True, 0, 0
'Geometry.ProjectEx pGCS_JGD_2000, esriTransformForward, pTKY2JGD2000_1, True, 0, 0
pGeometry.ProjectEx pGCS_JGD_2011, esriTransformForward, pJGD20002JGD2011_NTv2, True, 0, 0

'計測時間
Dim time2 As Double
time2 = Timer - time

DD2DMS pPoint.X, dd, mm, ss
Debug.Print "JGD2011 E", dd, mm, ss

DD2DMS pPoint.Y, dd, mm, ss
Debug.Print "JGD2011 N", dd, mm, ss

Debug.Print "処理時間:" & time2


End Sub

Function fncSpatialReference(valType As Long, blnGeographic As Boolean) As ISpatialReference

'DeveloperHelpから以下を参照

'地理座標系(pcsType)
'esriSRGeoCSType Contains
'esriSRGeoCS2Type Constants
'esriSRGeoCS3Type Constants

'投影座標系(gcsType)
'esriSRProjCSType Constants
'esriSRProjCS2Type Constants
'esriSRProjCS3Type Constants
'esriSRProjCS4Type Constants

'例
'valType = esriSRProjCS4Type.esriSRProjCS_JGD2000_Japan_Zone_9 'JGD2000 平面直角座標系第9系
'valType = esriSRGeoCS3Type.esriSRGeoCS_JapanGeodeticDatum2000 'JGD2000 地理座標系

Dim pSpatialReferenceFactory As ISpatialReferenceFactory
Set pSpatialReferenceFactory = New SpatialReferenceEnvironment

If blnGeographic = True Then
Dim pGeographicCoordinateSystem As IGeographicCoordinateSystem
Set pGeographicCoordinateSystem = pSpatialReferenceFactory.CreateGeographicCoordinateSystem(valType)

Set fncSpatialReference = pGeographicCoordinateSystem

Else
Dim pProjectedCoordinateSystem As IProjectedCoordinateSystem
Set pProjectedCoordinateSystem = pSpatialReferenceFactory.CreateProjectedCoordinateSystem(valType)

Set fncSpatialReference = pProjectedCoordinateSystem
End If

'XY座標精度の設定（9.2以降は設定必須）
Dim pSpatialReferenceResolution As ISpatialReferenceResolution
Set pSpatialReferenceResolution = fncSpatialReference
pSpatialReferenceResolution.SetDefaultXYResolution

'XY許容値の設定（9.2以降は設定必須）
Dim pSpatialReferenceTolerance As ISpatialReferenceTolerance
Set pSpatialReferenceTolerance = fncSpatialReference
pSpatialReferenceTolerance.SetDefaultXYTolerance

End Function


'NTv2Transformation(TKY2JGD)の作成
Public Function CreateTKY2JGD2000_NTv2Transformation() As IGeoTransformation
Dim pGridTransformation As IGridTransformation
Set pGridTransformation = New NTv2Transformation

pGridTransformation.Name = "Tokyo_To_JGD_2000_NTv2"
pGridTransformation.GridDatasetName = "Dataset_japan/tky2jgd.gsb"
pGridTransformation.Load

Set CreateTKY2JGD2000_NTv2Transformation = pGridTransformation

End Function

'NTv2Transformation(JGD2000 to JGD2011)の作成
Public Function CreateTKY2JGD2011_NTv2Transformation() As IGeoTransformation

Dim pGridTransformation As IGridTransformation
Set pGridTransformation = New NTv2Transformation

pGridTransformation.Name = "JGD_2000_To_JGD_2011_NTv2"
pGridTransformation.GridDatasetName = "Dataset_japan/touhokutaiheiyouoki2011.gsb"
pGridTransformation.Load

Dim pCompositeGeoTransformation As ICompositeGeoTransformation
Set pCompositeGeoTransformation = New CompositeGeoTransformation
pCompositeGeoTransformation.Add esriTransformForward, CreateTKY2JGD2000_NTv2Transformation
pCompositeGeoTransformation.Add esriTransformForward, pGridTransformation


Set CreateTKY2JGD2011_NTv2Transformation = pCompositeGeoTransformation

End Function


'NTv2Transformation(JGD2000 → JGD2011 )の作成
Public Function CreateJGD20002JGD2011_NTv2Transformation() As IGeoTransformation
Dim pGridTransformation As IGridTransformation
Set pGridTransformation = New NTv2Transformation

pGridTransformation.Name = "JGD_2000_To_JGD_2011_NTv2"
pGridTransformation.GridDatasetName = "Dataset_japan/tky2jgd.gsb"
pGridTransformation.Load


Set CreateTKY2JGD_NTv2Transformation = pGridTransformation

End Function

'GeocentricTranslationの作成
'3 パラメータ変換（地心変換）
Public Function CreateGeoTransformationTKY2JGD2000() As IGeoTransformation
Dim pSpatialReferenceFactory As ISpatialReferenceFactory
Set pSpatialReferenceFactory = New SpatialReferenceEnvironment

'定数による作成
Dim pGeocentricTranslation As IGeocentricTranslation
Set pGeocentricTranslation = pSpatialReferenceFactory.CreateGeoTransformation(esriSRGeoTransformation3Type.esriSRGeoTransformation_Tokyo_To_JGD_2000_1)

' Set pGeocentricTranslation = New GeocentricTranslation
' pGeocentricTranslation.Name = "Tokyo_To_JGD_2000_1"
' pGeocentricTranslation.PutParameters -146.414, 507.337, 680.507

Set CreateGeoTransformationTKY2JGD2000 = pGeocentricTranslation

End Function


Private Function DD2DMS(ByRef DecimalDegree As Double, ByRef Degree As Long, ByRef Minute As Long, ByRef Second As Double, Optional Reverce As Boolean = False)
'十進経緯から度分秒に変換
If Reverce = False Then
Dim pLatLonFormat As ILatLonFormat
Set pLatLonFormat = New LatLonFormat
pLatLonFormat.GetDMS DecimalDegree, Degree, Minute, Second

'度分秒から十進経緯度に変換
Else
DecimalDegree = Degree + Minute / 60 + Second / 3600

End If

End Function