#!/usr/bin/perl
# longlat2taiwanmapno -- 將經緯度換成台灣圖號
# Convert longitude/latitude to Taiwan map number
# Copyright       : http://www.fsf.org/copyleft/gpl.html
# Author          : Dan Jacobson -- http://jidanni.org/geo/taiwan_datums/
# Created On      : Thu Jan 17 01:28:05 2002
# Last Modified On: Tue Apr 11 22:16:50 2006
# Update Count    : 242
#Reports Taiwan government TWD67 based map numbers. Usage:
#$ echo 120.85788 24.18347|./longlat2taiwanmapno
##1> 120.85788 24.18347
#9521-2-025 9521-2-09 9521-2-NW
#Or $ ./longlat2taiwanmapno file... #one pair per line
#If point falls on a map boundary, this program reports the north
#and or east map number.
use strict;
use warnings;
while (<>) {
    die "$0: garbled input at line $.\n" unless /^([\d.]+)\s+([\d.]+)$/;
    my $x    = $1;
    my $y    = $2;
    my $n100 = ( $x * 40 ) % 10 + 91 - ( $y * 40 ) % 10 * 10;
    my $n25  = ( $x * 10 ) % 5 + 1 + 5 * ( 4 - ( $y * 20 ) % 5 );
    my $nesw = ( ( $y * 8 ) % 2 ? "N" : "S" ) . ( ( $x * 8 ) % 2 ? "E" : "W" );
    my $xtmp = ( $x * 4 ) % 2;
    my $n4   = ( $y * 4 ) % 2 ? ( $xtmp ? 1 : 4 ) : ( $xtmp ? 2 : 3 );
    my $nx2 = ( ( $x - 73 ) * 2 ) % 100;    #73: 新疆, %:釣魚台
    my $ny2 = int( ( $y - 13.5 ) * 2 );     #南海
    print "#$.> $_";
    printf "%02d%02d-%d-%03d ", $nx2, $ny2, $n4, $n100;    #.025 1:5000:
    printf "%02d%02d-%d-%02d ", $nx2, $ny2, $n4, $n25;     #.05  1:10000:
    printf "%02d%02d-%d-%s\n",   $nx2, $ny2, $n4, $nesw;    #.125 1:25000:
#   printf "%02d%02d-%d\n", $nx2, $ny2, $n4;                #.25  1:50000:
}

#當然不見得所列的圖都有出版，另有外島併入他圖
#$n100:政府設計：本來00至99就可以，豈至100，怕與25衝...
#$n25: the government should also have made it 00 to 05, 00 to 50, total 25.