Imagine, digital earth maps now just like real earth. Just look at Google Earth or Microsoft Virtual Earth, satellite photo result can show free for us. It will be nice if we have mobile gadget with GPS, we will always know where are we. Technicaly Map in software world has to be trend begin when GIS needs increase. Some developer confuse with coordinates what is UTM or what is Degree. Same with me, i have problem to build GIS system and I use SVG or XAML to make a map, generated from ArcView or mapinfo. The booth software support any coordinates type, but our client use degrees format to put the location. SVG or XAML are not support using degress, they both just support using coordinates X,Y. So I make the function to convert Degrees to UTM or UTM to Degrees.
Public Class Ellipsoid
Public id As Integer
Public ellipsoidName As String
Public EquatorialRadius As Double
Public eccentricitySquared As Double
Const PI As Double = 3.14159265
Const FOURTHPI As Double = PI / 4
Const deg2rad As Double = PI / 180
Const rad2deg As Double = 180.0 / PI
Const oneMiles = 1609.344
Const oneFoot = 0.3048
Const oneDegreeLat = 111044.736
Const oneMinuteLat = 1850.7456
Const oneSecondLat = 30.84576
Const oneDegreeLong = 67592.448
Const oneMinuteLong = 1126.5408
Const oneSecondLong = 18.77568
Public Sub New(ByVal iid As Integer, ByVal sname As String, ByVal radius As Double, ByVal ecc As Double)
id = iid
ellipsoidName = sname
EquatorialRadius = radius
eccentricitySquared = ecc
End Sub
Public Sub New(ByVal iid As Integer)
SetID(iid)
End Sub
Sub GPS2XY(ByVal gpspos As String, ByRef x As Double, ByRef y As Double)
If gpspos <> "" Then
Dim latdeg As String = ""
Dim longdeg As String = ""
Dim UTMZone As String = ""
GPS2LATLONGDeg(gpspos, latdeg, longdeg)
Dim lat As Double = DegMinSec2Deg(latdeg)
Dim longi As Double = DegMinSec2Deg(longdeg)
LLtoUTM(lat, longi, x, y, UTMZone)
If InStr(gpspos, "S") > 0 Then y = -1 * Math.Abs(y)
'If InStr(gpspos, "E") > 0 Then x = -1 * Math.Abs(x)
Else
x = 0
y = 0
End If
End Sub
Sub GPS2LATLONGDeg(ByVal gpspos As String, ByRef latdeg As String, ByRef longdeg As String)
gpspos = gpspos.Replace(";", "")
gpspos = gpspos.Replace(":", "")
Dim charx As String = IIf(InStr(gpspos, "N") > 0, "N", "S")
Dim chary As String = IIf(InStr(gpspos, "E") > 0, "E", "W")
latdeg = Trim(Mid(gpspos, InStr(gpspos, charx) + 1, InStr(gpspos, chary) - 2))
longdeg = Trim(Mid(gpspos, InStr(gpspos, chary) + 1))
If InStr(gpspos, "S") > 0 Then latdeg = "-" & latdeg
If InStr(gpspos, "W") > 0 Then longdeg = "-" & longdeg
End Sub
Function DegMinSec2Deg(ByVal degree As String) As Double
Dim ddegree As Double = 0
Dim dminute As Double = 0
Dim dsecond As Double = 0
Dim ismin As Boolean = (InStr(degree, "-") > 0)
Dim lat As Double = 0
If InStr(UCase(degree), "D") > 0 Then
ddegree = CDbl(Trim(Mid(degree, 1, InStr(UCase(degree), "D") - 1)))
degree = Trim(Mid(degree, InStr(UCase(degree), "D") + 1))
ElseIf InStr(degree, Chr(176)) > 0 Then
ddegree = CDbl(Trim(Mid(degree, 1, InStr(UCase(degree), Chr(176)) - 1)))
degree = Trim(Mid(degree, InStr(UCase(degree), Chr(176)) + 1))
ElseIf InStr(degree, Chr(186)) > 0 Then
ddegree = CDbl(Trim(Mid(degree, 1, InStr(UCase(degree), Chr(186)) - 1)))
degree = Trim(Mid(degree, InStr(UCase(degree), Chr(186)) + 1))
Else
ddegree = 0
End If
If InStr(UCase(degree), "'") > 0 Then
dminute = CDbl(Trim(Mid(degree, 1, InStr(UCase(degree), "'") - 1)))
Else
dminute = 0
End If
degree = Trim(Mid(degree, InStr(UCase(degree), "'") + 1))
If InStr(UCase(degree), """") > 0 Then
dsecond = CDbl(Trim(Mid(degree, 1, InStr(UCase(degree), """") - 1)))
ElseIf InStr(UCase(degree), "''") > 0 Then
dsecond = CDbl(Trim(Mid(degree, 1, InStr(UCase(degree), "''") - 1)))
Else
dsecond = 0
End If
lat = Math.Abs(ddegree) + (Math.Abs(dminute) / 60) + (Math.Abs(dsecond) / 3600)
If ismin Then
lat = lat * -1
End If
Return lat
End Function
Public Sub SetValue(ByVal iid As Integer, ByVal sname As String, ByVal radius As Double, ByVal ecc As Double)
id = iid
ellipsoidName = sname
EquatorialRadius = radius
eccentricitySquared = ecc
End Sub
Public Sub SetID(ByVal id)
If id = -1 Then
SetValue(-1, "Placeholder", 0, 0)
ElseIf id = 1 Then
SetValue(1, "Airy", 6377563, 0.00667054)
ElseIf id = 2 Then
SetValue(2, "Australian National", 6378160, 0.006694542)
ElseIf id = 3 Then
SetValue(3, "Bessel 1841", 6377397, 0.006674372)
ElseIf id = 4 Then
SetValue(4, "Bessel 1841 (Nambia) ", 6377484, 0.006674372)
ElseIf id = 5 Then
SetValue(5, "Clarke 1866", 6378206, 0.006768658)
ElseIf id = 6 Then
SetValue(6, "Clarke 1880", 6378249, 0.006803511)
ElseIf id = 7 Then
SetValue(7, "Everest", 6377276, 0.006637847)
ElseIf id = 8 Then
SetValue(8, "Fischer 1960 (Mercury) ", 6378166, 0.006693422)
ElseIf id = 9 Then
SetValue(9, "Fischer 1968", 6378150, 0.006693422)
ElseIf id = 10 Then
SetValue(10, "GRS 1967", 6378160, 0.006694605)
ElseIf id = 11 Then
SetValue(11, "GRS 1980", 6378137, 0.00669438)
ElseIf id = 12 Then
SetValue(12, "Helmert 1906", 6378200, 0.006693422)
ElseIf id = 13 Then
SetValue(13, "Hough", 6378270, 0.00672267)
ElseIf id = 14 Then
SetValue(14, "International", 6378388, 0.00672267)
ElseIf id = 15 Then
SetValue(15, "Krassovsky", 6378245, 0.006693422)
ElseIf id = 16 Then
SetValue(16, "Modified Airy", 6377340, 0.00667054)
ElseIf id = 17 Then
SetValue(17, "Modified Everest", 6377304, 0.006637847)
ElseIf id = 18 Then
SetValue(18, "Modified Fischer 1960", 6378155, 0.006693422)
ElseIf id = 19 Then
SetValue(19, "South American 1969", 6378160, 0.006694542)
ElseIf id = 20 Then
SetValue(20, "WGS 60", 6378165, 0.006693422)
ElseIf id = 21 Then
SetValue(21, "WGS 66", 6378145, 0.006694542)
ElseIf id = 22 Then
SetValue(22, "WGS-72", 6378135, 0.006694318)
ElseIf id = 23 Then
SetValue(23, "WGS-84", 6378137, 0.00669438)
End If
End Sub
Function UTMLetterDesignator(ByVal Lat As Double) As Char
Dim LetterDesignator As Char
If ((84 >= Lat) And (Lat >= 72)) Then
LetterDesignator = "X"
ElseIf ((72 > Lat) And (Lat >= 64)) Then
LetterDesignator = "W"
ElseIf ((64 > Lat) And (Lat >= 56)) Then
LetterDesignator = "V"
ElseIf ((56 > Lat) And (Lat >= 48)) Then
LetterDesignator = "U"
ElseIf ((48 > Lat) And (Lat >= 40)) Then
LetterDesignator = "T"
ElseIf ((40 > Lat) And (Lat >= 32)) Then
LetterDesignator = "S"
ElseIf ((32 > Lat) And (Lat >= 24)) Then
LetterDesignator = "R"
ElseIf ((24 > Lat) And (Lat >= 16)) Then
LetterDesignator = "Q"
ElseIf ((16 > Lat) And (Lat >= 8)) Then
LetterDesignator = "P"
ElseIf ((8 > Lat) And (Lat >= 0)) Then
LetterDesignator = "N"
ElseIf ((0 > Lat) And (Lat >= -8)) Then
LetterDesignator = "M"
ElseIf ((-8 > Lat) And (Lat >= -16)) Then
LetterDesignator = "L"
ElseIf ((-16 > Lat) And (Lat >= -24)) Then
LetterDesignator = "K"
ElseIf ((-24 > Lat) And (Lat >= -32)) Then
LetterDesignator = "J"
ElseIf ((-32 > Lat) And (Lat >= -40)) Then
LetterDesignator = "H"
ElseIf ((-40 > Lat) And (Lat >= -48)) Then
LetterDesignator = "G"
ElseIf ((-48 > Lat) And (Lat >= -56)) Then
LetterDesignator = "F"
ElseIf ((-56 > Lat) And (Lat >= -64)) Then
LetterDesignator = "E"
ElseIf ((-64 > Lat) And (Lat >= -72)) Then
LetterDesignator = "D"
ElseIf ((-72 > Lat) And (Lat >= -80)) Then
LetterDesignator = "C"
Else
LetterDesignator = "Z"
End If
Return LetterDesignator
End Function
Sub LLtoUTM(ByVal Lat As Double, ByVal Longi As Double, ByRef UTMeast As Double, ByRef UTMnorth As Double, ByVal UTMzone As String)
Dim eccSquared = eccentricitySquared
Dim k0 As Double = 0.9996
Dim longorigin As Double = 0
Dim eccPrimeSquared As Double = 0
Dim n, t, c, a, m
Dim longtemp As Double = (Longi + 180) - Int((Longi + 180) / 360) * 360 - 180
Dim LatRad As Double = Lat * deg2rad
Dim LongRad As Double = longtemp * deg2rad
Dim LongOriginRad As Double = 0
Dim ZoneNumber As Integer
ZoneNumber = Int((longtemp + 180) / 6) + 1
If (Lat >= 56.0 And Lat < 64.0 And longtemp >= 3.0 And longtemp < 12.0) Then ZoneNumber = 32
If (Lat >= 72.0 And Lat < 84.0) Then
If (longtemp >= 0.0 And longtemp < 9.0) Then
ZoneNumber = 31
ElseIf (longtemp >= 9.0 And longtemp < 21.0) Then
ZoneNumber = 33
ElseIf (longtemp >= 21.0 And longtemp < 33.0) Then
ZoneNumber = 35
ElseIf (longtemp >= 33.0 And longtemp < 42.0) Then
ZoneNumber = 37
End If
End If
longorigin = (ZoneNumber - 1) * 6 - 180 + 3 '+3 puts origin in middle of zone
LongOriginRad = longorigin * deg2rad
'sprintf(UTMZone, "%d%c", ZoneNumber, UTMLetterDesignator(Lat));
UTMzone = Str(ZoneNumber) & UTMLetterDesignator(Lat)
eccPrimeSquared = (eccentricitySquared) / (1 - eccentricitySquared)
n = EquatorialRadius / Math.Sqrt(1 - eccentricitySquared * Math.Sin(LatRad) * Math.Sin(LatRad))
t = Math.Tan(LatRad) * Math.Tan(LatRad)
c = eccPrimeSquared * Math.Cos(LatRad) * Math.Cos(LatRad)
a = Math.Cos(LatRad) * (LongRad - LongOriginRad)
m = EquatorialRadius * ((1 - eccentricitySquared / 4 - 3 * eccentricitySquared * eccentricitySquared / 64 - 5 * eccentricitySquared * eccentricitySquared * eccentricitySquared / 256) * LatRad _
- (3 * eccentricitySquared / 8 + 3 * eccentricitySquared * eccentricitySquared / 32 + 45 * eccentricitySquared * eccentricitySquared * eccentricitySquared / 1024) * Math.Sin(2 * LatRad) _
+ (15 * eccentricitySquared * eccentricitySquared / 256 + 45 * eccentricitySquared * eccentricitySquared * eccentricitySquared / 1024) * Math.Sin(4 * LatRad) _
- (35 * eccentricitySquared * eccentricitySquared * eccentricitySquared / 3072) * Math.Sin(6 * LatRad))
UTMeast = (k0 * n * (a + (1 - t + c) * a * a * a / 6 _
+ (5 - 18 * t + t * t + 72 * c - 58 * eccPrimeSquared) * a * a * a * a * a / 120) _
+ 500000.0)
UTMnorth = (k0 * (m + n * Math.Tan(LatRad) * (a * a / 2 + (5 - t + 9 * c + 4 * c * c) * a * a * a * a / 24 _
+ (61 - 58 * t + t * t + 600 * c - 330 * eccPrimeSquared) * a * a * a * a * a * a / 720)))
If (Lat < 0) Then
UTMnorth += 10000000.0 '10000000 meter offset for southern hemisphere
End If
End Sub
Sub UTMtoLL(ByVal ReferenceEllipsoid As Integer, ByVal UTMNorth As Double, ByVal UTMEast As Double, ByVal UTMZone As Char, _
ByRef Lat As Double, ByRef Longi As Double)
'converts UTM coords to lat/long. Equations from USGS Bulletin 1532
'East Longitudes are positive, West longitudes are negative.
'North latitudes are positive, South latitudes are negative
'Lat and Long are in decimal degrees.
'Written by Chuck Gantz- chuck.gantz@globalstar.com
Dim k0 As Double = 0.9996
Dim a As Double = EquatorialRadius
Dim eccSquared As Double = eccentricitySquared
Dim eccPrimeSquared As Decimal
Dim e1 As Double = (1 - Math.Sqrt(1 - eccSquared)) / (1 + Math.Sqrt(1 - eccSquared))
Dim N1, T1, C1, R1, D, M
Dim LongOrigin As Double = 0
Dim mu, phi1, phi1Rad
Dim x As Double
Dim y As Double
Dim ZoneNumber As Integer
Dim ZoneLetter As String
Dim NorthernHemisphere As Integer '1 for northern hemispher, 0 for southern
x = UTMEast - 500000.0 'remove 500,000 meter offset for longitude
y = UTMNorth
ZoneNumber = Mid(UTMZone, 1, Len(UTMZone) - 1)
ZoneLetter = Mid(UTMZone, Len(UTMZone) - 1)
' strtoul(UTMZone, &ZoneLetter, 10)
If ((Asc(ZoneLetter) - Asc("N")) >= 0) Then
NorthernHemisphere = 1 'point is in northern hemisphere
Else
NorthernHemisphere = 0 'point is in southern hemisphere
y -= 10000000.0 'remove 10,000,000 meter offset used for southern hemisphere
End If
LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3 '+3 puts origin in middle of zone
eccPrimeSquared = (eccSquared) / (1 - eccSquared)
M = y / k0
mu = M / (a * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256))
phi1Rad = mu + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.Sin(2 * mu) _
+ (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.Sin(4 * mu) _
+ (151 * e1 * e1 * e1 / 96) * Math.Sin(6 * mu)
phi1 = phi1Rad * rad2deg
N1 = a / Math.Sqrt(1 - eccSquared * Math.Sin(phi1Rad) * Math.Sin(phi1Rad))
T1 = Math.Tan(phi1Rad) * Math.Tan(phi1Rad)
C1 = eccPrimeSquared * Math.Cos(phi1Rad) * Math.Cos(phi1Rad)
R1 = a * (1 - eccSquared) / Math.Pow(1 - eccSquared * Math.Sin(phi1Rad) * Math.Sin(phi1Rad), 1.5)
D = x / (N1 * k0)
Lat = phi1Rad - (N1 * Math.Tan(phi1Rad) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * eccPrimeSquared) * D * D * D * D / 24 _
+ (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * eccPrimeSquared - 3 * C1 * C1) * D * D * D * D * D * D / 720)
Lat = Lat * rad2deg
Longi = (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * eccPrimeSquared + 24 * T1 * T1) _
* D * D * D * D * D / 120) / Math.Cos(phi1Rad)
Longi = LongOrigin + Longi * rad2deg
End Sub
End Class
This is VB code and run using .NET Fremework 2.0 and tested in google earth and the result is same.