Android Question how to calculate a NAD83(CSRS) to WGS84

yfleury

yfleury

Active Member
Licensed User
Longtime User
Hi!
I have a database with GPS coordinates in NAD83(CSRS) format.
It's for a google map app
I want to convert it to WGS84 (latitude, longitude)
How to do it, some link
 
Solution
roumei
Omar, thank you for the API explanation!

In case you want to calculate it directly, you can use the following snippet. The result matches the API result. Be careful with the MTM zone, check the prj file for the description of the projection.
B4X: 
Sub Class_Globals
    Private Root As B4XView
    Private xui As XUI
    
    Private const MathPIdiv180 As Double = cPI / 180
    Private Const WGS84_a As Double = 6378137
    Private Const WGS84_b As Double = 6356752.3142
    
    Type eGPSLocation(Latitude As Double, Longitude As Double)
    Type eProjCoordinate(Easting As Double, Northing As Double)
    
End Sub

Public Sub Initialize
'    B4XPages.GetManager.LogEvents = True
End Sub

'This event will be called once, before the page becomes...
Sort by date Sort by votes
TILogistic

TILogistic

Expert
Licensed User
Longtime User
An option can be:

Create a routine that does the following

Read your geojson file with JsonParse.


1658184374542.png


And for each location call the api (OkHttpUtils2):


Parse the resulting Json

JSON: 
{"x":"-72.20326051","y":"48.89612131","z":"0.00000000"}

1658184546775.png


Sample MAPS:
1658184767862.png


See API:
github.com

GitHub - maptiler/epsg.io: EPSG.io: Coordinate Systems Worldwide

EPSG.io: Coordinate Systems Worldwide. Contribute to maptiler/epsg.io development by creating an account on GitHub.
github.com github.com
1658185617990.png
 
Last edited:
Upvote 0
roumei

roumei

Active Member
Licensed User
Omar, thank you for the API explanation!

In case you want to calculate it directly, you can use the following snippet. The result matches the API result. Be careful with the MTM zone, check the prj file for the description of the projection.
B4X: 
Sub Class_Globals
    Private Root As B4XView
    Private xui As XUI
    
    Private const MathPIdiv180 As Double = cPI / 180
    Private Const WGS84_a As Double = 6378137
    Private Const WGS84_b As Double = 6356752.3142
    
    Type eGPSLocation(Latitude As Double, Longitude As Double)
    Type eProjCoordinate(Easting As Double, Northing As Double)
    
End Sub

Public Sub Initialize
'    B4XPages.GetManager.LogEvents = True
End Sub

'This event will be called once, before the page becomes visible.
Private Sub B4XPage_Created (Root1 As B4XView)
    Root = Root1
    Root.LoadLayout("MainPage")
End Sub

'You can see the list of page related events in the B4XPagesManager object. The event name is B4XPage.

Private Sub Button1_Click
    
    ' EPSG.IO API, thanks to Omar Parra A.
    ' https://epsg.io/trans?x=399871.321219788631424&y=5418344.341168709099293&z=0&s_srs=2950&t_srs=4326
    ' Result: {"x":"-72.20326051","y":"48.89612131","z":"0.00000000"}
    
    ' Calculation
    Dim iMTMZone As Int = 8
    Dim r As eGPSLocation = GetLatLonFromMTMNorthCoordinate(399871.321219788631424, 5418344.341168709099293, iMTMZone)
    
    Log(NumberFormat(r.Latitude, 1, 8) & " " & NumberFormat(r.Longitude, 1, 8))
    ' 48.89612131 -72.20326051
    
End Sub

Public Sub GetMTMNorthCoordinateFromLatLon(Latitude As Double, Longitude As Double, iMTMZone As Int) As eProjCoordinate
    If iMTMZone < 4 Or iMTMZone > 17 Then Return Null
    Return GPS_GetTransverseMercator_FromLatLon(Latitude, Longitude, 304800, 0, GetCentralMeridianForMTMZone(iMTMZone), 0.9999, 0, WGS84_a, WGS84_b)
End Sub

Public Sub GetLatLonFromMTMNorthCoordinate(Easting As Double, Northing As Double, iMTMZone As Int) As eGPSLocation
    If iMTMZone < 4 Or iMTMZone > 17 Then Return Null
    Return GPS_GetLatLon_FromTransverseMercator(Easting, Northing, 304800, 0, GetCentralMeridianForMTMZone(iMTMZone), 0.9999, 0, WGS84_a, WGS84_b)
End Sub

