擷取全港三角網測站及導線點 WGS84 座標列表

Sat, 07 Apr 2012 14:17 警告及聲明
以下資料有已知誤差，只作一般參考用途，而非定位測量指引。本人並不保証以下內容的真確性及安全程度，如引致任何意外或任何問題，本人並無責任承擔。

方法
我想製作一個可配合 GPS 接收機所顯示之 WGS84 經緯度的全港三角網測站及導線點列表，以作一些應用，所以我嘗試找這類列表。在找之前，我聯絡過地政總署，但對方稱未能提供，那惟有自己搜集。

於地政總署測繪處大地測量網站，公眾可下載每一個三角網測站及導線點的紀錄。它可以以該站的編號或在地圖上進行搜尋。

但如此的方法，我並不容易擷取列表，原因是：
一、該處並沒提供完整列表
二、網頁上並沒有記錄某站的 WGS84 或 HK1980 之座標，所有紀錄於 PDF 檔案上出現，但該等紀錄並不一定以文字形式儲存，不易以程式擷取。

為了應付以上問題，我嘗試過以下方案；
首先將每個列有每一站編號的 HTML 拆解出來，例如這一頁。
再按編號，下載所有 PDF 檔案，例如這個。
然後再用 ImageMacgick 的 convert 將所有 PDF 檔案轉換成 tiff 檔案，再以 tesseract 將圖像辨識成文字，之後再用 grep 和 sed 製作列表。

但這個方法有著一些問題：
一、不是所有 PDF 檔案都是標準格式而可以轉換成 tiff 檔案。
二、就算能夠轉換，都不一定成功辨識成文字。
三、圖像辨識仍需肉眼作檢查，難以短時間內完成。

故此我放棄了這方案。之後我發現其地圖搜索系統的 HTML 有 Image Map 的座標，讓用戶點選的時候，可以提供 PDF 檔案的連結。我利用這些 HTML Image Map 座標，再以兩點已知的站座標，整合起來，以一次多元方程計算，可以換算出找出原點的座標、每一像素對應的距離，引伸就可以網址的格式和 HTML Image Map 座標找出每一站 HK1980 座標。

這種方案完全避開了辨識 PDF 檔案的部分，但問題在於我並不是取用其檔案內的數據，而且經過計算，所以存有一定誤差，我並沒有計算誤差的量，但大概為兩米。

這種換算要以它提出最精細的版本（例子），以得出較精確的結果。由於它每一幅圖之間有重疊的部分，所以寫程式擷取時需作考慮。

我寫了以下 Python 程式，它可自動擷取包含全港各站的地圖，分析其 HTML 檔案，計算出其 HK1980 座標，再儲存成一個包括所有座標的檔案。

#!/usr/bin/env python

import urllib
import re

pixelLengthRatioN = 1.7700131
pixelLengthRatioE = 1.76896
topLeftN = 848148.741
topLeftE = 799809.53356

pixelOffsetN = 508-169
pixelOffsetE = 638-214

jSequence = ['NW','NE','SW','SE']
lSequence = ['A','B','C','D']

resultList = {}

for i in range(2, 17):
	i_n = int((i-1)/4)
	i_e = int((i-1)%4)
	for j in jSequence:
		j_n = int((jSequence.index(j)) / 2)
		j_e = int((jSequence.index(j)) % 2)
		for k in range(1,26):
			k_n = int((k-1)/5)
			k_e = int((k-1)%5)
			for l in lSequence:
				l_n = int((lSequence.index(l)) / 2)
				l_e = int((lSequence.index(l)) % 2)
				urlSegment = str(i)+"-"+j+"-"+str(k)+l
				print urlSegment
			
				content = ""
				try:
					f = urllib.urlopen("http://www.geodetic.gov.hk/jpgkey/1/HTM/" + urlSegment + "_tr.asp?L=T")
					content = f.read()
					f.close()
				except Exception:
					content = ""
				
				
				for line in content.splitlines():
					if "/jpgkey/GetSummary2.asp?stn=" in line:
						m = re.search(r"coords='([0-9]*)\,([0-9]*)\,.*stn=(.*)\&this=", line)
						if m:
							point = m.group(3)
							pointE = m.group(1)
							pointN = m.group(2)
							offsetN = (i_n*pixelOffsetN*20 + j_n*pixelOffsetN*10 + k_n*pixelOffsetN*2 + l_n*pixelOffsetN + int(pointN))* pixelLengthRatioN
							offsetE = (i_e*pixelOffsetE*20 + j_e*pixelOffsetE*10 + k_e*pixelOffsetE*2 + l_e*pixelOffsetE + int(pointE))* pixelLengthRatioE
							coordinatesN = topLeftN - offsetN
							coordinatesE = topLeftE + offsetE
							print point + " N: " + str(coordinatesN) + " E: " + str(coordinatesE)
							resultList[point] = (coordinatesN, coordinatesE)
							
