#!/usr/bin/perl -w
#
# 集成図
# Nowral
# 00/2/28
#
# 初期化
$etime = time;
use GD;
require "GUSI.ph";
$0 =~ /^(.*):/;
$cwd = $1; # カレントディレクトリ
$datadir = "Macintosh HD:地図&GPS:数値地図"; # 画像ディレクトリ
while(! -f "$datadir:KANRI.CSV") {
$datadir = MacPerl::Choose(&GUSI'AF_FILE, 0,
"画像データと管理ファイルはどこにありますか", "", &GUSI'CHOOSE_DIR) || die;
}
$kanricsv = "$datadir:KANRI.CSV";
&GeoTIFF::utminit();
$mpp = 20; # 出来上がり地図の縮尺(ピクセル当り長さ)
# 全体の大きさ
# 張り合わせる経緯度の範囲
# lon0
# lat0 +------+
# | |
# | |
# +------+ lat0
# lon0
#
MacPerl::Ask("緯度の下限、南限? (北緯：入力例 43/10)", "") =~ m#(\d+)[^\d](\d+)# || die;
$elat2 = $1 + $2/60;
MacPerl::Ask("緯度の上限、北限? (北緯：入力例 43/30)", "") =~ m#(\d+)[^\d](\d+)# || die;
$elat0 = $1 + $2/60;
MacPerl::Ask("経度の下限、西限? (東経：入力例 144/50)", "") =~ m#(\d+)[^\d](\d+)# || die;
$elon0 = $1 + $2/60;
MacPerl::Ask("経度の上限、東限? (東経：入力例 145/10)", "") =~ m#(\d+)[^\d](\d+)# || die;
$elon2 = $1 + $2/60;
die unless($elat0>$elat2 && $elon0<$elon2);
# 頭文字 <> エリア
# e -> f -> g -> h
# 四隅番号付け
# 0 - 3
# | |
# 1 - 2
($flat2, $flat0) = ($elat2, $elat0);
for($flonb=int($elon0/6)*6; $flonb<$elon2; $flonb+=6) {
# 図郭四隅経緯度
$flon0 = ($flonb >$elon0) ? $flonb : $elon0;
$flon2 = ($flonb+6<$elon2) ? $flonb+6 : $elon2;
$newargv = sprintf("%s:N%d%02dE%d%02d-N%d%02dE%d%02d.GIF", $cwd,
int($flat2), $flat2*60%60, int($flon0), $flon0*60%60,
int($flat0), $flat0*60%60, int($flon2), $flon2*60%60);
print "# $newargv\n";
print "# 図郭四隅経緯度\n";
printf "左上 : (%.02f, %.02f)\n", $flat0, $flon0;
printf "左下 : (%.02f, %.02f)\n", $flat2, $flon0;
printf "右下 : (%.02f, %.02f)\n", $flat2, $flon2;
printf "右上 : (%.02f, %.02f)\n", $flat0, $flon2;
# UTM座標
# x
# |
# |
# +---- y
$fzn = ($flon0+180) / 6 % 60 + 1; # ゾーン < 西境界含む
print "# ゾーン : $fzn\n";
($fx0, $fy0) = &GeoTIFF::deg2utm($flat0, $flon0, $fzn);
($fx1, $fy1) = &GeoTIFF::deg2utm($flat2, $flon0, $fzn);
($fx2, $fy2) = &GeoTIFF::deg2utm($flat2, $flon2, $fzn);
($fx3, $fy3) = &GeoTIFF::deg2utm($flat0, $flon2, $fzn);
print "# 図郭四隅UTM座標\n";
printf "左上 : (%7.0f, %7.0f)\n", $fx0, $fy0;
printf "左下 : (%7.0f, %7.0f)\n", $fx1, $fy1;
printf "右下 : (%7.0f, %7.0f)\n", $fx2, $fy2;
printf "右上 : (%7.0f, %7.0f)\n", $fx3, $fy3;
# ファイル座標
$fymin = ($fy0<$fy1) ? $fy0 : $fy1;
$fxmin = ($fx1<$fx2) ? $fx1 : $fx2;
$fymax = ($fy2>$fy3) ? $fy2 : $fy3;
$fxmax = ($fx3>$fx0) ? $fx3 : $fx0;
print "# 図郭四隅ファイル座標\n";
printf "左上 : (%7.1f, %7.1f)\n", ($fxmax-$fx0)/$mpp, ($fy0-$fymin)/$mpp;
printf "左下 : (%7.1f, %7.1f)\n", ($fxmax-$fx1)/$mpp, ($fy1-$fymin)/$mpp;
printf "右下 : (%7.1f, %7.1f)\n", ($fxmax-$fx2)/$mpp, ($fy2-$fymin)/$mpp;
printf "右上 : (%7.1f, %7.1f)\n", ($fxmax-$fx3)/$mpp, ($fy3-$fymin)/$mpp;
# イメージのサイズ
$fhght = int(($fxmax-$fxmin) / $mpp) + 1; # イメージの高さ
$fwdth = int(($fymax-$fymin) / $mpp) + 1; # イメージの幅
$fsize = int($fhght * $fwdth / 1024 / 1024 * 2 + 0.5); # 画像メモリ*2
print "# 図郭ピクセル数\n";
printf "縦方向 : %d\n", $fhght;
printf "横方向 : %d\n", $fwdth;
printf "メモリ : %d [MB]\n", $fsize;
MacPerl::Answer("MacPerlの割当メモリは、$fsize MB以上ありますか？",
"大丈夫", "不安かも") || die;
# 必要となる地図の一次メッシュコード
print "# 一次メッシュコード\n";
undef @mc1s;
undef @mc1e;
$i = 0;
for($b1=int($flat2*1.5); $b1<$flat0*1.5; ++$b1) {
for($l1=int($flon0); $l1<$flon2; ++$l1) {
$i++;
$mc1 = sprintf("%02d%02d", $b1, $l1-100); # 一次メッシュコード
if(-f "$datadir:$mc1.TIF") {
push(@mc1s, $mc1);
print "$i : $mc1\n";
}
else {
push(@mc1e, $mc1);
print "$i : $mc1 X\n";
}
}
}
if(@mc1e) { MacPerl::Answer("以下の１次メッシュコードの画像データが見つかりません。かまわず続けますか？
" . join(", ", sort @mc1e), "続行", "中止") || die; }
# 新規イメージ
$im = new GD::Image($fwdth, $fhght);
# カラーパレット < 版の色指定、上位ビットにある版を優先
# 注記版・・・・・・居住地名、自然物等の名称、または、建物記号等
# 墨版・・・・・・・鉄道、行政界
# 赤版・・・・・・・市街地、国道番号
# 褐版・・・・・・・道路、橋及び高架部、建物、堤防、がけ、岩礁等
# 緑版・・・・・・・等高線、、畑記号
# 藍版・・・・・・・水涯線、水田記号、湿地
# 藍マスク版・・・・河川水面、海水面、湖水面
foreach $i (0 .. 255) {
if($i & 0x80) { $im->colorAllocate( 0, 0, 0); } # 注記版
elsif($i & 0x40) { $im->colorAllocate(128,128,128); } # 墨版
elsif($i & 0x20) { $im->colorAllocate(255, 0, 0); } # 赤版
elsif($i & 0x10) { $im->colorAllocate(255,128, 0); } # 褐版
elsif($i & 0x08) { $im->colorAllocate( 0,255, 0); } # 緑版
elsif($i & 0x04) { $im->colorAllocate( 0, 0,255); } # 藍版
elsif($i & 0x01) { $im->colorAllocate(128,128,255); } # 未使用
elsif($i & 0x02) { $im->colorAllocate( 0,255,255); } # 藍マスク版
else { $im->colorAllocate(255,255,255); }
}
# 個々の範囲切り出し
#undef @mc1s;
foreach $mc1 (sort @mc1s) {
$gtime = time;
($glat0, $glon0, $glat2, $glon2) = &GeoTIFF::convinit($kanricsv, $mc1);
# 境界 <> マスク?
$hlat0 = ($flat0 > $glat0) ? $glat0 : $flat0;
$hlon0 = ($flon0 < $glon0) ? $glon0 : $flon0;
$hlat2 = ($flat2 < $glat2) ? $glat2 : $flat2;
$hlon2 = ($flon2 > $glon2) ? $glon2 : $flon2;
($hp0, $hl0) = &GeoTIFF::ll2pix($hlat0, $hlon0);
($hp1, $hl1) = &GeoTIFF::ll2pix($hlat2, $hlon0);
# ($hp2, $hl2) = &GeoTIFF::ll2pix($hlat2, $hlon2);
($hp3, $hl3) = &GeoTIFF::ll2pix($hlat0, $hlon2);
undef %hp;
# 左辺
($hpp, $hlp) = ($hp0, $hl0);
for($b=$hlat0*7200; $b>=$hlat2*7200; --$b) {
($hpn, $hln) = &GeoTIFF::ll2pix($b/7200, $hlon0);
$dhl = $hln - $hlp;
next unless($dhl);
die if($dhl**2>1);
$hp{$hlp} = "$hpp";
($hpp, $hlp) = ($hpn, $hln);
}
$hp{$hln} = $hpn;
# 上辺
($hpp, $hlp) = ($hp0, $hl0);
for($l=$hlon0*7200; $l<=$hlon2*7200; ++$l) {
($hpn, $hln) = &GeoTIFF::ll2pix($hlat0, $l/7200);
$dhl = $hln - $hlp;
if ($dhl== 1) { $hp{$hlp} .= " $hpp"; }
elsif($dhl==-1) { $hp{$hln} .= " $hpn"; }
elsif($dhl**2>1) { die; }
($hpp, $hlp) = ($hpn, $hln);
# print "$hpp , $hlp\n" if($dhl**2);
}
# 下辺
($hpp, $hlp) = ($hp1, $hl1);
for($l=$hlon0*7200; $l<=$hlon2*7200; ++$l) {
($hpn, $hln) = &GeoTIFF::ll2pix($hlat2, $l/7200);
$dhl = $hln - $hlp;
if ($dhl== 1) { $hp{$hln} .= " $hpn"; }
elsif($dhl==-1) { $hp{$hlp} .= " $hpp"; }
elsif($dhl**2>1) { die; }
($hpp, $hlp) = ($hpn, $hln);
# print "$hpp , $hlp\n" if($dhl**2);
}
# 右辺
($hpp, $hlp) = ($hp3, $hl3);
for($b=$hlat0*7200; $b>=$hlat2*7200; --$b) {
($hpn, $hln) = &GeoTIFF::ll2pix($b/7200, $hlon2);
$dhl = $hln - $hlp;
next unless($dhl);
die if($dhl**2>1);
$hp{$hlp} .= " $hpp";
($hpp, $hlp) = ($hpn, $hln);
}
$hp{$hln} .= " $hpn";
# 図葉ファイル
open(IN, "$datadir:$mc1.TIF") || next;
# ファイルヘッダ
read(IN, $_, 2); # バイトオーダー
die unless(unpack("A2", $_) eq "MM");
read(IN, $_, 2); # バージョン番号
read(IN, $_, 4); # 最初のイメージファイルディレクトリの位置
seek(IN, unpack("N", $_), 0);
# イメージファイルディレクトリ（IFD）
read(IN, $_, 2); # タグの数
die unless(unpack("n", $_) == 14);
read(IN, $_, 12); # ピクセル数
($tag, $type, $ndat, $npix) = unpack("nnNn", $_);
die unless($tag==256 && $type==3 && $ndat==1);
read(IN, $_, 12); # ライン数
($tag, $type, $ndat, $nlin) = unpack("nnNn", $_);
die unless($tag==257 && $type==3 && $ndat==1);
read(IN, $_, 12); # 1ピクセル当たりのビット数
read(IN, $_, 12); # 圧縮形式
read(IN, $_, 12); # 色表現法
read(IN, $_, 12); # 画像ブロックの開始位置
read(IN, $_, 12); # カラープレーンの数
read(IN, $_, 12); # 画像ブロックのライン数
read(IN, $_, 12); # 画像ブロックのデータ量
($tag, $type, $ndat, $data) = unpack("nnNN", $_);
die unless($tag==279 && $type==4);
$pos = tell(IN);
seek(IN, $data, 0);
undef @bsiz;
foreach $i (1 .. $ndat) { # 画像ブロックのデータ量（タグ279に対応）
read(IN, $_, 4);
push(@bsiz, unpack("N", $_));
}
seek(IN, $pos, 0);
read(IN, $_, 12); # ピクセル方向解像度
read(IN, $_, 12); # ライン方向解像度
read(IN, $_, 12); # データ保存様式
read(IN, $_, 12); # 解像度の単位
read(IN, $_, 12); # カラーパレット
read(IN, $_, 4); # "0000"
die unless(unpack("N", $_) == 0);
# データ読み込み
print "# 処理済みライン数\n";
seek(IN, 8, 0);
($gp, $gl) = (1, 1);
($gpb, $gpe) = ($npix+1, 0); # unless exists $hp{0}
foreach $i (@bsiz) {
read(IN, $buf, $i);
for($j=0; $j<$i; $j+=2) {
$n = -unpack("c", substr($buf, $j, 1));
$d = unpack("C", substr($buf, $j+1, 1));
if($d==0) { $gp += $n + 1; }
elsif($gp>=$gpb && $gp+$n<=$gpe) { # 範囲内
foreach $k (0 .. $n) {
($hx, $hy) = &GeoTIFF::pix2utm($gp, $gl);
$im->setPixel(int(($hy-$fymin)/$mpp), int(($fxmax-$hx)/$mpp), $d);
$gp++;
}
}
elsif($gp+$n<$gpb) { $gp += $n + 1; } # 範囲手前
elsif($gp>$gpe && !exists $hp{$gl}) { $gp += $n + 1; } # 範囲越え
else { # 境界上
foreach $k (0 .. $n) {
if(exists $hp{$gl} && $gp>$gpe) {
@ed = split(" ", $hp{$gl});
$gpb = shift @ed;
$gpe = shift @ed;
if(@ed) { $hp{$gl} = join(" ", @ed); }
else { delete $hp{$gl}; }
}
if($gp>=$gpb && $gp<=$gpe) {
($hx, $hy) = &GeoTIFF::pix2utm($gp, $gl);
$im->setPixel(int(($hy-$fymin)/$mpp), int(($fxmax-$hx)/$mpp), $d);
}
$gp++;
}
}
next if($gp <= $npix); # 1ライン終り
$gp = 1;
++$gl;
if(exists $hp{$gl}) {
@ed = split(" ", $hp{$gl});
$gpb = shift @ed;
$gpe = shift @ed;
if(@ed) { $hp{$gl} = join(" ", @ed); }
else { delete $hp{$gl}; }
}
else { ($gpb, $gpe) = ($npix+1, 0); }
}
print $gl-1, " / $nlin\n"; # 1画像ブロック終り
}
close(IN);
# 1ファイル終了
$sec = time - $gtime;
printf "# %s time %d:%02d:%02d\n", $mc1, int($sec/3600), $sec/60%60, $sec%60;
}
# 経線
print "# 経線描画...\n";
$dd = 5 / 60; # 5分刻み
for($l=int($flon0/$dd)+1; $l<$flon2/$dd; ++$l) {
for($b=$flat0*7200; $b>=$flat2*7200; --$b) {
($x, $y) = &GeoTIFF::deg2utm($b/7200, $l*$dd, $fzn);
($u, $v) = ( int(($y-$fymin)/$mpp), int(($fxmax-$x)/$mpp) );
$d = $im->getPixel($u, $v);
$im->setPixel($u, $v, $d | 0x40);
}
}
# 緯線
print "# 緯線描画...\n";
$dd = 5 / 60; # 5分刻み
for($b=int($flat2/$dd)+1; $b<$flat0/$dd; ++$b) {
for($l=$flon0*7200; $l<=$flon2*7200; ++$l) {
($x, $y) = &GeoTIFF::deg2utm($b*$dd, $l/7200, $fzn);
($u, $v) = ( int(($y-$fymin)/$mpp), int(($fxmax-$x)/$mpp) );
$d = $im->getPixel($u, $v);
$im->setPixel($u, $v, $d | 0x40);
}
}
# スケール、km方眼
print "# スケール描画...\n";
$dd = 10000 / $mpp; # 10kmおき
# 上辺、下辺
for($ub=0; $ub<$fwdth; $ub+=2*$dd) {
for($u=$ub; $u<$ub+$dd && $u<$fwdth; ++$u) {
foreach $v (0, 1, $fhght-2, $fhght-1) {
$d = $im->getPixel($u, $v);
$im->setPixel($u, $v, $d | 0x40);
}
}
}
# 左辺、右辺
for($vb=0; $vb<$fhght; $vb+=2*$dd) {
for($v=$vb; $v<$vb+$dd && $v<$fhght; ++$v) {
foreach $u (0, 1, $fwdth-2, $fwdth-1) {
$d = $im->getPixel($u, $v);
$im->setPixel($u, $v, $d | 0x40);
}
}
}
# 磁北線
print "# 磁北線描画...\n";
$dd = 5 / 60;
$h = 1 / 7200;
for($li=int($flon0/$dd)+1; $li<=$flon2/$dd; ++$li) {
$l = $li * $dd;
for($i=0; ($b=$flat2+$i*$h)<$flat0 && $l>$flon0; ++$i) {
$d = &GeoTIFF::dldb($b, $l);
$lc = $l + $h * $d;
$b += $h;
$l += $h * ($d + &GeoTIFF::dldb($b, $lc)) / 2;
($x, $y) = &GeoTIFF::deg2utm($b, $l, $fzn);
($u, $v) = ( int(($y-$fymin)/$mpp), int(($fxmax-$x)/$mpp) );
$d = $im->getPixel($u, $v);
$im->setPixel($u, $v, $d | 0x01);
}
}
# GIF出力
print "# GIF出力...\n";
open(OUT, ">$newargv") || die "Couldn't open $newargv :$^E";
MacPerl::SetFileInfo("GKON", "GIFf", $newargv);
print OUT $im->gif;
close(OUT);
undef $im;
# 1ゾーン終了
print "# $newargv\n";
print " ... done\n\n";
}
$sec = time - $etime;
printf "# total time %d:%02d:%02d\n", int($sec/3600), $sec/60%60, $sec%60;
#MacPerl::Quit(2);
die "The Unhappy End";
#!/usr/bin/perl -w
#
# 数値地図(地図画像)キャリブレーション用パッケージ
# Nowral
# 00/2/17
#
package GeoTIFF; # グローバル変数隔離のため
# $pi, $rd, $a, $f, $e2, $ed2, $cad, $cbd, $ccd, $cdd, $ced
# $zn, $x0, $y0, $c, $s, $mpph, $mppv, $p0, $l0
sub dldb { # 偏角
my($b, $l) = @_;
my($bd, $ld) = ($b - 37, $l - 138);
my($d) = 442.82 + 21.01*$bd - 7.36*$ld
- 0.197*$bd*$bd + 0.587*$bd*$ld - 0.961*$ld*$ld;
$d *= -$rd / 60;
return sin($d) / cos($d);
}
sub utminit { # グローバル変数準備
#定数
$pi = 4 * atan2(1, 1); # 円周率
$rd = $pi / 180; # radian / degree
# 日本測地系
$a = 6377397.155; # 赤道半径 (Tokyo)
$f = 1 / 299.152813; # 扁平率 (Tokyo)
$e2 = 2*$f - $f*$f; # 第1離心率
$ed2 = $e2 / (1-$e2); # 第2離心率
# 子午線弧長計算式係数
my($e4, $e6, $e8);
$e4 = $e2 * $e2;
$e6 = $e4 * $e2;
$e8 = $e6 * $e2;
$cad = $a * (1-$e2) * ( 1 + 3/4*$e2 + 45/64*$e4 + 175/256*$e6 + 11025/16384*$e8 );
$cbd = $a * (1-$e2) * ( 3/4*$e2 + 15/16*$e4 + 525/512*$e6 + 2205/2048*$e8 );
$ccd = $a * (1-$e2) * ( 15/64*$e4 + 105/256*$e6 + 2205/4096*$e8 );
$cdd = $a * (1-$e2) * ( 35/512*$e6 + 315/2048*$e8 );
$ced = $a * (1-$e2) * ( 315/16384*$e8 );
return 1;
}
sub convinit { # グローバル変数準備
($kanricsv, $mc1) = @_; # 管理ファイル、一次メッシュコード
# 管理ファイル
$_ = "";
open(IN, $kanricsv) || die;
$_ = <IN> until(/^\r$mc1,/);
close(IN);
chomp;
my(@ed) = split(",");
#１次ﾒｯｼｭｺｰﾄﾞ,図名(漢字),図名(よみ),地勢図ｺｰﾄﾞ,編集,修正,要部修正,発行年月日,図郭四隅緯経度左上緯度,分,秒,左上経度,分,秒,左下緯度,分,秒,左下経度,分,秒,右下緯度,分,秒,右下経度,分,秒,右上緯度,分,秒,右上経度,分,秒,,,,,,,,,図郭四隅ﾌｧｲﾙ座標左上P,左上L,左下P,左下L,右下P,右下L,右上P,右上L
# 図郭四隅経緯度
# 0 - 3
# | |
# 1 - 2
my($lat0, $lon0, $lat1, $lon1, $lat2, $lon2, $lat3, $lon3);
$lat0 = $ed[ 8] + $ed[ 9]/60 + $ed[10]/3600;
$lon0 = $ed[11] + $ed[12]/60 + $ed[13]/3600;
$lat1 = $ed[14] + $ed[15]/60 + $ed[16]/3600;
$lon1 = $ed[17] + $ed[18]/60 + $ed[19]/3600;
$lat2 = $ed[20] + $ed[21]/60 + $ed[22]/3600;
$lon2 = $ed[23] + $ed[24]/60 + $ed[25]/3600;
$lat3 = $ed[26] + $ed[27]/60 + $ed[28]/3600;
$lon3 = $ed[29] + $ed[30]/60 + $ed[31]/3600;
print "# 図郭四隅経緯度\n";
printf "左上 : (%.02f, %.02f)\n", $lat0, $lon0;
printf "左下 : (%.02f, %.02f)\n", $lat1, $lon1;
printf "右下 : (%.02f, %.02f)\n", $lat2, $lon2;
printf "右上 : (%.02f, %.02f)\n", $lat3, $lon3;
die unless($lat0==$lat3 && $lon0==$lon1 && $lat1==$lat2 && $lon2==$lon3);
# 図郭四隅ファイル座標 < 1から勘定
my($p1, $l1, $p2, $l2, $p3, $l3);
($p0, $l0, $p1, $l1, $p2, $l2, $p3, $l3) = @ed[40 .. 47];
print "# 図郭四隅ファイル座標\n";
printf "左上 : (%4d, %4d)\n", $p0, $l0;
printf "左下 : (%4d, %4d)\n", $p1, $l1;
printf "右下 : (%4d, %4d)\n", $p2, $l2;
printf "右上 : (%4d, %4d)\n", $p3, $l3;
# 左上原点
print "# 図郭四隅ファイル座標\n";
printf "左上 : (%4d, %4d)\n", 0, 0;
printf "左下 : (%4d, %4d)\n", $p1 - $p0, $l1 - $l0;
printf "右下 : (%4d, %4d)\n", $p2 - $p0, $l2 - $l0;
printf "右上 : (%4d, %4d)\n", $p3 - $p0, $l3 - $l0;
# UTM座標
# x
# |
# |
# +---- y
$zn = ($lon0+180) / 6 % 60 + 1; # ゾーン
print "# ゾーン : $zn\n";
my($x1, $y1, $x2, $y2, $x3, $y3);
($x0, $y0) = °2utm($lat0, $lon0, $zn);
($x1, $y1) = °2utm($lat1, $lon1, $zn);
($x2, $y2) = °2utm($lat2, $lon2, $zn);
($x3, $y3) = °2utm($lat3, $lon3, $zn);
print "# 図郭四隅UTM座標\n";
printf "左上 : (%7.0f, %7.0f)\n", $x0, $y0;
printf "左下 : (%7.0f, %7.0f)\n", $x1, $y1;
printf "右下 : (%7.0f, %7.0f)\n", $x2, $y2;
printf "右上 : (%7.0f, %7.0f)\n", $x3, $y3;
# 原点移動、軸名付け替え
# yyt
# |
# |
# +---- xxt
my($xxt1, $yyt1, $xxt2, $yyt2, $xxt3, $yyt3);
($xxt1, $yyt1) = ($y1 - $y0, $x1 - $x0);
($xxt2, $yyt2) = ($y2 - $y0, $x2 - $x0);
($xxt3, $yyt3) = ($y3 - $y0, $x3 - $x0);
my($ang) = atan2($yyt3, $xxt3);
print "# 回転角\n";
printf "%.2f\n", $ang / $rd;
$c = cos(-$ang);
$s = sin(-$ang);
# 回転、y軸反転
# +---- xx
# |
# |
# yy
my($xx1, $yy1, $xx2, $yy2, $xx3, $yy3);
($xx1, $yy1) = ($c*$xxt1 - $s*$yyt1, -$s*$xxt1 - $c*$yyt1);
($xx2, $yy2) = ($c*$xxt2 - $s*$yyt2, -$s*$xxt2 - $c*$yyt2);
($xx3, $yy3) = ($c*$xxt3 - $s*$yyt3, -$s*$xxt3 - $c*$yyt3);
# 縮尺(ピクセル当り長さ)
# $mpp = 2e5 * 1e-4; # [m/pixel] (= 20万 x 100 [μm / pixel])
my($c01, $c12, $c23, $c30);
$c01 = $yy1 / ($l1-$l0);
$c12 = ($xx2-$xx1) / ($p2-$p1);
$c23 = ($yy2-$yy3) / ($l2-$l3);
$c30 = $xx3 / ($p3-$p0);
($mpph, $mppv) = (($c12+$c30) / 2, ($c01+$c23) / 2);
print "# 図郭縮尺\n";
printf "縦方向 : %f\n", $mppv;
printf "横方向 : %f\n", $mpph;
print "# 図郭四隅ファイル座標\n"; # 大はずししてないかチェック
printf "左上 : (%7.1f, %7.1f)\n", 0, 0;
printf "左下 : (%7.1f, %7.1f)\n", $xx1/$mpph, $yy1/$mppv;
printf "右下 : (%7.1f, %7.1f)\n", $xx2/$mpph, $yy2/$mppv;
printf "右上 : (%7.1f, %7.1f)\n", $xx3/$mpph, $yy3/$mppv;
return ($lat0, $lon0, $lat2, $lon2);
}
sub ll2pix { # 経緯度 -> ファイル座標
my($lat, $lon) = @_;
my($xx, $yy, $xxt, $yyt, $x, $y, $p, $l);
($x, $y) = °2utm($lat, $lon, $zn);
($xxt, $yyt) = ($y - $y0, $x - $x0);
($xx, $yy) = ($c*$xxt - $s*$yyt, -$s*$xxt - $c*$yyt);
($p, $l) = ($xx/$mpph + $p0, $yy/$mppv + $l0);
return (int($p+0.5), int($l+0.5));
}
sub pix2ll { # ファイル座標 -> 経緯度
my($p, $l) = @_;
my($xx, $yy, $xxt, $yyt, $x, $y, $lat, $lon);
($xx, $yy) = (($p-$p0)*$mpph, ($l-$l0)*$mppv);
($xxt, $yyt) = ($c*$xx - $s*$yy, -$s*$xx - $c*$yy);
($x, $y) = ($yyt + $x0, $xxt + $y0);
($lat, $lon) = &utm2deg($x, $y);
return ($lat, $lon);
}
sub pix2utm { # ファイル座標 -> UTM座標
my($p, $l) = @_;
my($xx, $yy, $xxt, $yyt, $x, $y);
($xx, $yy) = (($p-$p0)*$mpph, ($l-$l0)*$mppv);
($xxt, $yyt) = ($c*$xx - $s*$yy, -$s*$xx - $c*$yy);
($x, $y) = ($yyt + $x0, $xxt + $y0);
return ($x, $y);
}
sub deg2utm { # 経緯度 -> UTM座標
my($b, $l, $zn) =@_;
my($x, $y) = &ll2xy($b*$rd, ($l-$zn*6+183)*$rd);
return ($x * 0.9996, $y * 0.9996 + 500000);
}
sub utm2deg { # UTM座標 -> 経緯度 # 逐次解法のため遅い!
my($x, $y, $zn) = @_;
my($b, $l, $db, $dl, $eps, $cb, $eta2, $rn, $xd, $yd);
$y -= 500000;
$x /= 0.9996;
$y /= 0.9996;
$eps = 0.001 / 3600 * $rd; # 収束条件 < 1000分の1秒
($db, $dl) = (1, 1);
$b = $x / $cad;
$cb = cos($b);
$eta2 = $ed2 * $cb * $cb;
$rn = $a / (1-$f) / sqrt(1+$eta2);
$l = $y / $rn / $cb;
while($db>$eps || $dl>$eps) {
($xd, $yd) = &ll2xy($b, $l);
$db = ($x-$xd) / $cad;
$cb = cos($b);
$eta2 = $ed2 * $cb * $cb;
$rn = $a / (1-$f) / sqrt(1+$eta2);
$dl = ($y-$yd) / $rn / $cb;
$b += $db;
$l += $dl;
}
return ($b/$rd, $l/$rd+$zn*6-183);
}
sub ll2xy { # Gauss-Klueger
my($b, $dl) = @_;
my($dl2, $dl4, $sb, $cb, $cb3, $cb5, $tb2, $tb4, $eta2, $rn, $x, $y);
$dl2 = $dl * $dl;
$dl4 = $dl2 * $dl2;
$sb = sin($b);
$cb = cos($b);
$cb3 = $cb * $cb * $cb;
$cb5 = $cb3 * $cb * $cb;
$tb2 = $sb * $sb / $cb / $cb;
$tb4 = $tb2 * $tb2;
$eta2 = $ed2 * $cb * $cb;
$rn = $a / (1-$f) / sqrt(1+$eta2); # 卯酉線曲率半径
$x = $cad*$b - $cbd/2*sin(2*$b) + $ccd/4*sin(4*$b)
- $cdd/6*sin(6*$b) + $ced/8*sin(8*$b); # 赤道から緯度Bまでの子午線弧長
$x += $rn/2*$sb*$cb*$dl2;
$x += $rn/24*$sb*$cb3*(5-$tb2+9*$eta2+4*$eta2*$eta2)*$dl4;
$x += $rn/720*$sb*$cb5*(61-58*$tb2+$tb4)*$dl4*$dl2;
$y = $rn*$cb*$dl;
$y += $rn/6*$cb3*(1-$tb2+$eta2)*$dl*$dl2;
$y += $rn/120*$cb5*(5-18*$tb2+$tb4)*$dl*$dl4;
return ($x, $y);
}