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

yfleury · Jul 15, 2022

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

roumei · Jul 19, 2022

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...

TILogistic · Jul 18, 2022

An option can be:

Create a routine that does the following

Read your geojson file with JsonParse.

And for each location call the api (OkHttpUtils2):

https://epsg.io/trans?x=399871.321219788631424&y=5418344.341168709099293&z=0&s_srs=2950&t_srs=4326

Parse the resulting Json

JSON:

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

Sample MAPS:

See API:

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

yfleury · Jul 19, 2022

It look great. Tomorow morning I will work on it.

roumei · Jul 19, 2022

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

yfleury · Jul 20, 2022

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

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.

roumei · Jul 20, 2022

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.

yfleury · Jul 20, 2022

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]]

roumei · Jul 21, 2022

That looks correct. It's EPSG:2950. The calculation should be as precise as possible.

yfleury · Jul 21, 2022

It's seem the wgs84 have 2m accuracy.

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

yfleury

Active Member

roumei

TILogistic

Expert

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

yfleury

Active Member

roumei

Active Member

yfleury

Active Member

roumei

Active Member

yfleury

Active Member

roumei

Active Member

yfleury

Active Member

Similar Threads

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

Active Member

Expert

Active Member

Active Member

Active Member

Active Member

Active Member

Active Member

Active Member

Similar Threads

Privacy & Transparency

Privacy & Transparency