Private Sub GetCentralMeridianForMTMZone(iZone As Int) As Double
    If iZone >= 4 And iZone <= 11 Then
        Return -61.5 - (iZone - 4) * 3
    End If
    If iZone >= 12 And iZone <= 17 Then
        Return -81.0 - (iZone - 12) * 3
    End If
    Return 0
End Sub

Private Sub GPS_GetTransverseMercator_FromLatLon(Latitude As Double, _
                                                         Longitude As Double, _
                                                         FalseEasting As Double, _
                                                         FalseNorthing As Double, _
                                                         CentralMeridian As Double, _
                                                         CentralMeridianScaleFactor As Double, _
                                                         LatitudeOfOrigin As Double, _
                                                         Ellipsoid_A As Double, _
                                                         Ellipsoid_B As Double) As eProjCoordinate
    
    Dim Northing, Easting, eccPrimeSquared, nu, T, C, A, M, M0 As Double
    Dim eccSquared As Double = Ellipsoid_B / Ellipsoid_A
    eccSquared = 1 - eccSquared * eccSquared

    If Longitude < -180 Then Longitude = Longitude + 360
    If Longitude > 180 Then Longitude = Longitude - 360
    Dim LatRad As Double = Latitude * MathPIdiv180
    Dim LongRad As Double = Longitude * MathPIdiv180
    Dim LongOriginRad As Double = CentralMeridian * MathPIdiv180
    Dim LatOriginRad As Double = LatitudeOfOrigin * MathPIdiv180
    eccPrimeSquared = eccSquared / (1 - eccSquared)
    nu = Ellipsoid_A / Sqrt(1 - eccSquared * Sin(LatRad) * Sin(LatRad))
    T = Tan(LatRad) * Tan(LatRad)
    C = eccPrimeSquared * Cos(LatRad) * Cos(LatRad)
    A = Cos(LatRad) * (LongRad - LongOriginRad)
    M = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatRad))
    M0 = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatOriginRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatOriginRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatOriginRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatOriginRad))
    Easting = (CentralMeridianScaleFactor * nu * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120))
    Northing = (CentralMeridianScaleFactor * (M - M0 + nu * 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)))

    Dim Coord As eProjCoordinate
    Coord.Easting = Easting + FalseEasting
    Coord.Northing = Northing + FalseNorthing
    Return Coord

End Sub

Private Sub GPS_GetLatLon_FromTransverseMercator(Easting As Double, _
                                                          Northing As Double, _
                                                          FalseEasting As Double, _
                                                          FalseNorthing As Double, _
                                                          CentralMeridian As Double, _
                                                          CentralMeridianScaleFactor As Double, _
                                                          LatitudeOfOrigin As Double, _
                                                          Ellipsoid_A As Double, _
                                                          Ellipsoid_B As Double) As eGPSLocation

    Dim Lat, Lon, N1, T1, C1, R1, D, M1, mu1, phi1Rad As Double
    Dim LongOriginRad As Double = CentralMeridian * MathPIdiv180
    Dim LatOriginRad As Double = LatitudeOfOrigin * MathPIdiv180
    Dim eccSquared As Double = Ellipsoid_B / Ellipsoid_A
    eccSquared = 1 - eccSquared * eccSquared
    Dim e1 As Double = (1 - Sqrt(1 - eccSquared)) / (1 + Sqrt(1 - eccSquared))
    Dim eccPrimeSquared As Double = (eccSquared) / (1 - eccSquared)
    Dim M0 As Double = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatOriginRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatOriginRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatOriginRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatOriginRad))
    M1 = M0 + (Northing - FalseNorthing) / CentralMeridianScaleFactor
    mu1 = M1 / (Ellipsoid_A * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256))
    phi1Rad = mu1 + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Sin(2 * mu1) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Sin(4 * mu1) + (151 * e1 * e1 * e1 / 96) * Sin(6 * mu1)
    N1 = Ellipsoid_A / Sqrt(1 - eccSquared * Sin(phi1Rad) * Sin(phi1Rad))
    T1 = Tan(phi1Rad) * Tan(phi1Rad)
    C1 = eccPrimeSquared * Cos(phi1Rad) * Cos(phi1Rad)
    R1 = Ellipsoid_A * (1 - eccSquared) / Power(1 - eccSquared * Sin(phi1Rad) * Sin(phi1Rad), 1.5)
    D = (Easting - FalseEasting) / (N1 * CentralMeridianScaleFactor)
    Lat = phi1Rad - (N1 * 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)
    Lon = (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) / Cos(phi1Rad)

    Dim e As eGPSLocation
    e.Latitude = Lat * (180.0 / cPI)
    e.Longitude = (Lon + LongOriginRad) * (180.0 / cPI)
    Return e

