直角平面座標CSVファイルから緯度経度CSVファイルへの一括変換プログラム(VBS)

　csvで作成された大量の直角座標XYを緯度経度に変換する必要があり、国土地理院の「緯度、経度への換算」サイトの一括変換も考えたたが、geojsonへ変換することもあり,Windowsなら他に何もプログラムがなくても動くようにvbsで自作することにしました。
　自作といっても国土地理院の「緯度、経度への換算」サイトの計算式をVBSで書き直しただけですが。。。

ここで公開していますのでお使いください。
【使用方法】
・XY2LatLng.zipを解凍するとXY2LatLng.txtができる
・XY2LatLng.txtをXY2LatLng.vbsに名前変更
・XY2LatLng.vbsにcsvファイルDrag&Dropするだけ
・変換元csvがあったフォルダーに緯度経度に変換された"元ファイル名_LatLng.csv"
　ファイルが作成されます。

【ソース XYLatLng.vbs】

Option Explicit
Dim args, arg, k, FName

Set args = WScript.Arguments
For Each arg In args
	Call MsgBox(arg & "を処理中")
	Do
		k = InputBox("平面直角座標系:" & vbCr & "系番号入力(1から19)")
	Loop While (k < 1 or k > 19)
	
	genLatLng(arg)
	
	Call MsgBox(FName & "_LatLng.csv" & " 書込み完了")
Next

