幾個地圖投影測試

 其實是為了製作拼接的大地圖

WGS84 Latitude Longitude: (用等經緯度去繪製)
緯度越高，X越小Y越大，且單位面積也越來越小
當然單位格子不可能是方格
(應該都是梯形)

Transverse Mercator: (橫麥卡托)
平面單位可以是公尺了
需定義中央經線(lon_0) 與從西向東起算(x_0)，定義 scale factor (k)
(ex. 二度 120~122，中央=121；六度 120~126，中央=123)
預期定義的經線邊界(度分帶邊界) 處會有變形(distort)，各別分帶預期無法拼接
低緯度有間隙，高緯度有重疊
每一格是方格，等周等積，但除中央經線以外，每格都是方格北(grid north)
和正北(南北向的經線) 有夾角，離越遠差越大
proj4js: +proj=tmerc

UTM:
平面單位可以是公尺了
經度60等分，每份6度
緯度-80 ~ 84，分21區，從C開始，沒有 I,O(因為避免與1,0混淆)
    每區8度，X區再多4度(72~80~84)
經度判別分區後設定(zone)，緯度從 0 起算，6度分帶
平面XY格子展繪到經緯時得注意分帶邊界
proj4js: +proj=utm

2度、6度分帶計算:
6度: 經度+180度，除每帶6度，除60分區取餘，取整數，+1
2度: 經度+180度，除每帶2度，除180分區取餘，取整數，+1 
UTM緯度分區: 緯度+80度，除每帶8度，除21分區取餘，取整數 (0起算)
再用 "cdefghjklmnpqrstuvwxx" (沒有 i,o) 跑 charAt
+180 +80 看來是抵銷度數 (-180~0~+180 ; -80~0~+80,+84)

分區算回度:
分區*分帶度數-抵銷度數
例: 51區*6度-180度 = 126

Web Mercator: (麥卡托)
平面單位似乎「不是公尺」
算是覆蓋全球，據說北極熊和南極企鵝無法管轄
雖然網格在 google map 上繪製近似方格(X:Y=1:1?)
但緯度越高，方格周長與面積等比例的縮小了
(所以每個格子也是正北了?)
proj4js: +proj=merc

Equal Area Cylindrical: (等積圓柱)
平面單位可以是公尺了
可以覆蓋全球了
等積就是每格面積都一樣了
(看起來也是正北了)
然而這裡等積的代價是，緯度越高，網格越窄長
緯度(lat_ts)可以設定 0, 30, 45 分別是 Lambert, Behrmann, Gall-Peters
Gall-Peters 網格在台灣看起來較扁
Behrmann 即使用到遼寧也還可以接受(X:Y 約 1:1.3)
proj4js: +proj=cea

JavaScript 四捨五入小數位到符合五的倍數

 javascript round:

引用 Math.round 僅能處理整數，且有瑕疵(參見: javascript round bug)

於是實際應用時會包裝如下:

function round(value, decimals) {

  return Number(Math.round(value+'e'+decimals)+'e-'+decimals);

}

這樣小數也能處理了。

另對於四捨五入到某個倍數的問題: (此例為5)

(js round to nearest multiple of ) 

解法是: Math.round(x/5)*5

同樣的，改以包裝過的 round 來處理，則是:

v = 24.056

round(round(v/5,2)*5,2) // 做兩次

就能正確得到 24.05 了

因為第一次是: round(24.056/5,2)*5 = 24.04999...7 (小數位被乘數放大了)

所以再 round 一次約束回小數兩位

round(24.04999...7,2) = 24.05

--

應用: 為何是五的倍數?

用於製作地圖時間隔錯開用。(even or odd)

例如緯經度方格 latitude & longitude grid ( unit 0.1)

選定 longitude 經度於偶數時(ex. 121.2) ，緯度向下平移 0.05  

則在此偶數 long 上，lat 則為 (0.0-> 0.95, 0.1->0.05, 0.2->0.15......) 

coortranscheck

 轉換檢核

應該可基於測試點來做檢查

https://wiki.osgeo.org/index.php?title=Taiwan_datums/Test_points&uselang=zh-tw

而所用的轉換參數，是參考這篇

http://mutolisp.logdown.com/posts/207563-taiwan-geodetic-coordinate-system-conversion

其中 MOI、OSG1、OSG4 (其他不認得) 得注意，分為

A: MOI 與 OSG4 是用於 proj4js ，上河根據結果應該是用OSG4

MOI的東北角是OSG4 約 0.68m,朝向28.4度(北部)

B: PCTrans 與我的則是用 OSG1的參數(Miller Liu, 2008)

而 A的2.8m，朝向30度則是B

在B附近則有台電圖號座標定位系統(台灣EMAP圖資)與成大

分別是我的、成大、台電 在 0.6m內，朝向98度(東南東)

--

不過奇怪的是，proj4js 查看原始碼發現 towgs84 應該是

m,m,m,asec,asec,asec,ppm

其中 asec 讀入後是一律乘以 SEC_TO_RAD ，換言之讀入的是 asec

故 MOI 與 OSG4 不知道為何記載是這樣

過小的 rxyz 幾乎與 0 無異，sf 過小給下去則等同於1

最終 OSG1 若用在 proj4js ，其中 rxyz 負負得正，則可轉換到跟B一樣的位置

應該是 PCTrans/我的 是 PVR，proj4js 是 CFR 的緣故(來源待查)

PVR: Position Vector Rotation
CFR: Coordinate Frame Rotation

參考: 

https://www.bluemarblegeo.com/knowledgebase/calculator-2020sp1/datum_shifts/Seven_Parameter_CFR.htm 

--

原來以前就已經寫過找最近基石點的方式來轉換

將資料納入頁面只需100kb

--

平面轉經緯應該公式沒有問題。

經緯轉平面可能得參考 (U)TM Projection 。

https://github.com/jjimenezshaw/Leaflet.UTM

proj4js 平面轉經緯轉平面，結果是 0

我的和 PCTrans 結果一樣??? 97平面 dx = 0.003m, dy = -0.005m

--

其他:

PCTrans 4.2.8 (自4.2.5版後新增了 Miller Liu 的 TWD67 7參數定義)

PCTrans官方下載頁面 (目前版本5.1)