End Sub
 
Upvote 1
Solution
yfleury

yfleury

Active Member
Licensed User
Longtime User
roumei said:
Omar, thank you for the API explanation!

In case you want to calculate it directly, you can use the following snippet. The result matches the API result. Be careful with the MTM zone, check the prj file for the description of the projection.
B4X: 
Sub Class_Globals
    Private Root As B4XView
    Private xui As XUI
   
    Private const MathPIdiv180 As Double = cPI / 180
    Private Const WGS84_a As Double = 6378137
    Private Const WGS84_b As Double = 6356752.3142
   
    Type eGPSLocation(Latitude As Double, Longitude As Double)
    Type eProjCoordinate(Easting As Double, Northing As Double)
   
End Sub

Public Sub Initialize
'    B4XPages.GetManager.LogEvents = True
End Sub

'This event will be called once, before the page becomes visible.
Private Sub B4XPage_Created (Root1 As B4XView)
    Root = Root1
    Root.LoadLayout("MainPage")
End Sub

'You can see the list of page related events in the B4XPagesManager object. The event name is B4XPage.

Private Sub Button1_Click
   
    ' EPSG.IO API, thanks to Omar Parra A.
    ' https://epsg.io/trans?x=399871.321219788631424&y=5418344.341168709099293&z=0&s_srs=2950&t_srs=4326
    ' Result: {"x":"-72.20326051","y":"48.89612131","z":"0.00000000"}
   
    ' Calculation
    Dim iMTMZone As Int = 8
    Dim r As eGPSLocation = GetLatLonFromMTMNorthCoordinate(399871.321219788631424, 5418344.341168709099293, iMTMZone)
   
    Log(NumberFormat(r.Latitude, 1, 8) & " " & NumberFormat(r.Longitude, 1, 8))
    ' 48.89612131 -72.20326051
   
End Sub

Public Sub GetMTMNorthCoordinateFromLatLon(Latitude As Double, Longitude As Double, iMTMZone As Int) As eProjCoordinate
    If iMTMZone < 4 Or iMTMZone > 17 Then Return Null
    Return GPS_GetTransverseMercator_FromLatLon(Latitude, Longitude, 304800, 0, GetCentralMeridianForMTMZone(iMTMZone), 0.9999, 0, WGS84_a, WGS84_b)
End Sub

Public Sub GetLatLonFromMTMNorthCoordinate(Easting As Double, Northing As Double, iMTMZone As Int) As eGPSLocation
    If iMTMZone < 4 Or iMTMZone > 17 Then Return Null
    Return GPS_GetLatLon_FromTransverseMercator(Easting, Northing, 304800, 0, GetCentralMeridianForMTMZone(iMTMZone), 0.9999, 0, WGS84_a, WGS84_b)
End Sub

Private Sub GetCentralMeridianForMTMZone(iZone As Int) As Double
    If iZone >= 4 And iZone <= 11 Then
        Return -61.5 - (iZone - 4) * 3
    End If
    If iZone >= 12 And iZone <= 17 Then
        Return -81.0 - (iZone - 12) * 3
    End If
    Return 0
End Sub

Private Sub GPS_GetTransverseMercator_FromLatLon(Latitude As Double, _
                                                         Longitude As Double, _
                                                         FalseEasting As Double, _
                                                         FalseNorthing As Double, _
                                                         CentralMeridian As Double, _
                                                         CentralMeridianScaleFactor As Double, _
                                                         LatitudeOfOrigin As Double, _
                                                         Ellipsoid_A As Double, _
                                                         Ellipsoid_B As Double) As eProjCoordinate
   
    Dim Northing, Easting, eccPrimeSquared, nu, T, C, A, M, M0 As Double
    Dim eccSquared As Double = Ellipsoid_B / Ellipsoid_A
    eccSquared = 1 - eccSquared * eccSquared

    If Longitude < -180 Then Longitude = Longitude + 360
    If Longitude > 180 Then Longitude = Longitude - 360
    Dim LatRad As Double = Latitude * MathPIdiv180
    Dim LongRad As Double = Longitude * MathPIdiv180
    Dim LongOriginRad As Double = CentralMeridian * MathPIdiv180
    Dim LatOriginRad As Double = LatitudeOfOrigin * MathPIdiv180
    eccPrimeSquared = eccSquared / (1 - eccSquared)
    nu = Ellipsoid_A / Sqrt(1 - eccSquared * Sin(LatRad) * Sin(LatRad))
    T = Tan(LatRad) * Tan(LatRad)
    C = eccPrimeSquared * Cos(LatRad) * Cos(LatRad)
    A = Cos(LatRad) * (LongRad - LongOriginRad)
    M = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatRad))
    M0 = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatOriginRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatOriginRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatOriginRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatOriginRad))
    Easting = (CentralMeridianScaleFactor * nu * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120))
    Northing = (CentralMeridianScaleFactor * (M - M0 + nu * 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)))

    Dim Coord As eProjCoordinate
    Coord.Easting = Easting + FalseEasting
    Coord.Northing = Northing + FalseNorthing
    Return Coord

