Taiwan datums/cad2twd67

#!/usr/bin/perl
# cad2twd67 -- # 地籍座標轉換成TWD67
# Taiwan "Cadastral" to TWD67 coordinate transformation
# Copyright       : http://www.fsf.org/copyleft/gpl.html
# Created On      : Fri May 13 23:28:16 2005

# 地籍座標即日據時代沿用至今地籍圖之常見座標。
# 單位為「間」，等於 20/11 米，即六台尺。
# 原點為台中公園的 89 號三角點。

# 只是抄李弘洲﹑朱杏修﹑曾德福 在
# 1999年第十八屆測量學術及應用研討會論文集之公式。
use strict;
use warnings;

# 「六十九年成果」之二度分帶坐標：
my $XCtm69 = 227361.634 + 0.0;
my $YCtm69 = 2632574.582 + 0.0;

# 大約本程式是針對全省。
# 故若只針對單一鄉鎮不遠處，
# 亦已有一已知點，
# 可以試細調以上的 0.0，
# 否則誤差可達幾十公尺。

# 轉換參數
my $XCcad = 5750;
my $YCcad = -21300;
my $A     = 1.8182516286522;
my $B     = -0.004167109289753;

# 每行輸入之地籍坐標（單位為「間」）
my $Xcad;
my $Ycad;

while (<>) {
    die "Enter Xcad, Ycad" if ( ( $Xcad, $Ycad, undef ) = split ) != 3;

    # $Xtmtrn, $Ytmtrn 是由地籍坐標($Xcad, $Ycad)經四參數轉換套合至
    #「六十九年成果」而成之近似二度分帶坐標。
    my $Xtmtrn = $A * ( $Xcad - $XCcad ) - $B * ( $Ycad - $YCcad ) + $XCtm69;
    my $Ytmtrn = $B * ( $Xcad - $XCcad ) + $A * ( $Ycad - $YCcad ) + $YCtm69;
    printf "%.3f\t %.3f\n", $Xtmtrn, $Ytmtrn;
}

# 以上 "69" 大概類似 TWD67 我猜。
# 轉回去？ 要深讀他們的文章...

# Usage example:
# $ more cadfile
# -5000 1250
# -14100 -23500
# $ perl cad2twd67 < cadfile
# 207909.397       2673620.953
# 191260.172       2628657.146