Ngn som har lite kod (helst i vb.net) med en funktion som tar RT90-koordinater som inparametrar och ger WGS84decimalt som utdata?Koordinatkonvertering RT90>>WGS84 (dec)
Sv: Koordinatkonvertering RT90>>WGS84 (dec)
Public MustInherit Class Position
Public Property Latitude As Double
Public Property Longitude As Double
Public Property GridFormat As Grid
Public Sub New(ByVal lat As Double, ByVal lon As Double, ByVal format As Grid)
Latitude = lat
Longitude = lon
GridFormat = format
End Sub
Public Sub New(ByVal format As Grid)
GridFormat = format
End Sub
End Class
Public Class WGS84Position
Inherits Position
Public Enum WGS84Format
Degrees = 0
DegreesMinutes = 1
DegreesMinutesSeconds = 2
End Enum
Public Sub New()
MyBase.New(Grid.WGS84)
End Sub
Public Sub New(ByVal latitude As Double, ByVal longitude As Double)
MyBase.New(Grid.WGS84)
latitude = latitude
longitude = longitude
End Sub
Public Sub New(ByVal positionString As String, ByVal format As WGS84Format)
MyBase.New(Grid.WGS84)
If format = WGS84Format.Degrees Then
positionString = positionString.Trim()
Dim lat_lon As String() = positionString.Split(" "c)
If lat_lon.Length = 2 Then
Latitude = Double.Parse(lat_lon(0).Replace(".", ","))
Longitude = Double.Parse(lat_lon(1).Replace(".", ","))
Else
Throw New FormatException("The position string is invalid")
End If
ElseIf format = WGS84Format.DegreesMinutes OrElse format = WGS84Format.DegreesMinutesSeconds Then
Dim firstValueEndPos As Integer = 0
If format = WGS84Format.DegreesMinutes Then
firstValueEndPos = positionString.IndexOf("'")
ElseIf format = WGS84Format.DegreesMinutesSeconds Then
firstValueEndPos = positionString.IndexOf("""")
End If
Dim lat As String = positionString.Substring(0, firstValueEndPos + 1).Trim()
Dim lon As String = positionString.Substring(firstValueEndPos + 1).Trim()
SetLatitudeFromString(lat, format)
SetLongitudeFromString(lon, format)
End If
End Sub
Public Sub SetLatitudeFromString(ByVal value As String, ByVal format As WGS84Format)
value = value.Trim()
If format = WGS84Format.DegreesMinutes Then
Latitude = ParseValueFromDmString(value, "S")
ElseIf format = WGS84Format.DegreesMinutesSeconds Then
Latitude = ParseValueFromDmsString(value, "S")
ElseIf format = WGS84Format.Degrees Then
Latitude = Double.Parse(value)
End If
End Sub
Public Sub SetLongitudeFromString(ByVal value As String, ByVal format As WGS84Format)
If format = WGS84Format.DegreesMinutes Then
Longitude = ParseValueFromDmString(value, "W")
ElseIf format = WGS84Format.DegreesMinutesSeconds Then
Longitude = ParseValueFromDmsString(value, "W")
ElseIf format = WGS84Format.Degrees Then
Longitude = Double.Parse(value)
End If
End Sub
Public Function LatitudeToString(ByVal format As WGS84Format) As String
If format = WGS84Format.DegreesMinutes Then
Return ConvToDmString(Latitude, "N"c, "S"c)
ElseIf format = WGS84Format.DegreesMinutesSeconds Then
Return ConvToDmsString(Latitude, "N"c, "S"c)
Else
Return Latitude.ToString()
End If
End Function
Public Function LongitudeToString(ByVal format As WGS84Format) As String
If format = WGS84Format.DegreesMinutes Then
Return ConvToDmString(Longitude, "E"c, "W"c)
ElseIf format = WGS84Format.DegreesMinutesSeconds Then
Return ConvToDmsString(Longitude, "E"c, "W"c)
Else
Return Longitude.ToString()
End If
End Function
Private Function ConvToDmString(ByVal value As Double, ByVal positiveValue As Char, ByVal negativeValue As Char) As String
If value = Double.MinValue Then
Return ""
End If
Dim degrees = Math.Floor(Math.Abs(value))
Dim minutes = (Math.Abs(value) - degrees) * 60
Return String.Format("{0} {1}º {2}'", If(value >= 0, positiveValue, negativeValue), degrees, (Math.Floor(minutes * 10000) / 10000))
End Function
Private Function ConvToDmsString(ByVal value As Double, ByVal positiveValue As Char, ByVal negativeValue As Char) As String
If value = Double.MinValue Then
Return ""
End If
Dim degrees = Math.Floor(Math.Abs(value))
Dim minutes = Math.Floor((Math.Abs(value) - degrees) * 60)
Dim seconds = (Math.Abs(value) - degrees - minutes / 60) * 3600
Return String.Format("{0} {1}º {2}' {3}""", If(value >= 0, positiveValue, negativeValue), degrees, minutes, Math.Round(seconds, 5))
End Function
Private Function ParseValueFromDmString(ByVal value As String, ByVal positiveChar As String) As Double
Dim retVal As Double = 0.0
If Not String.IsNullOrEmpty(value) Then
Dim direction As String = value(0).ToString()
value = value.Substring(1).Trim()
Dim degree As String = value.Substring(0, value.IndexOf("º"))
value = value.Substring(value.IndexOf("º") + 1).Trim()
Dim minutes As String = value.Substring(0, value.IndexOf("'"))
value = value.Substring(value.IndexOf("'") + 1).Trim()
retVal = Double.Parse(degree)
retVal += Double.Parse(minutes.Replace(".", ",")) / 60
If retVal > 90 Then
retVal = Double.MinValue
End If
If direction = positiveChar OrElse direction = "-" Then
' Cannot convert ExpressionStatementSyntax, System.ArgumentOutOfRangeException: Exception of type 'System.ArgumentOutOfRangeException' was thrown.
' Parameter name: op
' Actual value was MultiplyAssignmentStatement.
' at ICSharpCode.CodeConverter.Util.VBUtil.GetExpressionOperatorTokenKind(SyntaxKind op)
' at ICSharpCode.CodeConverter.VB.NodesVisitor.MakeAssignmentStatement(AssignmentExpressionSyntax node)
' at ICSharpCode.CodeConverter.VB.NodesVisitor.VisitAssignmentExpression(AssignmentExpressionSyntax node)
' at Microsoft.CodeAnalysis.CSharp.Syntax.AssignmentExpressionSyntax.Accept[TResult](CSharpSyntaxVisitor`1 visitor)
' at Microsoft.CodeAnalysis.CSharp.CSharpSyntaxVisitor`1.Visit(SyntaxNode node)
' at ICSharpCode.CodeConverter.VB.CommentConvertingNodesVisitor.DefaultVisit(SyntaxNode node)
' at Microsoft.CodeAnalysis.CSharp.CSharpSyntaxVisitor`1.VisitAssignmentExpression(AssignmentExpressionSyntax node)
' at Microsoft.CodeAnalysis.CSharp.Syntax.AssignmentExpressionSyntax.Accept[TResult](CSharpSyntaxVisitor`1 visitor)
' at ICSharpCode.CodeConverter.VB.MethodBodyVisitor.ConvertSingleExpression(ExpressionSyntax node)
' at ICSharpCode.CodeConverter.VB.MethodBodyVisitor.VisitExpressionStatement(ExpressionStatementSyntax node)
' at Microsoft.CodeAnalysis.CSharp.Syntax.ExpressionStatementSyntax.Accept[TResult](CSharpSyntaxVisitor`1 visitor)
' at Microsoft.CodeAnalysis.CSharp.CSharpSyntaxVisitor`1.Visit(SyntaxNode node)
' at ICSharpCode.CodeConverter.VB.CommentConvertingMethodBodyVisitor.ConvertWithTrivia(SyntaxNode node)
' at ICSharpCode.CodeConverter.VB.CommentConvertingMethodBodyVisitor.DefaultVisit(SyntaxNode node)
'
' Input:
' retVal *= -1;
'
End If
Else
retVal = Double.MinValue
End If
Return retVal
End Function
Private Function ParseValueFromDmsString(ByVal value As String, ByVal positiveChar As String) As Double
Dim retVal As Double = 0.0
If Not String.IsNullOrEmpty(value) Then
Dim direction As String = value(0).ToString()
value = value.Substring(1).Trim()
Dim degree As String = value.Substring(0, value.IndexOf("º"))
value = value.Substring(value.IndexOf("º") + 1).Trim()
Dim minutes As String = value.Substring(0, value.IndexOf("'"))
value = value.Substring(value.IndexOf("'") + 1).Trim()
Dim seconds As String = value.Substring(0, value.IndexOf(""""))
retVal = Double.Parse(degree)
retVal += Double.Parse(minutes) / 60
retVal += Double.Parse(seconds.Replace(".", ",")) / 3600
If retVal > 90 Then
retVal = Double.MinValue
Return retVal
End If
If direction = positiveChar OrElse direction = "-" Then
' Cannot convert ExpressionStatementSyntax, System.ArgumentOutOfRangeException: Exception of type 'System.ArgumentOutOfRangeException' was thrown.
' Parameter name: op
' Actual value was MultiplyAssignmentStatement.
' at ICSharpCode.CodeConverter.Util.VBUtil.GetExpressionOperatorTokenKind(SyntaxKind op)
' at ICSharpCode.CodeConverter.VB.NodesVisitor.MakeAssignmentStatement(AssignmentExpressionSyntax node)
' at ICSharpCode.CodeConverter.VB.NodesVisitor.VisitAssignmentExpression(AssignmentExpressionSyntax node)
' at Microsoft.CodeAnalysis.CSharp.Syntax.AssignmentExpressionSyntax.Accept[TResult](CSharpSyntaxVisitor`1 visitor)
' at Microsoft.CodeAnalysis.CSharp.CSharpSyntaxVisitor`1.Visit(SyntaxNode node)
' at ICSharpCode.CodeConverter.VB.CommentConvertingNodesVisitor.DefaultVisit(SyntaxNode node)
' at Microsoft.CodeAnalysis.CSharp.CSharpSyntaxVisitor`1.VisitAssignmentExpression(AssignmentExpressionSyntax node)
' at Microsoft.CodeAnalysis.CSharp.Syntax.AssignmentExpressionSyntax.Accept[TResult](CSharpSyntaxVisitor`1 visitor)
' at ICSharpCode.CodeConverter.VB.MethodBodyVisitor.ConvertSingleExpression(ExpressionSyntax node)
' at ICSharpCode.CodeConverter.VB.MethodBodyVisitor.VisitExpressionStatement(ExpressionStatementSyntax node)
' at Microsoft.CodeAnalysis.CSharp.Syntax.ExpressionStatementSyntax.Accept[TResult](CSharpSyntaxVisitor`1 visitor)
' at Microsoft.CodeAnalysis.CSharp.CSharpSyntaxVisitor`1.Visit(SyntaxNode node)
' at ICSharpCode.CodeConverter.VB.CommentConvertingMethodBodyVisitor.ConvertWithTrivia(SyntaxNode node)
' at ICSharpCode.CodeConverter.VB.CommentConvertingMethodBodyVisitor.DefaultVisit(SyntaxNode node)
'
' Input:
' retVal *= -1;
'
End If
Else
retVal = Double.MinValue
End If
Return retVal
End Function
Public Overrides Function ToString() As String
Return String.Format("Latitude: {0} Longitude: {1}", LatitudeToString(WGS84Format.DegreesMinutesSeconds), LongitudeToString(WGS84Format.DegreesMinutesSeconds))
End Function
End Class
Public Class RT90Position
Inherits Position
Public Enum RT90Projection
rt90_7_5_gon_v = 0
rt90_5_0_gon_v = 1
rt90_2_5_gon_v = 2
rt90_0_0_gon_v = 3
rt90_2_5_gon_o = 4
rt90_5_0_gon_o = 5
End Enum
Public Sub New(ByVal x As Double, ByVal y As Double)
MyBase.New(x, y, Grid.RT90)
Projection = RT90Projection.rt90_2_5_gon_v
End Sub
Public Sub New(ByVal x As Double, ByVal y As Double, ByVal projection As RT90Projection)
MyBase.New(x, y, Grid.RT90)
projection = projection
End Sub
Public Sub New(ByVal position As WGS84Position, ByVal rt90projection As RT90Projection)
MyBase.New(Grid.RT90)
Dim gkProjection As GaussKreuger = New GaussKreuger()
gkProjection.swedish_params(GetProjectionString(rt90projection))
Dim lat_lon = gkProjection.geodetic_to_grid(position.Latitude, position.Longitude)
Latitude = lat_lon(0)
Longitude = lat_lon(1)
Projection = rt90projection
End Sub
Public Function ToWGS84() As WGS84Position
Dim gkProjection As GaussKreuger = New GaussKreuger()
gkProjection.swedish_params(ProjectionString)
Dim lat_lon = gkProjection.grid_to_geodetic(Latitude, Longitude)
Dim newPos As WGS84Position = New WGS84Position() With {
.Latitude = lat_lon(0),
.Longitude = lat_lon(1),
.GridFormat = Grid.WGS84
}
Return newPos
End Function
Private Function GetProjectionString(ByVal projection As RT90Projection) As String
Dim retVal As String = String.Empty
Select Case projection
Case RT90Projection.rt90_7_5_gon_v
retVal = "rt90_7.5_gon_v"
Case RT90Projection.rt90_5_0_gon_v
retVal = "rt90_5.0_gon_v"
Case RT90Projection.rt90_2_5_gon_v
retVal = "rt90_2.5_gon_v"
Case RT90Projection.rt90_0_0_gon_v
retVal = "rt90_0.0_gon_v"
Case RT90Projection.rt90_2_5_gon_o
retVal = "rt90_2.5_gon_o"
Case RT90Projection.rt90_5_0_gon_o
retVal = "rt90_5.0_gon_o"
Case Else
retVal = "rt90_2.5_gon_v"
End Select
Return retVal
End Function
Public Property Projection As RT90Projection
Public ReadOnly Property ProjectionString As String
Get
Return GetProjectionString(Projection)
End Get
End Property
Public Overrides Function ToString() As String
Return String.Format("X: {0} Y: {1} Projection: {2}", Latitude, Longitude, ProjectionString)
End Function
End Class
Public Class SWEREF99Position
Inherits Position
Public Enum SWEREFProjection
sweref_99_tm = 0
sweref_99_12_00 = 1
sweref_99_13_30 = 2
sweref_99_15_00 = 3
sweref_99_16_30 = 4
sweref_99_18_00 = 5
sweref_99_14_15 = 6
sweref_99_15_45 = 7
sweref_99_17_15 = 8
sweref_99_18_45 = 9
sweref_99_20_15 = 10
sweref_99_21_45 = 11
sweref_99_23_15 = 12
End Enum
Public Sub New(ByVal n As Double, ByVal e As Double)
MyBase.New(n, e, Grid.SWEREF99)
Projection = SWEREFProjection.sweref_99_tm
End Sub
Public Sub New(ByVal n As Double, ByVal e As Double, ByVal projection As SWEREFProjection)
MyBase.New(n, e, Grid.SWEREF99)
projection = projection
End Sub
Public Sub New(ByVal position As WGS84Position, ByVal projection As SWEREFProjection)
MyBase.New(Grid.SWEREF99)
Dim gkProjection As GaussKreuger = New GaussKreuger()
gkProjection.swedish_params(GetProjectionString(projection))
Dim lat_lon = gkProjection.geodetic_to_grid(position.Latitude, position.Longitude)
Latitude = lat_lon(0)
Longitude = lat_lon(1)
projection = projection
End Sub
Public Function ToWGS84() As WGS84Position
Dim gkProjection As GaussKreuger = New GaussKreuger()
gkProjection.swedish_params(ProjectionString)
Dim lat_lon = gkProjection.grid_to_geodetic(Latitude, Longitude)
Dim newPos As WGS84Position = New WGS84Position() With {
.Latitude = lat_lon(0),
.Longitude = lat_lon(1),
.GridFormat = Grid.WGS84
}
Return newPos
End Function
Private Function GetProjectionString(ByVal projection As SWEREFProjection) As String
Dim retVal As String = String.Empty
Select Case projection
Case SWEREFProjection.sweref_99_tm
retVal = "sweref_99_tm"
Case SWEREFProjection.sweref_99_12_00
retVal = "sweref_99_1200"
Case SWEREFProjection.sweref_99_13_30
retVal = "sweref_99_1330"
Case SWEREFProjection.sweref_99_14_15
retVal = "sweref_99_1415"
Case SWEREFProjection.sweref_99_15_00
retVal = "sweref_99_1500"
Case SWEREFProjection.sweref_99_15_45
retVal = "sweref_99_1545"
Case SWEREFProjection.sweref_99_16_30
retVal = "sweref_99_1630"
Case SWEREFProjection.sweref_99_17_15
retVal = "sweref_99_1715"
Case SWEREFProjection.sweref_99_18_00
retVal = "sweref_99_1800"
Case SWEREFProjection.sweref_99_18_45
retVal = "sweref_99_1845"
Case SWEREFProjection.sweref_99_20_15
retVal = "sweref_99_2015"
Case SWEREFProjection.sweref_99_21_45
retVal = "sweref_99_2145"
Case SWEREFProjection.sweref_99_23_15
retVal = "sweref_99_2315"
Case Else
retVal = "sweref_99_tm"
End Select
Return retVal
End Function
Public Property Projection As SWEREFProjection
Public ReadOnly Property ProjectionString As String
Get
Return GetProjectionString(Projection)
End Get
End Property
Public Overrides Function ToString() As String
Return String.Format("N: {0} E: {1} Projection: {2}", Latitude, Longitude, ProjectionString)
End Function
End Class
Public Class GaussKreuger
Private axis As Double
Private flattening As Double
Private central_meridian As Double
Private scale As Double
Private false_northing As Double
Private false_easting As Double
Public Sub swedish_params(ByVal projection As String)
If projection = "rt90_7.5_gon_v" Then
grs80_params()
central_meridian = 11.0 + 18.375 / 60.0
scale = 1.000006
false_northing = -667.282
false_easting = 1500025.141
ElseIf projection = "rt90_5.0_gon_v" Then
grs80_params()
central_meridian = 13.0 + 33.376 / 60.0
scale = 1.0000058
false_northing = -667.13
false_easting = 1500044.695
ElseIf projection = "rt90_2.5_gon_v" Then
grs80_params()
central_meridian = 15.0 + 48.0 / 60.0 + 22.624306 / 3600.0
scale = 1.00000561024
false_northing = -667.711
false_easting = 1500064.274
ElseIf projection = "rt90_0.0_gon_v" Then
grs80_params()
central_meridian = 18.0 + 3.378 / 60.0
scale = 1.0000054
false_northing = -668.844
false_easting = 1500083.521
ElseIf projection = "rt90_2.5_gon_o" Then
grs80_params()
central_meridian = 20.0 + 18.379 / 60.0
scale = 1.0000052
false_northing = -670.706
false_easting = 1500102.765
ElseIf projection = "rt90_5.0_gon_o" Then
grs80_params()
central_meridian = 22.0 + 33.38 / 60.0
scale = 1.0000049
false_northing = -672.557
false_easting = 1500121.846
ElseIf projection = "bessel_rt90_7.5_gon_v" Then
bessel_params()
central_meridian = 11.0 + 18.0 / 60.0 + 29.8 / 3600.0
ElseIf projection = "bessel_rt90_5.0_gon_v" Then
bessel_params()
central_meridian = 13.0 + 33.0 / 60.0 + 29.8 / 3600.0
ElseIf projection = "bessel_rt90_2.5_gon_v" Then
bessel_params()
central_meridian = 15.0 + 48.0 / 60.0 + 29.8 / 3600.0
ElseIf projection = "bessel_rt90_0.0_gon_v" Then
bessel_params()
central_meridian = 18.0 + 3.0 / 60.0 + 29.8 / 3600.0
ElseIf projection = "bessel_rt90_2.5_gon_o" Then
bessel_params()
central_meridian = 20.0 + 18.0 / 60.0 + 29.8 / 3600.0
ElseIf projection = "bessel_rt90_5.0_gon_o" Then
bessel_params()
central_meridian = 22.0 + 33.0 / 60.0 + 29.8 / 3600.0
ElseIf projection = "sweref_99_tm" Then
sweref99_params()
central_meridian = 15.0
scale = 0.9996
false_northing = 0.0
false_easting = 500000.0
ElseIf projection = "sweref_99_1200" Then
sweref99_params()
central_meridian = 12.0
ElseIf projection = "sweref_99_1330" Then
sweref99_params()
central_meridian = 13.5
ElseIf projection = "sweref_99_1500" Then
sweref99_params()
central_meridian = 15.0
ElseIf projection = "sweref_99_1630" Then
sweref99_params()
central_meridian = 16.5
ElseIf projection = "sweref_99_1800" Then
sweref99_params()
central_meridian = 18.0
ElseIf projection = "sweref_99_1415" Then
sweref99_params()
central_meridian = 14.25
ElseIf projection = "sweref_99_1545" Then
sweref99_params()
central_meridian = 15.75
ElseIf projection = "sweref_99_1715" Then
sweref99_params()
central_meridian = 17.25
ElseIf projection = "sweref_99_1845" Then
sweref99_params()
central_meridian = 18.75
ElseIf projection = "sweref_99_2015" Then
sweref99_params()
central_meridian = 20.25
ElseIf projection = "sweref_99_2145" Then
sweref99_params()
central_meridian = 21.75
ElseIf projection = "sweref_99_2315" Then
sweref99_params()
central_meridian = 23.25
Else
central_meridian = Double.MinValue
End If
End Sub
Private Sub grs80_params()
axis = 6378137.0
flattening = 1.0 / 298.257222101
central_meridian = Double.MinValue
End Sub
Private Sub bessel_params()
axis = 6377397.155
flattening = 1.0 / 299.1528128
central_meridian = Double.MinValue
scale = 1.0
false_northing = 0.0
false_easting = 1500000.0
End Sub
Private Sub sweref99_params()
axis = 6378137.0
flattening = 1.0 / 298.257222101
central_meridian = Double.MinValue
scale = 1.0
false_northing = 0.0
false_easting = 150000.0
End Sub
Public Function geodetic_to_grid(ByVal latitude As Double, ByVal longitude As Double) As Double()
Dim x_y As Double() = New Double(1) {}
Dim e2 As Double = flattening * (2.0 - flattening)
Dim n As Double = flattening / (2.0 - flattening)
Dim a_roof As Double = axis / (1.0 + n) * (1.0 + n * n / 4.0 + n * n * n * n / 64.0)
Dim A As Double = e2
Dim B As Double = (5.0 * e2 * e2 - e2 * e2 * e2) / 6.0
Dim C As Double = (104.0 * e2 * e2 * e2 - 45.0 * e2 * e2 * e2 * e2) / 120.0
Dim D As Double = (1237.0 * e2 * e2 * e2 * e2) / 1260.0
Dim beta1 As Double = n / 2.0 - 2.0 * n * n / 3.0 + 5.0 * n * n * n / 16.0 + 41.0 * n * n * n * n / 180.0
Dim beta2 As Double = 13.0 * n * n / 48.0 - 3.0 * n * n * n / 5.0 + 557.0 * n * n * n * n / 1440.0
Dim beta3 As Double = 61.0 * n * n * n / 240.0 - 103.0 * n * n * n * n / 140.0
Dim beta4 As Double = 49561.0 * n * n * n * n / 161280.0
Dim deg_to_rad As Double = Math.PI / 180.0
Dim phi As Double = latitude * deg_to_rad
Dim lambda As Double = longitude * deg_to_rad
Dim lambda_zero As Double = central_meridian * deg_to_rad
Dim phi_star As Double = phi - Math.Sin(phi) * Math.Cos(phi) * (A + B * Math.Pow(Math.Sin(phi), 2) + C * Math.Pow(Math.Sin(phi), 4) + D * Math.Pow(Math.Sin(phi), 6))
Dim delta_lambda As Double = lambda - lambda_zero
Dim xi_prim As Double = Math.Atan(Math.Tan(phi_star) / Math.Cos(delta_lambda))
Dim eta_prim As Double = math_atanh(Math.Cos(phi_star) * Math.Sin(delta_lambda))
Dim x As Double = scale * a_roof * (xi_prim + beta1 * Math.Sin(2.0 * xi_prim) * math_cosh(2.0 * eta_prim) + beta2 * Math.Sin(4.0 * xi_prim) * math_cosh(4.0 * eta_prim) + beta3 * Math.Sin(6.0 * xi_prim) * math_cosh(6.0 * eta_prim) + beta4 * Math.Sin(8.0 * xi_prim) * math_cosh(8.0 * eta_prim)) + false_northing
Dim y As Double = scale * a_roof * (eta_prim + beta1 * Math.Cos(2.0 * xi_prim) * math_sinh(2.0 * eta_prim) + beta2 * Math.Cos(4.0 * xi_prim) * math_sinh(4.0 * eta_prim) + beta3 * Math.Cos(6.0 * xi_prim) * math_sinh(6.0 * eta_prim) + beta4 * Math.Cos(8.0 * xi_prim) * math_sinh(8.0 * eta_prim)) + false_easting
x_y(0) = Math.Round(x * 1000.0) / 1000.0
x_y(1) = Math.Round(y * 1000.0) / 1000.0
Return x_y
End Function
Public Function grid_to_geodetic(ByVal x As Double, ByVal y As Double) As Double()
Dim lat_lon As Double() = New Double(1) {}
If central_meridian = Double.MinValue Then
Return lat_lon
End If
Dim e2 As Double = flattening * (2.0 - flattening)
Dim n As Double = flattening / (2.0 - flattening)
Dim a_roof As Double = axis / (1.0 + n) * (1.0 + n * n / 4.0 + n * n * n * n / 64.0)
Dim delta1 As Double = n / 2.0 - 2.0 * n * n / 3.0 + 37.0 * n * n * n / 96.0 - n * n * n * n / 360.0
Dim delta2 As Double = n * n / 48.0 + n * n * n / 15.0 - 437.0 * n * n * n * n / 1440.0
Dim delta3 As Double = 17.0 * n * n * n / 480.0 - 37 * n * n * n * n / 840.0
Dim delta4 As Double = 4397.0 * n * n * n * n / 161280.0
Dim Astar As Double = e2 + e2 * e2 + e2 * e2 * e2 + e2 * e2 * e2 * e2
Dim Bstar As Double = -(7.0 * e2 * e2 + 17.0 * e2 * e2 * e2 + 30.0 * e2 * e2 * e2 * e2) / 6.0
Dim Cstar As Double = (224.0 * e2 * e2 * e2 + 889.0 * e2 * e2 * e2 * e2) / 120.0
Dim Dstar As Double = -(4279.0 * e2 * e2 * e2 * e2) / 1260.0
Dim deg_to_rad As Double = Math.PI / 180
Dim lambda_zero As Double = central_meridian * deg_to_rad
Dim xi As Double = (x - false_northing) / (scale * a_roof)
Dim eta As Double = (y - false_easting) / (scale * a_roof)
Dim xi_prim As Double = xi - delta1 * Math.Sin(2.0 * xi) * math_cosh(2.0 * eta) - delta2 * Math.Sin(4.0 * xi) * math_cosh(4.0 * eta) - delta3 * Math.Sin(6.0 * xi) * math_cosh(6.0 * eta) - delta4 * Math.Sin(8.0 * xi) * math_cosh(8.0 * eta)
Dim eta_prim As Double = eta - delta1 * Math.Cos(2.0 * xi) * math_sinh(2.0 * eta) - delta2 * Math.Cos(4.0 * xi) * math_sinh(4.0 * eta) - delta3 * Math.Cos(6.0 * xi) * math_sinh(6.0 * eta) - delta4 * Math.Cos(8.0 * xi) * math_sinh(8.0 * eta)
Dim phi_star As Double = Math.Asin(Math.Sin(xi_prim) / math_cosh(eta_prim))
Dim delta_lambda As Double = Math.Atan(math_sinh(eta_prim) / Math.Cos(xi_prim))
Dim lon_radian As Double = lambda_zero + delta_lambda
Dim lat_radian As Double = phi_star + Math.Sin(phi_star) * Math.Cos(phi_star) * (Astar + Bstar * Math.Pow(Math.Sin(phi_star), 2) + Cstar * Math.Pow(Math.Sin(phi_star), 4) + Dstar * Math.Pow(Math.Sin(phi_star), 6))
lat_lon(0) = lat_radian * 180.0 / Math.PI
lat_lon(1) = lon_radian * 180.0 / Math.PI
Return lat_lon
End Function
Private Function math_sinh(ByVal value As Double) As Double
Return 0.5 * (Math.Exp(value) - Math.Exp(-value))
End Function
Private Function math_cosh(ByVal value As Double) As Double
Return 0.5 * (Math.Exp(value) + Math.Exp(-value))
End Function
Private Function math_atanh(ByVal value As Double) As Double
Return 0.5 * Math.Log((1.0 + value) / (1.0 - value))
End Function
End Class