End Sub

Private Sub GPS_GetLatLon_FromTransverseMercator(Easting As Double, _
                                                          Northing As Double, _
                                                          FalseEasting As Double, _
                                                          FalseNorthing As Double, _
                                                          CentralMeridian As Double, _
                                                          CentralMeridianScaleFactor As Double, _
                                                          LatitudeOfOrigin As Double, _
                                                          Ellipsoid_A As Double, _
                                                          Ellipsoid_B As Double) As eGPSLocation

    Dim Lat, Lon, N1, T1, C1, R1, D, M1, mu1, phi1Rad As Double
    Dim LongOriginRad As Double = CentralMeridian * MathPIdiv180
    Dim LatOriginRad As Double = LatitudeOfOrigin * MathPIdiv180
    Dim eccSquared As Double = Ellipsoid_B / Ellipsoid_A
    eccSquared = 1 - eccSquared * eccSquared
    Dim e1 As Double = (1 - Sqrt(1 - eccSquared)) / (1 + Sqrt(1 - eccSquared))
    Dim eccPrimeSquared As Double = (eccSquared) / (1 - eccSquared)
    Dim M0 As Double = Ellipsoid_A * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatOriginRad - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(2 * LatOriginRad) + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * Sin(4 * LatOriginRad) - (35 * eccSquared * eccSquared * eccSquared / 3072) * Sin(6 * LatOriginRad))
    M1 = M0 + (Northing - FalseNorthing) / CentralMeridianScaleFactor
    mu1 = M1 / (Ellipsoid_A * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256))
    phi1Rad = mu1 + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Sin(2 * mu1) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Sin(4 * mu1) + (151 * e1 * e1 * e1 / 96) * Sin(6 * mu1)
    N1 = Ellipsoid_A / Sqrt(1 - eccSquared * Sin(phi1Rad) * Sin(phi1Rad))
    T1 = Tan(phi1Rad) * Tan(phi1Rad)
    C1 = eccPrimeSquared * Cos(phi1Rad) * Cos(phi1Rad)
    R1 = Ellipsoid_A * (1 - eccSquared) / Power(1 - eccSquared * Sin(phi1Rad) * Sin(phi1Rad), 1.5)
    D = (Easting - FalseEasting) / (N1 * CentralMeridianScaleFactor)
    Lat = phi1Rad - (N1 * 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)
    Lon = (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) / Cos(phi1Rad)

    Dim e As eGPSLocation
    e.Latitude = Lat * (180.0 / cPI)
    e.Longitude = (Lon + LongOriginRad) * (180.0 / cPI)
    Return e

End Sub
Click to expand...
Your code it's great and work like I want.
But some time the accuracy is not good like 16 meters away. Also, I need to check the x,y of nad83 to a web site if it is same position. Feedback will come later.
 
Upvote 0
roumei

roumei

Active Member
Licensed User
How do you determine that the accuracy is not good?
Can you upload the prj file? It should contain the info about the projection. Maybe it's something different than the assumed EPSG:2950.
 
Upvote 0
yfleury

yfleury

Active Member
Licensed User
Longtime User
In the .prf file, I have this
B4X: 
PROJCS["NAD_1983_CSRS_MTM_8",GEOGCS["GCS_North_American_1983_CSRS",DATUM["D_North_American_1983_CSRS",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",304800.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-73.5],PARAMETER["Scale_Factor",0.9999],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
 
Upvote 0
You must log in or register to reply here.

Similar Threads

S
  • Question
Android Question Possible bug with warwound Spatialite Library
Replies
3
Views
834
synasir
S
I
  • Article
B4A Library FusedLocationProviderGMS (Latest)
2
Replies
21
Views
6K
Ivica Golubovic
I
Terradrones
  • Question
Android Question Enum
Replies
8
Views
2K
ddk1
D
booch34
  • Article
Android Code Snippet WGS 84 reader
Replies
0
Views
1K
booch34
booch34
Erel
  • Article
B4A Library [B4X] [class] Geodesic - Convert Lat Lon / UTM coordinates
Replies
10
Views
13K
jamesnz
J
Top