f = open('result.txt', 'w')							
for key, value in resultList.iteritems():
	f.write(key + ';' + str(value[0]) + ';' + str(value[1]) + '\n')
f.close()

執行此 Python 程式後再於 Shell 內執行以下指令以排序：

cat result.txt | sort > resultSorted.txt

打開 resultSorted.txt，可見以下的格式：

該三角網測站或導線點的編號;N;E

例子為：

1120.02;815922.112488;816722.56012
1120.03;815442.438938;816858.77004
119;815987.602973;840920.16396
124;821844.576321;837030.22092
128;817203.601973;843168.51212
129;819042.645584;841549.91372

由於這只是 HK1980 的格網座標，而非 UTM 座標或經緯度，亦非 WGS84 的 UTM 座標或經緯度，並不能配合以 WGS84 為基準的 GPS 接收器運作，故需要先進行轉換。

香港大地測量基準說明及 Schematic Diagram Showing Transformation of Coordinates between WGS84 Geographic Coordinates and HK 1980 Grid Coordinate 提供了轉換的方式及說明。這個方式似是複雜，我在 github 找到了轉換的函數，我就引入使用。

以下程式作者以 GNU AFFERO GENERAL PUBLIC LICENSE 釋出：

"""
   Convertion between HK1980 grid reference system and WGS84
 
   (Ported from the PHP code write by Abel Cheung)
"""
import math

M_PI_180 = 0.0174532925199432957692369076849

def wgs84_to_hk1980 (lat,lon,mode=2):
	"""
	  Converts HK1980 data to WGS84 data
	 
	  @param int mode 1 = Hayford, 2 = WGS84
	  @return A tuple of HK1980 grid reference system (North , East )
	"""

	assert(mode == 1 or mode == 2 )

	assert( lat >= 22 and lat < 23)

	assert( lon >= 113 and lon < 115) 

	if mode == 1:
		d4 = 22.312133333333335 * M_PI_180
		d5 = 114.17855555555556 * M_PI_180
	else:
		d4 = 22.310602777777778 * M_PI_180;
		d5 = 114.1810138888889 * M_PI_180;

	d6 = lat * M_PI_180;
	d7 = lon * M_PI_180;

	d8 = SMER (mode, 0, d4);
	d9 = SMER (mode, 0, d6);

	(d34, d35) = RADIUS (mode, d6);

	d10 = (d7 - d5) * math.cos (d6);
	d11 = math.tan (d6);
	d12 = pow (d11, 2);
	d13 = pow (d11, 4);
	d14 = pow (d11, 6);

	d15 = d35 / d34;
	d16 = pow (d15, 2);
	d17 = pow (d15, 3);
	d18 = pow (d15, 4);

	d19 = d9 - d8;

	d20 = (d35 / 2) * pow (d10, 2) * d11;
	d21 = (d20 / 12) * pow (d10, 2) * (d16 * 4 + d15 - d12);
	d22 = (d21 / 30) * pow (d10, 2);
	d22 *=   d18 * (88 - d12 * 192) \
		- d17 * (28 - d12 * 168) \
		+ d16 * (1 - d12 * 32) \
		- d15 * d12 * 2 \
		+ d13;
	d23 = d22 / 56 * pow (d10, 2) \
		* (1385 - 3111 * d12 + d13 * 543 - d14);
	d24 = d35 * d10;
	d25 = d24 / 6  * pow (d10, 2);
	d26 = d25 / 20 * pow (d10, 2);

	d27 = d26 / 42 * pow (d10, 2); 
	# The original line from Abel's code should be 
	#  d27 = d27 / 42 * pow (d10, 2);

	d25 *= d15 - d12;
	d26 *= d17 * 4 * (1 - d12 * 6) \
		+ d16 * (1 + d12 * 8) \
		- d15 * d12 * 2 + d13;
	d27 *= 61 - d12 * 479 + d13 * 179 - d14;

	d36 = d19 + d20 + d21 + d22 + d23 + 819069.80000000005;
	d37 = d24 + d25 + d26 + d27 + 836694.05000000005;

	if mode == 2:
		d28 = d36;
		d29 = d37;
		d30 = 0.99999983729999997;
		d31 = -2.7858E-005;

		d36 = (d28 * d30 - d29 * d31) - 23.098331000000002;
		d37 = (d28 * d31 + d29 * d30) + 23.149764999999999;

	north = d36;
	east = d37;

	return (north , east)

