標題 座標轉換 時間 Wed Oct 21 02:53:22 2009 ─────────────────────────────────────── txt 編輯軟體： Crimson Editor 本文用到了 big5 符號表。
*排版亂掉了
*所有來源都失效了
*(U)TM轉換應可採用:
https://github.com/jjimenezshaw/Leaflet.UTM
學習 電力座標 到 座標轉換 0. a, b , 1/f , e^2 , llh , ECEF-XYZ a: Semi-major axis 參考橢球長半徑 (長軸/短軸) b: Polar Radius 參考橢球短半徑, 兩極半徑 1/f: 扁率倒數 e^2: 偏心率平方 first eccentricity squared e is the eccentricity of the earth's elliptical cross-section 偏心率 地球 橢圓狀 橫截面 (from uwgb) (Geodetic) LLh : Latitude , Longitude , height (Ellipsoid) 經度 , 緯度 , 橢球高 phi,φ , Lambda λ 單位: degree , 度 橢球高相關資料： 橢球高 h: Ellipsoid Height 大地水準面高、大地起伏 N: Geoid(al) Height / Geoid(al) undulation 正高、水準高 H: Orthometric Height h: Ellipsoid --along Ellipsoid normal--> Topography (Earth's Surface) N: Ellipsoid --along Ellipsoid normal--> Geoid H: Geoid --Plumb Line--> Topography ＊三條都是直線。其中 H 是垂直水準面到真實地表 h 是垂直橢球面到真實地表(GPS) ，N 是垂直橢球面到水準面。 ＊故 h,N 不一定在 H 上，h 與 N 的延線也未必在同一條線上（？） 固定公式為： h = H + N 以七星山為例，計算 by PCTrans 4.2.5 EGM 96 Geoid source:NGA 地球重力模型1996 大地水準面 來源: National Geospatial-Intelligence Agency （美國國家地理空間情報局） Topography ----↑ ↑ Geoid ----│ └ H=1120.36 Ellipsoid -----└ h=1139.589(gps) └ N=19.23 ----------------------------------------------------- ECEF-XYZ : Earth Centered, Earth Fixed Cartesian coordinates X Y Z 地 心 地 固 卡氏 座標系 X_Prime Meridian : 本初子午線 經過 英國/格林威治 Y_equator : 赤道 單位: m , 公尺 Z | |___ Y (y=0, x=90) / / X (=0) ＊地心地固參考面，適合用來表示空間中的位置和向量。 ------------------------------------------------------ 計算時全部轉成徑度量。 radian(徑"度"量) = (弧度"秒"/3600)×PI()/180 換成度 ; 3.1415926... 另， degree(弧度度量) = radian*180/PI() ------------------------------------------------------ 1. LLh -> xyz 需要準備的有... Lat1, Long1 , h1 -> rLat1 , rLong1 (a,b) 或 (a,1/f) 1/f = a/(a-b) ; b = a*(1-f) ＊所以你必須先知道：待處理的[經緯度]所使用的[參考橢球]為何。 ------------------------------------------------------- 再算得: e^2 ; chi , N(lat) e^2 = 2*f-f^2 = 1- b^2 / a^2 (numerical eccentricity) ∵ ellipse , first eccentricity = SQRT(a^2 - b^2) / a ∴ squared = (a^2 - b^2) / a^2 , (a^2/a^2 提出) , 1 - b^2 / a^2 χ,chi = SQRT( 1 - (e^2)*(SIN(rLat1)^2) ) ν, N(lat) = a / chi , Normal , Radius of curvature in prime vertical 垂線 半徑 曲率 卯酉圈 ＊從橢球表面作垂線，然後交於地心地固圓球Z軸的距離。 ＊比如說我站在地表某點，往四周水平處看即地平面，與地平面垂直者當然就是卯酉圈。 ＊兩個講的是一樣的東西。 ---------------------------------------------------------- 帶入公式： lat <-> phi,φ ; long <-> Lambda, λ ; nu,ν X1 = (N + h)*COS(rLat1)*COS(rLong1) Y1 = (N + h)*COS(rLat1)*SIN(rLong1) Z1 = ( (N)*(1 - e^2) + h )*SIN(rLat1) ---------------------------------------------------------- 各地經緯度座標，轉換至WGS84基準橢球經緯度座標的過程，以七參數為例： 選擇參考橢球(1) -> 輸入 LLh(1) -> 轉成 XYZ(1) -> 代入七參數轉換引擎 使用參考橢球(1) 轉換至 WGS84 零度基準所需的七參數 dx dy dz rx ry rz s ＊以 PCTrans 為例，其中所列的轉換參數都是相對於 WGS84 基準。 ＊也就是 dx dy dz = 0 m , rx ry rz = 0 sec , s = 0 ppm ＊雖然是荷蘭海軍的程式，採用的轉換慣例卻是 Coordinate Frame Rotation (by Help) 得 XYZ(0)，基準是 WGS84。 然後，再從 WGS84基準，轉換至當地新基準橢球經緯度。 輸入先前 XYZ(0) 選定目標地參考橢球(2)，查詢轉換至目的地橢球基準的七參數 ＊查表都是查得 "當地轉WGS84" 所需的七參數，這邊加工一下，通通乘(-1) ＊其中注意尺度系數 s：是單位為 ppm 時乘(-1) 最後得 XYZ(2)，目標地的橢球基準，就可轉為 LLh(2) 整理： LLh(1) -> XYZ(1) -> XYZ(0) -> XYZ(2) -> LLh(2) ------------------------------------------------------ 七參數轉換：這裡採用 Bursa-Wolf s:source,1 t: target,2 r: rotation,ε s: scale factor (ppm) ,κ 目的 平移 尺度 旋轉 來源 [X2] [dX] [ 1 rz -ry ] [X1] [Y2] = [dY] + (1+s)．[-rz 1 rx ]．[Y1] [Z2] [dZ] [ ry -rx 1 ] [Z1] ＊POSC 的圖八成錯了。因為旋轉矩陣的第一排是 rz rx ，第二排卻也是 rz rx 可能是第一排最右的 r"x" 是 r"y" 誤標的。 http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs35.html http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/helm1.gif ＊注意這裡公式的旋轉矩陣慣例，採用的是 Coordinate Frame Rotation(CFR) 像是美系探勘生產工業。 歐系則是 Position Vector Transformation (PVR) ＊Position Vector Transformation ：從各軸原點，沿軸正方向去看，是順時針旋轉。 ╱ Z + ╱ ╱ ↑<╮ ╱ 右手定則：四指彎曲方向為旋轉正值。 ∠__╰>│___ ╱ │ ｜ positive end , 搭便車(hitch-hiking pose)大拇指所在方位 ├───→+ Y ╭╭╭╮ ╭> ╮ ↙ ｜ ↓ ↓ ↓ ↑╭→ ↑⊙↓ 假設有顆子彈沿著Y軸正方向跑... ----- -- ｜ ╯ ╯ ╰ <╯ 看到子彈的屁股。
+↙ | ｜ ╱ X | ｜╱
| ---------| 從三軸來看其旋轉面： Position Vector Rotation (PVR)：各軸右手定則遵循方向為 PVR 正值。 Coordinate Frame Rotation (CFR)： 1st 2nd 3rd dφ/dt dθ/dt dψ/dt Z 0,1 Y 1,2 X 2,3 ↑<╮ 。 <-Taiwan-> 。 ↑<╮ ↑<╮ ╰>│ ╰>│ ╰>│ ├──→ Y 0-1,2 ├──→ X 1-2,3 ├──→ Z 2-3 ╱ <╯ ╱ <╯ ╱ <╯ ↙ φ+ ↙ θ+ ↙ ψ+ X 0-1 Z 0,1-2 Y 1,2-3 。 ↖台灣大概在這 (sequence X=1 Y=2 Z=3 ; 3-2-1 ; 0= 1=' 2='' 3=''') ＊ CFR 負值方向就跟右手定則方向一樣了。 ＊如果看到的公式，其中的旋轉矩陣符號如此者，則是 PVR [X2] [dX] [ 1 -rz ry ] [X1] [Y2] = [dY] + (1+s)．[ rz 1 -rx ]．[Y1] [Z2] [dZ] [-ry rx 1 ] [Z1] ＊差別在於，在使用七參數時，rx ry rz 差一負號。 ＊在 PVR 為正值，到了 CFR 為負值。反之亦然。 如果使用上述 PVR 慣例的轉換公式 以 TWD67 7para 為例：如果查得的參數是 dx = -730.16 rx = -7.968 rx = 7.968 dy = -346.212 ry = -3.5498 ----------> ry = 3.5498 dz = -472.186 rz = -0.4063 rz = 0.4063 已知此參數慣例為 CFR ，要用於 PVC 慣例的公式，則 rx ry rz 乘 (-1) 如果不知道手上的參數慣例是 PVR 還是 CFR ？（例如 Transformaton NATO.pdf :P ） 1.查閱標準技術資料 2.利用已有控制點查核 3. 試誤 七參數建議得註明所使用的演算法，轉換來源與目標，與轉換慣例！ ------------------------------------------------------ [X2] [dX] [ 1 rz -ry ] [X1] [Y2] = [dY] + (1+s)．[-rz 1 rx ]．[Y1] [Z2] [dZ] [ ry -rx 1 ] [Z1] 首先先解釋一下 (1+s) 這部份： s 單位為 ppm ，一般可見公式寫為 (1+ds) 或寫在矩陣內對角線上。[R'] ＊這是個很長的推導... @@"。簡單的說 [R] = [U]+[I] （[U] 對角線皆為 0 ， [I] 為單位矩陣對角線皆為 1 ，其餘皆為 0） ＊個人想法是，(1+ds)‧([U] + [I]) = (1)‧[U] + ds‧[U] + 1‧[I] + ds‧[I]，其中 ds‧[U] 可省略。（？） ，最後 [U] + (1+ds)‧[I] = [R'] 得再換算為 scale factor M = (1 + s*10^-6 ) ，有些稱為 k ＊ scale factor 皆為正值。當 ppm 為 0 時則 scale 為 1 ； ppm scale -1000 0.999 -100 0.9999 -10 0.99999 -1 0.999999 +10 1.00001 +100 1.0001 +1000 1.001 ＊個人解釋部份： ＊如同消去 ppm 為乘以 10^-6 ，同樣的消去 scale 或許也是 -1 (100% = 1) 於是公式解得： X2 = dX + M*( (1)* X1 + (1) *Y1* rz + (-1)* Z1*ry ) Y2 = dY + M*((-1)* X1*rz + (1) *Y1 + (1) * Z1*rx ) Z2 = dZ + M*( (1)* X1*ry + (-1)*Y1*rx + (1) * Z1 ) 還記得這個嗎？ LLh(1) -> XYZ(1) -> XYZ(0) -> XYZ(2) -> LLh(2) XYZ(1) -> XYZ(0) 所使用的是 "轉換到 WGS84" 的參數 XYZ(0) -> XYZ(2) 已經是 WGS84 了，轉換至目的地則七參數 乘 -1 例如：1. TWD67 轉 WGS84 2. WGS84 轉 TWD67 dX dY dZ rx(sec) ry rz ds(ppm) 1. -730.16 -346.212 -472.186 -7.968 -3.5498 -0.4063 -18.2 2. 730.16 346.212 472.186 7.968 3.5498 0.4063 18.2 注意代入公式計算時， rx ry rz 要換成徑度量，ds要再經過轉換得 1. rx ry rz k -0.00003863 -0.00001721 -0.00000197 0.9999818 2. 0.00003863 0.00001721 0.00000197 1.0000182 由於是弧度"秒"，要先換成"度"，再換成徑度量，故為 (-7.968 / 3600) *PI()/180 連續解兩次的公式則是： X0 = dX1 + M1*( (1)* X1 + (1) *Y1*rz1 + (-1)* Z1*ry1 ) Y0 = dY1 + M1*((-1)* X1*rz1 + (1) *Y1 + (1) * Z1*rx1 ) Z0 = dZ1 + M1*( (1)* X1*ry1 + (-1)*Y1*rx1 + (1) * Z1 ) X2 = dX2 + M2*( (1)* X0 + (1) *Y0*rz2 + (-1)* Z0*ry2 ) Y2 = dY2 + M2*((-1)* X0*rz2 + (1) *Y0 + (1) * Z0*rx2 ) Z2 = dZ2 + M2*( (1)* X0*ry2 + (-1)*Y0*rx2 + (1) * Z0 ) ------------------------------------------------------ 終於，得到 X2 Y2 Z2 目標地的 XYZ 了，就可以轉回 LLh 需要準備的有： X2 Y2 Z2 -> 目標地橢球定義 -> a2 , b2 or f2 -> e^2 , e'^2 這次多了二次離心率平方 (Second eccentricity) ：e'^2 e'^2 = f*(2-f)/(1-f)^2 = a^2/b^2 - 1 e^2 = 2*f-f^2 = 1- b^2 / a^2 http://en.wikipedia.org/wiki/Geodetic_system ------------------------------------------------------ by Bowring, B. 1976 http://www.colorado.edu/geography/gcraft/notes/datum/gif/xyzllh.gif atan: arctan , arc tangent the: theta , θ phi,φ=rLat , Lambda,λ = rLong , PI() = π 以下一些 2 會省略，像是 X2,Y2,Z2,a2,b2 Lambda = atan2_function( X , Y ) in Excel Form* p = r(x,y) = SQRT( X^2 + Y^2 ) the = atan( (Z*a)/(p*b) ) or = atan2_function(p*b , Z*a) in Excel Form* phi = atan( (Z + e'^2 * b * SIN(the)^3) / (p - e^2 * a * COS(the)^3) ) N(phi) = a / ( (1 - e^2 * SIN(phi)^2)^0.5 ) h = ( p / COS(phi)) - N(phi) , phi =< PI()/4 ; or = ( Z / SIN(phi)) - N(phi)*(1 - e^2) , phi > PI()/4 ＊求橢球高的公式，會因緯度而有所不同(?) ，台灣緯度小於 45 度。 ＊ atan2 function : 儘管一般論文/網頁資料上看到的是 lambda = atan2(Y,X) 注意在 excel 使用函數時，得是 atan2(X,Y) atan2 是個跟四象限有關的反三角函數，使得傳回值在 -π ～ π 之間。 (詳細見 atan2 wiki @@") ＊個人經驗上認為，似乎所有 atan 都改用 atan2 來處理比較妥當？ ＊（當然還是得考量到傳回範圍的問題。） 像是 Lambda 的 atan2 如果改用 atan 算，會有座標落於四象限的問題， 得再換算才能得到正確的經度。而用 atan2 就直接是正確的結果。(-180 ~ +180) 用 atan2 也為了解決除以 0 的問題。即使是 atan2 in Excel 也不能兩者皆為 0 (js x,y=0 atan2(y,x) 傳回 0 ,atan 傳回 NaN ; atan(y/0) = atan2(y,x) ) 理論上 x,y 為 0 會出現在 lat=90 , long=0 ,然而 LLh 轉 xyz 時 lat=90 , rlat=1.570796 , 最後 X = 0.000000000392 得到一個 ^-10 次方極小的值 如果是 X,Y 是使用者給定，寫程式時要設定判別檢查，兩者絕對值相加 仍小於 ^-10 就直接給予 long = 0.0 ＊此為逼近解法。(approximated closed-form) ＊ h 與 phi 的精度為： εh < 1.5 * 10^-3 m , εphi < 1.2 * 10^-3 m , 當 h < 400 km 時 εh < 2 * 10^-5 m , εphi < 1.5 * 10^-5 m , 當 h < 40 km 時 εh < 2 * 10^-7 m , εphi < 1.5 * 10^-7 m , 當 h < 4 km 時 ＊ so, 台灣地表物沒有大於 4000 公尺的問題...(?) ------------------------------------------------------ 電力座標轉換： 詳見：上河文化 - 大地座標系統漫談 大地座標系統與二度分帶座標 http://www.sunriver.com.tw/grid_tm2.htm 電力座標系統解讀 http://www.sunriver.com.tw/grid_taipower.htm 始祖：積丹尼 Dan Jacobson http://jidanni.org/geo/taipower/howto.html wiki OSGeo http://wiki.osgeo.org/wiki/Taiwan_Power_Company_grid 電力座標格式：例: B6030CE63 G5522AC1234 ; AB25-K3388（有些會先標後面） [分區號][圖編號x][圖編號y][百分位區號][百分位區號][十分位x][十分位y] B 60 30 C E 6 3 A~W* 00~99* 00~99 A~H A~E 0~9 0~9 [個分位x][個分位y] ＊10m 內還有其他電力設備，定位精度到 1m。 0~9 0~9 ＊單位：公尺 (m) ＊有些分區號會被省略，例如台北地區的"B"。但百分位英文區號一定會有。 ＊個人：如果只有一組四位數字，先預設為圖編號。 ＊如果看到英文標字超出上述範圍，也無法依照以下規則繼續計算。 [分區號]：每一區 80 km * 50km （8萬 * 5萬） 起點值查詢如下表： Area BottomX BottomY 90000 170000 250000 330000 2750000 A B C 2700000 D E F 2650000 G H 2600000 J K L 2550000 M N O 2500000 P Q R 2450000 T U 2400000 V W ＊excel 用 vlookup ，A欄英文編號 ，B,C欄對應 x,y ＊ 250000m 處，也就是 TWD67 二度分帶中央經線 121 度。 [圖編號x,y] 99---------99,99 每一單位為 800m * 500m | | 共切成 100x100 = 10,000 份 , 0~ 9999 | | ＊所以圖編號xx也可是視為 0000 ~ 9900 yy...▆ | xxyy 合起來就表示是該區第xxyy號圖。 | : | (循序了話yy優先。↑) 0----xx----99 每一單位▆，再分為： 4 E □□□□□□□□ [百分位圖區號A~H,A~E] 3 D □□□□□□□□ 每一單位 100m * 100m 2 C □□■□□□□□ 1 B □□□□□□□□ 0 A □□□□□□□□ A B C D E F G H 0 1 2 3 4 5 6 7 [十分位x,y]：■ [個分位x,y] 9--------9,9(*10) 9--------9,9 | | | | 5..■........|..............→ 5.... # | | : | | : | 0--3-----9 0----5---9 舉例 K7335DB2406 則計算為 K -> x = 170000 , y = 2600000 (base) x = 73 * 800 y = 35 * 500 D,B -> x = 3*100 , y = 1*100 2,4 -> x = 2*10 , y = 4*10 0,6 -> x = 0 , y = 6 sum(x) = 170000 + 73*800 + 3*100 + 2*10 + 0 = 228720 sum(y) = 2600000 + 35*500 + 1*100 + 4*10 + 6 = 2617646 ----------------------------------------------------------------- 台灣本島反算電力座標： ＊有時候走在路上，忘了圖號大區要對應多少XY？ 即使記得， 一定要換算成平面座標X6位Y7位又顯得過於麻煩。 於是可以將想要前往的目的地先查出 TWD67-TM2 平面座標，再換算回電力座標 出門後只要記得目的地電力座標，沿途看到電杆電箱上的座標便可得知相對位置關係。 例如，目的地已知是：K7335DB16，路過一個電杆號為 K7336FC07 可知大概還要再往南走約 500m 再往西走約 200m 其中主要看 K7336 與 K7335 ，y 差 1 為 500m ，DB FC 分別換為 (3,1)(5,2) x 差 2 = 200m 或是習慣看成 (D1,B6),(F0,C7) 可換成 (310,160),(500,270) 例如使用 GoogleMap , GoogleEarth 定位得 WGS84 經緯座標 透過一些轉換程式得 TWD67 經緯再轉 TWD67-TM2 X,Y 平面座標。 單位是 m ，只取 X,Y 座標整數部份。 ＊-- 直接分區判別，非台灣區不適合以下經驗公式-------------------------- Y 方向剛好是 五萬 的倍數。
根據之前的起點值查詢列表，反查得所在分區 Area
＊excel 可用 index(分區表,match(Y,Y列),match(X,X列)) 的方式來求得。
＊--------------------20100816----------------------------------------
*2021.02.04:
大分區為 X*Y = 8萬*5萬
8的倍數: 0,8,16,24,32,40
由於中央經線121 為 250000m
故可說X方向分區起點向東偏移 10000m
(ex. B分區X起點是 240000+10000=250000)
故
X商 = INT(X/80000)，分區起點Area.X = X商*80000+10000
然而加上偏移量會有大於 X的情形，則 (X商-1)*80000+10000
例如分區中X再10000就要跨下一區 ex. 319999(剩10001),320000(剩10000)
INT(319999/80000)=3, 3*80000+10000=250000
但INT(320000/80000)=4, 4*80000+10000=330000 就超過 320000 了
Y商 = INT(Y/50000)，分區起點Area.Y = Y商*50000
x_mod = X - area.X
y_mod = Y - area.Y (落點座標 - 所在分區起點值 X,Y)
剩餘部份做商餘商餘(800 -> 100 -> 10 -> 1)
800xx , 500yy ; 100x , 100y ; 10x , 10y
xx_quo = INT( x_mod/800) --> 這裡求得圖號 xx (商)
xx_mod = MOD( x_mod,800) --> 剩下的部份， (餘)
x3_100 = INT(xx_mod/100) --> 再除以 100，得百分位分區 x: A~H , 0~7
x3_100n (name)
x3_mod = MOD(xx_mod,100) --> 剩下的部份，
x2_10 = INT(x3_mod/10) --> 再除以 10，得十分位
x1_ = INT(MOD(x3_mod,10))--> 最後剩下的部份就是個位數了
yy_quo = INT( y_mod/500) --> 這裡求得圖號 yy
yy_mod = MOD( y_mod/500) --> 剩下的部份，
y3_100 = INT(yy_mod/100) --> 再除以 100，得百分位分區 y: A~E , 0~4
y3_100n (name)
y3_mod = MOD(yy_mod,100) --> 剩下的部份，
y2_10 = INT(y3_mod/10) --> 再除以 10，得十分位
y1_ = INT(MOD(y3_mod,10))--> 最後剩下的部份就是個位數了
最後將 Area + xx_quo + yy_quo + x3_100n + y3_100n + x2_10 + y2_10 + x1_ + y1_
全部湊起來就是電力座標了。
---------------------------------------------------------------------------
(U)TM 轉 LLh
http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html
http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm
http://www.sunriver.com.tw/grid_tm2.htm
http://www.hydro.nl/pgs/en/pctrans_en.htm PCTrans4.2.5 2.1 3.1
*以下主要採用 uwgb 威斯康辛大學綠灣分校的資料來編寫
*相異處將以 ☆ 表示並說明。
(U)TM: (Universal) Transverse Mercator
LLh: Latitude , Longitude , height (Ellipsoid)
輸入值: 平面座標-in.x in.y (北半球對應為 xE , yN )
輸入值控制參數：
大地基準 datum: (a1, 1/f) or (a1,b1)
計算得： e^2 = 2*f-f^2 = 1-b^2/a^2
, e'^2 = f*(2-f)/(1-f)^2 = a^2/b^2 - 1
*簡化為 e2 , e_2 ; 再加算 e4 = e2^2 , e6 = e2^3
*之後的 a,b,f 沒標示成 a1 b1 f1 為了便於閱讀而省略。
再算得： e1 = (1 - (1 - e2)^(1/2)) / (1 + (1 - e2)^(1/2) )
*注意: e1 不是 e2 開平方根。
接著算： J1 J2 J3 J4
J1 = (3*(e1)/2 - 27*(e1^3)/32) +... *(...) 表示為序列公式
J2 = (21*(e1^2)/16 - 55*(e1^4)/32) +... *本來是無法縮減成閉合公式的。
J3 = (151*(e1^3)/96) +... *因為計算橢圓弧長的公式是為
J4 = (1097*(e1^4)/512) -... *橢圓積分 elliptic integrals
*可因為精度需求只計算到上列所示之公式即可。
TM 橫麥投影控制參數: TM2 TM3 TM6 ?
Central meridian 中央子午線(經線) 121 deg t_rLong0
Central parallel 中央緯線 0 deg t_rLat0
scale factor 比例係數 0.9999 t_k0
False Easting 橫座標平移(東移) 250000m t_FE
False Northing 0m t_FN
計算得：M0 ("t_rLat0" , a1 , e2 , e4 , e6)
M0 = a1*( (1 - (e2)/4 - 3*(e4)/64 - 5*(e6)/256 ...)*t_rLat0
- (3*(e2)/8 + 3*(e4)/32 + 45*(e6)/1024 ...)*SIN(2*t_rLat0)
+ (15*(e4)/256 + 45*(e6)/1024 ...)*SIN(4*t_rLat0)
- (35*(e6)/3072 ...)*SIN(6*t_rLat0) )
*(...) 本為序列公式型態。原因同上。
*USGS use M. Army use S.
#####以下為批次欄位，隨輸入值變化#####
Calculate the Meridional Arc
M1 = M0 + (in.y - t_FN) / t_k0
☆ 採用 POSC 的資料。因為 M0,FN = 0 所以 uwgb 寫 M = y(N) / k0
mu (a1 , e2 , e4 , e6 , M1)
= M1 / ( a1*(1 - (e2)/4 - 3*(e4)/64 - 5*(e6)/256 ) )
Calculate Footprint Latitude
fp,phi_1 (mu, J1,J2,J3,J4)
= mu + J1*SIN(2*mu) + J2*SIN(4*mu) + J3*SIN(6*mu) + J4*SIN(8*mu)
C1 (fp , e_2)
= e_2 * COS(fp)^2
T1 = TAN(fp)^2
R1,rho (a1 , e2 , fp)
= a1*(1-e2) / (1 - e2*SIN(fp)^2)^(3/2)
N1,nu = a1 / (1-(e2)*SIN(fp)^2)^(1/2)
D = (in.x - t_FE) / (N1*t_k0)
☆採用 POSC 的資料。 FE = 500,000 , 記載於該段第一行
儘管只有寫 x/N1*k0 ，加上第一行說明後便是 ( x - 500,000) / N1*k0
Q1 = N1*TAN(fp)/R1
Q2 = ((D^2)/2)
Q3 = (5 + 3*T1 + 10*C1 - 4*C1^2 - 9*e_2)*(D^4)/24
Q4 = (61 + 90*T1 + 298*C1 +45*T1^2 - 252*e_2 - 3*C1^2 )*(D^6)/720
out.rLat = fp - Q1*(Q2 - Q3 + Q4)
Q5 = D
Q6 = (1 + 2*T1 + C1)*(D^3)/6
Q7 = (5 - 2*C1 + 28*T1 - 3*C1^2 + 8*e_2 + 24*T1^2)*(D^5)/120
out.rLong = t_rLong0 + (Q5 - Q6 + Q7)/COS(fp)
最後 out(rLat,rLong) * 180 / PI()
換成度度量 degree.
*個人意見：從平面座標轉換到經緯，無法得到橢球高，h 從缺。
*或考慮另從正高+水準面來獲得。
---------------------------------------------------------------------------
＊Q:南半球了話怎麼辦？
A: lat0 long0 FE FN 設定為負值。
---------------------------------------------------------------------------
LLh 轉 TM
輸入值: 平面座標-in.Lat in.Long -> rLat , rLong
輸入值控制參數：
大地基準 datum: (a2 , 1/f) or (a2 , b2)
再算得: e^2 , e^4 , e^6 ; e'^2 , e'^4 (簡化為: e2,e4,e6 ; e_2 , e_4 )
b = a*(1-f)
TM 橫麥投影控制參數: 例.台灣二度分帶
t2_rLong0 121
t2_rLat0 0
t2_k0 0.9999
t2_FE 250000
t2_FN 0
再算得: M0 ("t2_rLong0" , a , e2 , e4 , e6 )
= a*( (1 - (e2)/4 - 3*(e4)/64 - 5*(e6)/256 ...)*t2_rLat0
- (3*(e2)/8 + 3*(e4)/32 + 45*(e6)/1024 ...)*SIN(2*t2_rLat0)
+ (15*(e4)/256 + 45*(e6)/1024 ...)*SIN(4*t2_rLat0)
- (35*(e6)/3072 ...)*SIN(6*t2_rLat0) )
*以上是 USGS 美國地質勘探局 的用法。
*另一種是 Army (美軍?) 用 S
得再加算 n = (a-b)/(a+b)
S = A'*lat - B'*sin(2*lat) + C'*sin(4*lat) - D'*sin(6*lat) + E'*sin(8*lat)
其中 lat ，算 S0 時用 rlat0 , 算 S1 時為 in.rlat ，同 M 法。
A' = a*(1 - n + (5/4)*(n^2 - n^3) + (81/64)*(n^4 - n^5) ...)
B' = (3*a*n/2)*(1 - n + (7/8)*(n^2 - n^3) + (55/64)*(n^4 - n^5) ...)
C' = (15*a*n^2 /16)*(1 - n + (3/4)*( n^2 - n^3) ...)
D' = (35*a*n^3 /48)*(1 - n + (11/16)*( n^2 - n^3) ...)
E' = (315*a*n^4 / 512)*(1-n) <--單位公尺，已到小數五位，影響甚小
*(...) 表示公式為序列型態。運算時不用
☆ uwgb 的資料應該是打錯字了, "t"an . word&pdf 版 E' /51 應為 /512
另 S 法也記載於 POSC 3.4.4.1a Lambert Conic Near-Conformal
(把 n 乘入了. POSC 的 n 提出來就等於 uwgb 的公式。)
http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34e1.html
#####以下為批次欄位，隨輸入值變化#####
Nu = a / SQRT( 1-(e2)*SIN(in.rLat)^2 )
p = in.rLong - t2_rLong0
M1(USGS)
= a*( (1 - (e2)/4 - 3*(e4)/64 - 5*(e6)/256 ...)*in.rLat
- (3*(e2)/8 + 3*(e4)/32 + 45*(e6)/1024 ...)*SIN(2*in.rLat)
+ (15*(e4)/256 + 45*(e6)/1024 ...)*SIN(4*in.rLat)
- (35*(e6)/3072 ...)*SIN(6*in.rLat) )
dM = M1 - M0 (如果用 S 則是 dS = S1 - S0 )
☆ 參考 POSC 補上 dM . 因為 UTM lat0 = 0 所以 M0 = 0 , and S=M
K1 = dM*t2_k0
K2 = t2_k0*Nu*SIN(in.rLat)*COS(in.rLat)/2 = t2_k0*Nu*SIN(2*in.rLat)/4
K3 = (t2_k0*Nu*SIN(in.rLat)*COS(in.rLat)^3/24)*
(5-TAN(in.rLat)^2 + 9*e_2*COS(in.rLat)^2 + 4*e_4*COS(in.rLat)^4)
☆相對於 POSC 的資料，其實 K3 之後應該還有一個 K3' ...
out.y = (t2_FN) + K1 + K2*p^2 + K3*p^4
☆ uwgb UTM FN = 0 , 故省略。如 FN ≠ 0 於此加入。
K4 = t2_k0*Nu*COS(in.rLat)
K5 = (t2_k0*Nu*COS(in.rLat)^3 / 6)*(1 - TAN(in.rLat)^2 + e_2*COS(in.rLat)^2)
☆相對於 POSC 的資料，其實 K5 之後應該還有一個 K5' ...
tmp.x = K4*p + K5*p^3
out.x = FE + tmp.x
☆根據最後一段敘述，UTM FE = 500,000
---------------------------------------------------------------------------
外傳：關於 POSC 與 UWGB 公式的差異
POSC 石油技術開放標準聯盟。
Petrotechnical Open Standards Consortium
儘管 POSC 一次就把公式全擺在一起列出，POSC 是完整的 Transverse Mercator 公式。
UWGB 可能比較好入門？ 但要注意，UWGB 只轉換 UTM
是預設 FN = 0 , lat0 = 0
平面TM 轉 經緯LL 的部份是完全可對應的。
經過對照後可發現：
out.rlat = fp - Q1*[Q2 - Q3 + Q4]
Q1 = N1*tan(fp)/R1
out.rlong = rlong0 + [Q5 - Q6 + Q7]/cos(fp)
ν1 = N1 , φ1 = fp , ρ1 = R1 , μ1 = mu
fp = mu + (J1)*sin(2*mu)
+ (J2)*SIN(4*mu)
+ (J3)*SIN(6*mu)
+ (J4)*SIN(8*mu) + ...
e1 = (1 - SQRT(1 - e2) ) /
(1 + SQRT(1 - e2) )
之後, μ1 = mu , M1 , T1 皆相同。
C1 = e_2*cos(φ1)^2 *POSC 應該是漏標1 ，應為 φ1 也就是 fp 。
e'^2 二次離心率平方也可以由 e^2 一次離心率平方算出
e_2 = e2 / [1-e2]
D = (in.x - t_FE) / (N1*t_k0)
************************************************************************
LL 轉 TM 則 uwgb 略少一部份。
out.xE = FE + k0*ν[A + (1-T+C)*(A^3)/6 + (5-18*T+T^2+72*C-58*e_2)*(A^5)/120]
K4*p K5*p^3 [---------uwgb 少的部份-----------]
1.
A = (p)*cos(lat) , p = long - long0 , lat=φ,phi ; long=λ,lambda ; ν=Nu
uwgb: K4*p = k0 * Nu * cos(lat) * (p)
'A'
T = tan(lat)^2
C = e_2*cos(lat)^2
2.
k0*ν*(A^3)/6 *(1-T+C) = k0*Nu*(p^3)*cos(lat)^3/6 *(1-T+C)
' A^3 '
uwgb: K5*p^3 = (k0*Nu*COS(Lat)^3/6)*p^3
*(1 - TAN(Lat)^2 + e_2*COS(Lat)^2)
(1 - 'T' + 'C' )
uwgb 若參照 POSC 公式增補 K5'，則為：
K5' = k0*nu*(5 - 18*tan(lat)^2 + tan(lat)^4 + 72*(e_2*cos(lat)^2) - 58*e_2)
*( (p^5)*cos(lat)^5 ) / 120
（乘開了；或補算 T,C,A）
-----------------------------------------------------------------------
K1 K2 p^2 K3 p^4
out.yN = FN + k0*{ M - M0 + ν*tanφ* [ A^2/2 + (5-T+9*C+4*C^2)*A^4/24 +...] }
^^^^
[uwgb 少的部份-----------]
1.
k0*(M - M0)
uwgb: K1 = k0*S = k0*dM
2.
k0* nu*tan(lat)* A^2/2 = k0* nu*tan(lat)* (p^2)*cos(lat)^2 / 2
三角函數：tan(lat) = sin(lat)/cos(lat) , 所以乘 cos(lat)^2 消掉cos平方項
最後等於 sin()*cos()
= k0* nu* sin(lat)*cos(lat) /2 *(p^2)
uwgb: K2*p^2 = k0*Nu*SIN(Lat)*COS(Lat)/2 *(p^2)
2倍角公式：sin(2θ) = 2*sin(θ)*cos(θ)
∴ sin(θ)*cos(θ) = sin(2θ)/2
∴ = k0*Nu*SIN(2*Lat)/4
3.
k0* nu*tan(lat)* (5-T+9*C+4*C^2)* A^4/24
^^^^^^^^^^^^^^^ ^^^^^^
k0*nu*tan(lat)*(p^4)*cos(lat)^4 /24 , 三角函數
= k0*nu* sin(lat)*cos(lat)^3 /24 *(p^4)
uwgb: ↓ ↓ ↓ ↓ ↓
K3*p^4 = (k0*Nu* SIN(Lat)*COS(Lat)^3 /24)*(p^4)
*(5 - TAN(Lat)^2 + 9*e_2*COS(Lat)^2 + 4*e_4*COS(Lat)^4)
'T' '9*C' '4*C^2'
uwgb 若參照 POSC 公式增補 K3'，則為：
K3' = k0*nu*tan(lat)*
(61 - 58*T + T^2 + 600*C - 330*e_2)*A^6 /720
*算到 K3' K5' 時，y差小數 8 位 ， x差小數 6 位 ，實可忽略不計。
----------------END-----------------