Sub genLatLng(arg)
	Dim fso, pos, PathName
	Set fso = WScript.CreateObject("Scripting.FileSystemObject")

	Dim inputFile
	Set inputFile = fso.OpenTextFile(arg, 1, False, 0)
	
	' パス、ファイル名取得
	pos = InStrRev(arg,".")
	PathName = Left(arg, pos - 1)
	FName = Mid(arg, InStrRev(arg,"\") + 1, pos - InStrRev(arg,"\") - 1)
	
	' 出力ファイル
	Dim outputFile
	Set outputFile = fso.OpenTextFile(PathName & "_LatLng.csv", 2, True)

	Do Until inputFile.AtEndOfStream
		Dim lineStr,aryStrings
		lineStr = inputFile.ReadLine
		aryStrings = Split(lineStr, ",")
		If aryStrings(0) = "" Then
			Exit Do
		End If

		Dim latLngStr
		latLngStr = xy2LatLng(k, aryStrings(0), aryStrings(1))
		outputFile.WriteLine latLngStr
	Loop
	
	inputFile.Close
	outputFile.Close

End Sub

Function xy2LatLng(k, x, y)
	Dim b0, L0
	Select Case k
		case  1
			b0=33.0
			L0=129.5
		case  2
			b0=33.0
			L0=131.0
		case  3
			b0=36.0
			L0=132.1666666666666
		case  4
			b0=33.0
			L0=133.5
		case  5
			b0=36.0
			L0=134.3333333333333
		case  6
			b0=36.0
			L0=136.0
		case  7
			b0=36.0
			L0=137.1666666666666
		case  8
			b0=36.0
			L0=138.5
		case  9
			b0=36.0
			L0=139.833333333333
		case 10
			b0=40.0
			L0=140.833333333333
		case 11
			b0=44.0
			L0=140.25
		case 12
			b0=44.0
			L0=142.25
		case 13
			b0=44.0
			L0=144.25
		case 14
			b0=26.0
			L0=142.0
		case 15
			b0=26.0
			L0=127.5
		case 16
			b0=26.0
			L0=124.0
		case 17
			b0=26.0
			L0=131.0
		case 18
			b0=20.0
			L0=136.0
		case 19
			b0=26.0
			L0=154.0
	End Select

  Const a  = 6378137.0             ' 長半径
  Const F  = 298.257222101         ' 逆扁平率
  Const m0 = 0.9999                ' 平面直角座標系のX軸上における縮尺係数
	Const Pi = 3.1415926535897932

	' 緯度経度をラジアン変換
	Dim radLat0, radLng0
	radLat0 = deg2rad(b0)            ' 座標原点の緯度
	radLng0 = deg2rad(L0)            ' 座標原点の経度

	Dim n
	n  = 1/(2*F-1)
	
	Dim beta1, beta2, beta3, beta4, beta5
	beta1 = n/2 - 2/3*(n^2.0) + 37/96*(n^3.0) - 1/360*(n^4.0) - 81/512*(n^5.0)
	beta2 = 1/48*(n^2.0) + 1/15*(n^3.0) - 437/1440*(n^4.0) + 46/105*(n^5.0)
	beta3 = 17/480*n^3.0 - 37/840*(n^4.0) - 209/4480*(n^5.0)
	beta4 = 4397/161280*(n^4.0) - 11/504*(n^5.0)
	beta5 = 4583/161280*(n^5.0)

	Dim delta1, delta2, delta3, delta4, delta5, delta6
	delta1 = 2*n - 2/3*(n^2.0) - 2*(n^3.0) + 116/45*(n^4.0) + 26/45*(n^5.0) - 2854/675*(n^6.0)
	delta2 = 7/3*(n^2.0) - 8/5*(n^3.0) - 227/45*(n^4.0) + 2704/315*(n^5.0) + 2323/945*(n^6.0)
	delta3 = 56/15*(n^3.0) - 136/35*(n^4.0) - 1262/105*(n^5.0) + 73814/2835*(n^6.0)
	delta4 = 4279/630*(n^4.0) - 332/35*(n^5.0) - 399572/14175*(n^6.0)
	delta5 = 4174/315*(n^5.0) - 144838/6237*(n^6.0)
	delta6 = 601676/22275*(n^6.0)
	
	Dim A0, A1, A2, A3, A4, A5, A6
	A0 = 1 + n^2.0/4 + (n^4.0)/64
	A1 = -3/2*(n - n^3.0/8 - (n^5.0)/64)
	A2 = 15/16*(n^2.0 - (n^4.0)/4)
	A3 = -35/48*(n^3.0 - 5/16*(n^5.0))
	A4 = 315/512*(n^4.0)
	A5 = -693/1280*(n^5.0)

	Dim Sbar, Abar, xi, eta
	Sbar = (m0*a/(1 + n))*(A0*radLat0 + (A1*Sin(2*1*radLat0) + A2*Sin(2*2*radLat0)) + A3*Sin(2*3*radLat0) + A4*Sin(2*4*radLat0) + A5*Sin(2*5*radLat0))
	Abar = (m0*a/(1 + n))*A0
	xi   = (x + Sbar)/Abar
	eta  = y/Abar
	
	Dim xi_d, eta_d, chi, lat, lng
	xi_d = xi - (beta1*Sin(2*1*xi)*cosh(2*1*eta) + beta2*Sin(2*2*xi)*cosh(2*2*eta) + beta3*Sin(2*3*xi)*cosh(2*3*eta) + beta4*Sin(2*4*xi)*cosh(2*4*eta) + beta5*Sin(2*5*xi)*cosh(2*5*eta))
	eta_d = eta - (beta1*Cos(2*1*xi)*sinh(2*1*eta) + beta2*Cos(2*2*xi)*sinh(2*2*eta) + beta3*Cos(2*3*xi)*sinh(2*3*eta) + beta4*Cos(2*4*xi)*sinh(2*4*eta) + beta5*Cos(2*5*xi)*sinh(2*5*eta))
	chi = asin(Sin(xi_d)/cosh(eta_d))

	lat = chi + (delta1*Sin(2*1*chi) + delta2*sin(2*2*chi) + delta3*Sin(2*3*chi) + delta4*Sin(2*4*chi) + delta5*Sin(2*5*chi) + delta6*Sin(2*6*chi))
	lng = radLng0 + Atn(sinh(eta_d)/Cos(xi_d))
	
	lat = lat * (180/Pi)
	lng = lng * (180/Pi)
	xy2LatLng = lat & "," & lng
End Function

Function deg2rad(x)
	Const Pi = 3.1415926535897932
	deg2rad = (x * Pi)/180
End Function

Function asin(x)
	asin =  Atn(X / Sqr(-X * X + 1))
End Function

Function cosh(x)
	cosh = (Exp(X) + Exp(-X)) / 2
End Function

Function sinh(x)
	sinh = (Exp(X) - Exp(-X)) / 2
End Function

コメント 2

こちらのコードを利用したいのですが著作権をご教示願います。
by カトウダイジロウ (2022-01-18 21:59) 

変更・使用を含め全てフリーです。
VBAにはパスワードをかけていませんのでご自由に変更してください。
その他のExcel便利ツールを勤務先のHP
https://coresys.co.jp/system_rpa/excel/
で公開しています。ご利用ください。
by tyama0467 (2022-01-25 11:41)