def hk1980_to_wgs84 (north,east,mode=2):
	"""
	 Converts HK1980 data to WGS84 data
	 
	 @param mode 1 = Hayford, 2 = WGS84
	 @return A tuple of (lat,lon)
	"""
	assert(mode == 1 or mode == 2 )	

	d1 = north;
	d2 = east;

	if mode == 1:
		d4 = 22.312133333333335 * M_PI_180
		d5 = 114.17855555555556 * M_PI_180
	else:
		d4 = 22.310602777777778 * M_PI_180;
		d5 = 114.1810138888889 * M_PI_180;

	if (mode == 2):
		d6 = 1.0000001619000001;
		d7 = 2.7858E-05;
		d8 = 23.098979;
		d9 = -23.149125000000002;
		d10 = d6 * d1 - d7 * d2 + d8;
		d11 = d7 * d1 + d6 * d2 + d9;
		d1 = d10;
		d2 = d11;

	d28 = d1 - 819069.80000000005;
	d37 = d2 - 836694.05000000005;
	d12 = 6.8535615239999998;
	d13 = 110736.3925;
	d14 = (math.sqrt (d28 * d12 * 4 + pow (d13, 2)) - d13) / 2 / d12 * M_PI_180;
	d15 = d4 + d14;
	d16 = 0;
	d18 = 0;

	while True:
		d15 += d16;
		d17 = SMER (mode, d4, d15);
		d18 = d28 - d17;
		d39, d41 = RADIUS (mode, d15);
		d16 = d18 / d39;
		if not (abs(d18) > 1E-06):
			break;

	(d40, d42) = RADIUS (mode, d15);
	d19 = math.tan (d15);
	d20 = pow (d19, 2);
	d21 = pow (d19, 4);
	d23 = d42 / d40;
	d24 = pow (d23, 2);
	d25 = pow (d23, 3);
	d26 = pow (d23, 4);
	d28 = d37 / d42;
	d38 = pow (d28, 2);
	d29 = d37 / d40 * d28 * d19 / 2;
	d30 = d29 / 12 * d38 * ((9 * d23 * (1 - d20) - 4 * d24) + 12 * d20);
	d31 = d29 / 360 * pow (d38, 2);
	d31 *= 8 * d26 * (11 - 24 * d20) \
		- 12 * d25 * (21 - 71 * d20) \
		+ 15 * d24 * (15 - 98 * d20 + 15 * d21) \
		+ 180 * d23 * (5 * d20 - 3 * d21) \
		+ 360 * d21;
	d32 = d29 / 20160 * pow (d38, 3) * (1385 + 3633 * d20 + 4095 * d21 + 1575 * d20 * d21);
	d44 = d15 - d29 + d30 - d31 + d32;
	d33 = d28 / math.cos (d15);
	d34 = d33 * d38 / 6 * (d23 + 2 * d20);
	d35 = d33 * pow (d38, 2) / 120;
	d35 *= d24 * (9 - 68 * d20) - 4 * d25 * (1 - 6 * d20) + 72 * d23 * d20 + 24 * d21;
	d36 = d33 * pow (d38, 3) / 5040 * (61 + 662 * d20 + 1320 * d21 + 720 * d20 * d21);
	d43 = d5 + d33 - d34 + d35 - d36;
	d44 /= M_PI_180;
	d43 /= M_PI_180;

	lat = d44;
	lon = d43;

	return (lat,lon)

def SMER (mode, d1, d2) :
	"""
	 * @param int mode
	 * @param float $d1
	 * @param float $d2
	"""
	if mode == 1:
		d3 = 6378388
		d4 = 0.0033670033670033669
	else:
		d3 = 6378137
		d4 = 0.0033528106647429845

	d5 = 2.0 * d4 - pow (d4, 2);

	d6 = 1.0 \
		+ 0.75 * d5 \
		+ 0.703125 * pow (d5, 2) \
		+ 0.68359375 * pow (d5, 3);
	d7 = 0.75 * d5 \
		+ 0.9375 * pow (d5, 2) \
		+ 1.025390625 * pow (d5, 3);
	d8 = 0.234375 * pow (d5, 2) \
		+ 0.41015625 * pow (d5, 3);
	d9 = 0.068359375 * pow (d5, 3);

	d10 = d2 - d1;
	d11 = math.sin (d2 * 2) - math.sin (d1 * 2);
	d12 = math.sin (d2 * 4) - math.sin (d1 * 4);
	d13 = math.sin (d2 * 6) - math.sin (d1 * 6);

	d14 = d3 * (1.0 - d5);
	d14 *= (d6 * d10) \
		- (d7 * d11 / 2.0) \
		+ (d8 * d12 / 4.0) \
		- (d9 * d13 / 6.0);

	return d14;


def RADIUS (mode, d1):
	"""
	  @access private
	  @param mode
	  @param d1
	 """
	if (mode == 1):
		d2 = 6378388
		d3 = 0.0033670033670033669
	else:
		d2 = 6378137
		d3 = 0.0033528106647429845

	d4 = d3 * 2 - pow (d3, 2);
	d5 = 1.0 - d4 * pow (math.sin (d1), 2);
	d6 = (d2 * (1.0 - d4)) / pow (d5, 1.5);
	d7 = d2 / math.sqrt (d5);

	return (d6, d7);


if __name__ == "__main__":
	hk1980 = (838477.970,818097.267)

	wgs84 = hk1980_to_wgs84(*hk1980)
	print "%s => %s"  % (hk1980 , wgs84)

	print "%s => %s"  % (wgs84 , wgs84_to_hk1980(*wgs84))

而我則寫了以下程式配合使用，由於涉及轉換，誤差將增大：

#!/usr/bin/env python

from HK1980 import *

f = open('resultSorted.txt', 'r')
f2 = open('resultSortedConverted.txt', 'w')
for line in f:
	point, strN, strE = line.replace("\n", "").split(";")
	N = float(strN)
	E = float(strE)
	
	converted = hk1980_to_wgs84 (N,E,2)
	f2.write(point + ',' + strN + ',' + strE + ',' + str(converted[0]) + ',' + str(converted[1]) + '\n')
f.close()
f2.close()

它會輸出 resultSortedConverted.txt，格式為：

該三角網測站或導線點的編號,N,E,WGS84 北緯度,WGS84 東經度

例子為：

1120.02,815922.112488,816722.56012,22.2820658608,113.987220081
1120.03,815442.438938,816858.77004,22.2777356611,113.988547571
119,815987.602973,840920.16396,22.2827616686,114.222019021
124,821844.576321,837030.22092,22.3356598535,114.184275982
128,817203.601973,843168.51212,22.2937353657,114.243840972
129,819042.645584,841549.91372,22.3103487937,114.228139838

而 WGS84 北緯度及WGS84 東經度則可直接配合 GPS 接收器使用。

是次轉換的結果可在此下載：
下載網址

本文連結