Skip to content

Instantly share code, notes, and snippets.

@marsyang1

marsyang1/TWD97_change_to_WGS.vb

Last active July 6, 2022 02:39
Show Gist options
  • Save marsyang1/f5f87004b98d438159bf7c4acad7aef8 to your computer and use it in GitHub Desktop.
Save marsyang1/f5f87004b98d438159bf7c4acad7aef8 to your computer and use it in GitHub Desktop.
Download ZIP
TWD97_change_to_WGS VB測試
Raw
TWD97_change_to_WGS.vb
' Online VB Compiler.
' Code, Compile, Run and Debug VB program online.
' Write your code in this editor and press "Run" button to execute it.
'
' 參考來源
' https://github.com/Chao-wei-chu/TWD97_change_to_WGS
' 因為朋友需求 , 試試調整成vb結構
'
Module VBModule
Sub Main()
' // Target Validate converter page
' // https://wdms.epa.gov.tw/idms/coordtrans.htm
' //
' // 宣告調整參數
' // https://github.com/Chao-wei-chu/TWD97_change_to_WGS
Dim dx As Double = 250000.0
Dim dy As Double = 0.0
Dim lon0 As Double = 121 * Math.PI / 180
Dim k0 As Double = 0.9999
Dim a As Double = 6378137.0
Dim b As Double = 6356752.314245
Dim e As Double = Math.pow((1- Math.pow(b,2)/ Math.pow(a,2)), 0.5)
Console.WriteLine("dx:"+dx.ToString)
Console.WriteLine("dy:"+dy.ToString)
Console.WriteLine("lon0:"+lon0.ToString)
Console.WriteLine("k0:"+k0.ToString)
Console.WriteLine("a:"+a.ToString)
Console.WriteLine("b:"+b.ToString)
Console.WriteLine("e:"+e.ToString)
' // 宣告來源
Dim sourceX As Double = 250000.000
Dim sourceY As Double = 2544283.154
Dim x As Double = sourceX - dx
Dim y As Double = sourceY - dy
' // Calculate the Meridional Arc
Dim M As Double = y/k0
Dim mu As Double = M/(a*(1.0 - Math.pow(e, 2)/4.0 - 3*Math.pow(e, 4)/64.0 - 5*Math.pow(e, 6)/256.0))
Dim e1 As Double = (1.0 - Math.pow((1.0 - Math.pow(e, 2)), 0.5)) / (1.0 + Math.pow((1.0 - Math.pow(e, 2)), 0.5))
Dim J1 As Double = (3*e1/2 - 27*Math.pow(e1, 3)/32.0)
Dim J2 As Double = (21*Math.pow(e1, 2)/16 - 55*Math.pow(e1, 4)/32.0)
Dim J3 As Double = (151*Math.pow(e1, 3)/96.0)
Dim J4 As Double = (1097*Math.pow(e1, 4)/512.0)
Dim fp As Double = mu + J1*Math.sin(2*mu) + J2*Math.sin(4*mu) + J3*Math.sin(6*mu) + J4*Math.sin(8*mu)
' // Calculate Latitude and Longitude
Dim e2 As Double = Math.pow((e*a/b), 2)
Dim C1 As Double = Math.pow(e2*Math.cos(fp), 2)
Dim T1 As Double = Math.pow(Math.tan(fp), 2)
Dim R1 As Double = a*(1-Math.pow(e, 2))/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), (3.0/2.0))
Dim N1 As Double = a/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), 0.5)
Dim D As Double = x/(N1*k0)
' // lat
Dim Q1 As Double = N1*Math.tan(fp)/R1
Dim Q2 As Double = (Math.pow(D, 2)/2.0)
Dim Q3 As Double = (5 + 3*T1 + 10*C1 - 4*Math.pow(C1, 2) - 9*e2)*Math.pow(D, 4)/24.0
Dim Q4 As Double = (61 + 90*T1 + 298*C1 + 45*Math.pow(T1, 2) - 3*Math.pow(C1, 2) - 252*e2)*Math.pow(D, 6)/720.0
Dim lat As Double = fp - Q1*(Q2 - Q3 + Q4)
' // long
Dim Q5 As Double = D
Dim Q6 As Double = (1 + 2*T1 + C1)*Math.pow(D, 3)/6
Dim Q7 As Double = (5 - 2*C1 + 28*T1 - 3*Math.pow(C1, 2) + 8*e2 + 24*Math.pow(T1, 2))*Math.pow(D, 5)/120.0
Dim lon As Double = lon0 + (Q5 - Q6 + Q7)/Math.cos(fp)
Console.WriteLine("lat:"+lat.ToString)
Console.WriteLine("lon:"+lon.ToString)
Dim latitude As Double = radians_to_degrees(lat)
Dim longitude As Double = radians_to_degrees(lon)
Console.WriteLine("Final Result:=====================================")
Console.WriteLine("Source X:"+sourceX.ToString)
Console.WriteLine("Source Y:"+sourceY.ToString)
Console.WriteLine("-------------------------------")
Console.WriteLine("latitude:"+latitude.ToString)
Console.WriteLine("longitude:"+longitude.ToString)
End Sub
' 用到的參考
' https://theprogrammingexpert.com/vba-radians-to-degrees/
Function radians_to_degrees(radians as Double) as Double
radians_to_degrees = radians * (180 / Math.Pi)
End Function
End Module
Raw
TWD97_change_to_WGS_v2.java
package com.example.testgis;
import org.junit.jupiter.api.Test;
import org.springframework.boot.test.context.SpringBootTest;
import java.util.HashMap;
import java.util.Map;
/**
* 根據原本的code調整成單一執行Class
*/
@SpringBootTest
class TestGisApplicationTests {
@Test
void contextLoads() {
System.out.println("======================== Coding start ============================ ");
// Target Validate converter page
// https://wdms.epa.gov.tw/idms/coordtrans.htm
//
// 宣告調整參數
// https://github.com/Chao-wei-chu/TWD97_change_to_WGS
final Map<String,Double> tmParameter = new HashMap();
tmParameter.put("dx",250000.0);
tmParameter.put("dy",0.0);
tmParameter.put("lon0",121 * Math.PI / 180);
tmParameter.put("k0",0.9999);
tmParameter.put("a",6378137.0);
tmParameter.put("b",6356752.314245);
// 宣告來源
final double sourceX = 250000.000;
final double sourceY = 2544283.154;
final double dx = tmParameter.get("dx");
final double dy = tmParameter.get("dy");
final double lon0 = tmParameter.get("lon0");
final double k0 = tmParameter.get("k0");
final double a = tmParameter.get("a");
final double b = tmParameter.get("b");
final double e = Math.pow((1- Math.pow(b,2)/ Math.pow(a,2)), 0.5);
double x = sourceX - dx;
double y = sourceY - dy;
// Calculate the Meridional Arc
double M = y/k0;
System.out.println("x:"+x);
System.out.println("y:"+y);
System.out.println("M:"+M);
// Calculate Footprint Latitude
double mu = M/(a*(1.0 - Math.pow(e, 2)/4.0 - 3*Math.pow(e, 4)/64.0 - 5*Math.pow(e, 6)/256.0));
double e1 = (1.0 - Math.pow((1.0 - Math.pow(e, 2)), 0.5)) / (1.0 + Math.pow((1.0 - Math.pow(e, 2)), 0.5));
double J1 = (3*e1/2 - 27*Math.pow(e1, 3)/32.0);
double J2 = (21*Math.pow(e1, 2)/16 - 55*Math.pow(e1, 4)/32.0);
double J3 = (151*Math.pow(e1, 3)/96.0);
double J4 = (1097*Math.pow(e1, 4)/512.0);
double fp = mu + J1*Math.sin(2*mu) + J2*Math.sin(4*mu) + J3*Math.sin(6*mu) + J4*Math.sin(8*mu);
// Calculate Latitude and Longitude
double e2 = Math.pow((e*a/b), 2);
double C1 = Math.pow(e2*Math.cos(fp), 2);
double T1 = Math.pow(Math.tan(fp), 2);
double R1 = a*(1-Math.pow(e, 2))/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), (3.0/2.0));
double N1 = a/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), 0.5);
double D = x/(N1*k0);
// lat
double Q1 = N1*Math.tan(fp)/R1;
double Q2 = (Math.pow(D, 2)/2.0);
double Q3 = (5 + 3*T1 + 10*C1 - 4*Math.pow(C1, 2) - 9*e2)*Math.pow(D, 4)/24.0;
double Q4 = (61 + 90*T1 + 298*C1 + 45*Math.pow(T1, 2) - 3*Math.pow(C1, 2) - 252*e2)*Math.pow(D, 6)/720.0;
double lat = fp - Q1*(Q2 - Q3 + Q4);
System.out.println("lat:"+lat);
// long
double Q5 = D;
double Q6 = (1 + 2*T1 + C1)*Math.pow(D, 3)/6;
double Q7 = (5 - 2*C1 + 28*T1 - 3*Math.pow(C1, 2) + 8*e2 + 24*Math.pow(T1, 2))*Math.pow(D, 5)/120.0;
double lon = lon0 + (Q5 - Q6 + Q7)/Math.cos(fp);
System.out.println("lon:"+lon);
final String latitude = ""+Math.toDegrees(lat);
final String longitude = ""+Math.toDegrees(lon);
System.out.println("目標值經度:"+longitude+",緯度:"+latitude);
System.out.println("======================== Coding End ============================ ");
}
}
Raw
TWD97_change_to_WGS_v2.vb
' Online VB Compiler.
' Code, Compile, Run and Debug VB program online.
' Write your code in this editor and press "Run" button to execute it.
'
' 參考來源
' https://github.com/Chao-wei-chu/TWD97_change_to_WGS
' 因為朋友需求 , 試試調整成vb結構
'
' 線上測試用編輯器
' https://www.onlinegdb.com/online_vb_compiler
'
' 計算結果參考對照站台]
' https://wdms.epa.gov.tw/idms/coordtrans.htm
Module VBModule
' // 宣告測試用來源
Const sourceX As Double = 250000.000
Const sourceY As Double = 2544283.154
' // Target Validate converter page
' // https://wdms.epa.gov.tw/idms/coordtrans.htm
' //
' // 宣告調整參數
' // https://github.com/Chao-wei-chu/TWD97_change_to_WGS
Const dx As Double = 250000.0
Const dy As Double = 0.0
Const lon0 As Double = 121 * Math.PI / 180
Const k0 As Double = 0.9999
Const a As Double = 6378137.0
Const b As Double = 6356752.314245
' //初始化計算用相關參數
Dim e As Double = Math.pow((1- Math.pow(b,2)/ Math.pow(a,2)), 0.5)
Dim e1 As Double = (1.0 - Math.pow((1.0 - Math.pow(e, 2)), 0.5)) / (1.0 + Math.pow((1.0 - Math.pow(e, 2)), 0.5))
Dim e2 As Double = Math.pow((e*a/b), 2)
Dim J1 As Double = (3*e1/2 - 27*Math.pow(e1, 3)/32.0)
Dim J2 As Double = (21*Math.pow(e1, 2)/16 - 55*Math.pow(e1, 4)/32.0)
Dim J3 As Double = (151*Math.pow(e1, 3)/96.0)
Dim J4 As Double = (1097*Math.pow(e1, 4)/512.0)
Sub Main()
' Debug用
' printPrepareParam()
Dim wgsObj As Double() = convertToWGS(sourceX,sourceY)
Dim latitude As Double = wgsObj(1)
Dim longitude As Double =wgsObj(2)
Console.WriteLine("Debug用 , 計算結果================================")
Console.WriteLine("來源X:"+sourceX.ToString)
Console.WriteLine("來源Y:"+sourceY.ToString)
Console.WriteLine("------------------------")
Console.WriteLine("經度wgsObj.longitude:"+longitude.ToString)
Console.WriteLine("緯度wgsObj.latitude:"+latitude.ToString)
Console.WriteLine("Debug用 , 計算結束================================")
' 這段是為了要給Excel 用
' With Sheets("Main")
' Dim wgsObj As Double()
' For r = 3 To .Cells(.Rows.Count, 1).End(xlUp).Row
' x = .Cells(r, 1)
' y = .Cells(r, 2)
' Debug.Print "x=" & x
' Debug.Print "y=" & y
' wgsObj = convertToWGS(x,y)
' .Cells(r, 4) = wgsObj(1)
' .Cells(r, 5) = wgsObj(2)
' Next
End Sub
Function convertToWGS(twd97x as Double , twd97y as Double) As Double()
Dim x As Double = twd97x - dx
Dim y As Double = twd97y - dy
' // Calculate the Meridional Arc
Dim M As Double = y/k0
Dim mu As Double = M/(a*(1.0 - Math.pow(e, 2)/4.0 - 3*Math.pow(e, 4)/64.0 - 5*Math.pow(e, 6)/256.0))
Dim fp As Double = mu + J1*Math.sin(2*mu) + J2*Math.sin(4*mu) + J3*Math.sin(6*mu) + J4*Math.sin(8*mu)
' // Calculate Latitude and Longitude
Dim C1 As Double = Math.pow(e2*Math.cos(fp), 2)
Dim T1 As Double = Math.pow(Math.tan(fp), 2)
Dim R1 As Double = a*(1-Math.pow(e, 2))/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), (3.0/2.0))
Dim N1 As Double = a/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), 0.5)
Dim D As Double = x/(N1*k0)
' // lat
Dim Q1 As Double = N1*Math.tan(fp)/R1
Dim Q2 As Double = (Math.pow(D, 2)/2.0)
Dim Q3 As Double = (5 + 3*T1 + 10*C1 - 4*Math.pow(C1, 2) - 9*e2)*Math.pow(D, 4)/24.0
Dim Q4 As Double = (61 + 90*T1 + 298*C1 + 45*Math.pow(T1, 2) - 3*Math.pow(C1, 2) - 252*e2)*Math.pow(D, 6)/720.0
Dim lat As Double = fp - Q1*(Q2 - Q3 + Q4)
' // long
Dim Q5 As Double = D
Dim Q6 As Double = (1 + 2*T1 + C1)*Math.pow(D, 3)/6
Dim Q7 As Double = (5 - 2*C1 + 28*T1 - 3*Math.pow(C1, 2) + 8*e2 + 24*Math.pow(T1, 2))*Math.pow(D, 5)/120.0
Dim lon As Double = lon0 + (Q5 - Q6 + Q7)/Math.cos(fp)
Console.WriteLine("lat:"+lat.ToString)
Console.WriteLine("lon:"+lon.ToString)
Dim latitude As Double = radians_to_degrees(lat)
Dim longitude As Double = radians_to_degrees(lon)
Console.WriteLine("緯度latitude:"+latitude.ToString)
Console.WriteLine("經度longitude:"+longitude.ToString)
Dim sets(5) As Double
sets(1) = latitude
sets(2) = longitude
Return(sets)
End Function
' 用到的參考
' https://theprogrammingexpert.com/vba-radians-to-degrees/
Function radians_to_degrees(radians as Double) as Double
radians_to_degrees = radians * (180 / Math.Pi)
End Function
Sub printPrepareParam()
Console.WriteLine("dx:"+dx.ToString)
Console.WriteLine("dy:"+dy.ToString)
Console.WriteLine("lon0:"+lon0.ToString)
Console.WriteLine("k0:"+k0.ToString)
Console.WriteLine("a:"+a.ToString)
Console.WriteLine("b:"+b.ToString)
Console.WriteLine("e:"+e.ToString)
Console.WriteLine("e1:"+e1.ToString)
Console.WriteLine("e2:"+e2.ToString)
Console.WriteLine("J1:"+J1.ToString)
Console.WriteLine("J2:"+J2.ToString)
Console.WriteLine("J3:"+J3.ToString)
Console.WriteLine("J4:"+J4.ToString)
End Sub
End Module
Raw
TWD97_change_to_WGS_v3.vb
' Online VB Compiler.
' Code, Compile, Run and Debug VB program online.
' Write your code in this editor and press "Run" button to execute it.
'
' 參考來源
' https://github.com/Chao-wei-chu/TWD97_change_to_WGS
' 因為朋友需求 , 試試調整成vb結構
'
' 線上測試用編輯器
' https://www.onlinegdb.com/online_vb_compiler
'
' 計算結果參考對照站台]
' https://wdms.epa.gov.tw/idms/coordtrans.htm
Module VBModule
' // 宣告測試用來源
Const sourceX As Double = 250000.000
Const sourceY As Double = 2544283.154
' // Target Validate converter page
' // https://wdms.epa.gov.tw/idms/coordtrans.htm
' //
' // 宣告調整參數
' // https://github.com/Chao-wei-chu/TWD97_change_to_WGS
Const dx As Double = 250000.0
Const dy As Double = 0.0
Const lon0 As Double = 121 * Math.PI / 180
Const k0 As Double = 0.9999
Const a As Double = 6378137.0
Const b As Double = 6356752.314245
' //初始化計算用相關參數
Dim e As Double = 0.0
Dim e1 As Double = 0.0
Dim e2 As Double = 0.0
Dim J1 As Double = 0.0
Dim J2 As Double = 0.0
Dim J3 As Double = 0.0
Dim J4 As Double = 0.0
Sub Main()
' Debug用
prepareParam()
printPrepareParam()
Dim wgsObj As Double() = convertToWGS(sourceX,sourceY)
Dim latitude As Double = 0.0
Dim longitude As Double = 0.0
latitude = wgsObj(1)
longitude = wgsObj(2)
Console.WriteLine("Debug用 , 計算結果================================")
Console.WriteLine("來源X:"+sourceX.ToString)
Console.WriteLine("來源Y:"+sourceY.ToString)
Console.WriteLine("------------------------")
Console.WriteLine("經度wgsObj.longitude:"+longitude.ToString)
Console.WriteLine("緯度wgsObj.latitude:"+latitude.ToString)
Console.WriteLine("Debug用 , 計算結束================================")
' 這段是為了要給Excel 用
' With Sheets("Main")
' Dim wgsObj As Double()
' For r = 3 To .Cells(.Rows.Count, 1).End(xlUp).Row
' x = .Cells(r, 1)
' y = .Cells(r, 2)
' Debug.Print "x=" & x
' Debug.Print "y=" & y
' wgsObj = convertToWGS(x,y)
' .Cells(r, 4) = wgsObj(1)
' .Cells(r, 5) = wgsObj(2)
' Next
End Sub
Function convertToWGS(twd97x as Double , twd97y as Double) As Double()
Dim x As Double = twd97x - dx
Dim y As Double = twd97y - dy
' // Calculate the Meridional Arc
Dim M As Double = y/k0
Dim mu As Double = M/(a*(1.0 - Math.pow(e, 2)/4.0 - 3*Math.pow(e, 4)/64.0 - 5*Math.pow(e, 6)/256.0))
Dim fp As Double = mu + J1*Math.sin(2*mu) + J2*Math.sin(4*mu) + J3*Math.sin(6*mu) + J4*Math.sin(8*mu)
' // Calculate Latitude and Longitude
Dim C1 As Double = Math.pow(e2*Math.cos(fp), 2)
Dim T1 As Double = Math.pow(Math.tan(fp), 2)
Dim R1 As Double = a*(1-Math.pow(e, 2))/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), (3.0/2.0))
Dim N1 As Double = a/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), 0.5)
Dim D As Double = x/(N1*k0)
' // lat
Dim Q1 As Double = N1*Math.tan(fp)/R1
Dim Q2 As Double = (Math.pow(D, 2)/2.0)
Dim Q3 As Double = (5 + 3*T1 + 10*C1 - 4*Math.pow(C1, 2) - 9*e2)*Math.pow(D, 4)/24.0
Dim Q4 As Double = (61 + 90*T1 + 298*C1 + 45*Math.pow(T1, 2) - 3*Math.pow(C1, 2) - 252*e2)*Math.pow(D, 6)/720.0
Dim lat As Double = fp - Q1*(Q2 - Q3 + Q4)
' // long
Dim Q5 As Double = D
Dim Q6 As Double = (1 + 2*T1 + C1)*Math.pow(D, 3)/6
Dim Q7 As Double = (5 - 2*C1 + 28*T1 - 3*Math.pow(C1, 2) + 8*e2 + 24*Math.pow(T1, 2))*Math.pow(D, 5)/120.0
Dim lon As Double = lon0 + (Q5 - Q6 + Q7)/Math.cos(fp)
Console.WriteLine("lat:"+lat.ToString)
Console.WriteLine("lon:"+lon.ToString)
Dim latitude As Double = radians_to_degrees(lat)
Dim longitude As Double = radians_to_degrees(lon)
Console.WriteLine("緯度latitude:"+latitude.ToString)
Console.WriteLine("經度longitude:"+longitude.ToString)
Dim sets(5) As Double
sets(1) = latitude
sets(2) = longitude
convertToWGS = (sets)
End Function
' 用到的參考
' https://theprogrammingexpert.com/vba-radians-to-degrees/
Function radians_to_degrees(radians as Double) as Double
radians_to_degrees = radians * (180 / Math.Pi)
End Function
' //初始化計算用相關參數
Sub prepareParam()
' //初始化計算用相關參數
e = Math.pow((1- Math.pow(b,2)/ Math.pow(a,2)), 0.5)
e1 = (1.0 - Math.pow((1.0 - Math.pow(e, 2)), 0.5)) / (1.0 + Math.pow((1.0 - Math.pow(e, 2)), 0.5))
e2 = Math.pow((e*a/b), 2)
J1 = (3*e1/2 - 27*Math.pow(e1, 3)/32.0)
J2 = (21*Math.pow(e1, 2)/16 - 55*Math.pow(e1, 4)/32.0)
J3 = (151*Math.pow(e1, 3)/96.0)
J4 = (1097*Math.pow(e1, 4)/512.0)
End Sub
Sub printPrepareParam()
Console.WriteLine("dx:"+dx.ToString)
Console.WriteLine("dy:"+dy.ToString)
Console.WriteLine("lon0:"+lon0.ToString)
Console.WriteLine("k0:"+k0.ToString)
Console.WriteLine("a:"+a.ToString)
Console.WriteLine("b:"+b.ToString)
Console.WriteLine("e:"+e.ToString)
Console.WriteLine("e1:"+e1.ToString)
Console.WriteLine("e2:"+e2.ToString)
Console.WriteLine("J1:"+J1.ToString)
Console.WriteLine("J2:"+J2.ToString)
Console.WriteLine("J3:"+J3.ToString)
Console.WriteLine("J4:"+J4.ToString)
End Sub
End Module
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment