#!/usr/local/bin/perl # #なんでこうしないんだろう？、というようなところは、ただ単に #知らないからだけです。いい方法があれば教えてください。 # #glc-sb.pl -- subroutines for GPS log compress # #更新履歴 #2006/06/25 サブルーチンだけ取り出し # #カシミール３Ｄのpotファイル専用 #書式は以下 #DAT=WGS84/WGS84/0.000/0.00000000/0/0/0 #ORD=POI/tex=/ido=034'41'24'5/kei=135'33'05'7/adr=1/alt=-5/col=12/pda=22-SEP-05/pti=21:44:54 #ORD=POI/tex=/ido=034'41'26'3/kei=135'33'05'4/adr=1/alt=+9/col=12/pda=22-SEP-05/pti=21:45:23 #ORD=POI/tex=/ido=034'41'25'8/kei=135'33'05'3/adr=2/alt=+9/col=12/pda=22-SEP-05/pti=21:47:20 #ORD=POI/tex=/ido=034'41'25'8/kei=135'33'05'3/adr=2/alt=+11/col=12/pda=22-SEP-05/pti=21:47:24 # 緯度 経度 軌跡番号 高度 表示色 日付 時刻 #日付・時刻はUTC #ORD=POI と tex= は固定にしておく #出力は.trkデータ #書式は以下 #H のところは日付を入れる #H SOFTWARE NAME & VERSION #I PCX5 2.09 # #H COORDINATE SYSTEM #U LAT LON DEG # #H R DATUM IDX DA DF DX DY DZ #M G WGS 84 121 +0.000000e+00 +0.000000e+00 +0.000000e+00 +0.000000e+00 +0.000000e+00 # #H TH.POT(1) # #H LATITUDE LONGITUDE DATE TIME ALT #T N34.6901389 E135.5515833 22-SEP-05 21:44:54 -5 #T N34.6906389 E135.5515000 22-SEP-05 21:45:23 9 # #H TH.POT(2) # #H LATITUDE LONGITUDE DATE TIME ALT #T N34.6905000 E135.5514722 22-SEP-05 21:47:20 9 #T N34.6905000 E135.5514722 22-SEP-05 21:47:24 11 $PI = 3.14159265358979; ################################################# #trkデータ一行書き出し #&trk_write($lat,$long,$pda,$pti,$alt) sub trk_write{ # my($lat,$long,$pda,$pti,$alt,$rea) = @_; my($lat,$long,$pda,$pti,$alt) = @_; if($lat < 0){$lat_trk = "S";}else{$lat_trk = "N";} my $lat_trk = $lat_trk . abs($lat); if($long< 0){$long_trk= "W";}else{$long_trk= "E";} my $long_trk= $long_trk. abs($long); # print "T $lat_trk $long_trk $pda $pti $alt $rea\n"; print "T $lat_trk $long_trk $pda $pti $alt\n"; } 1; #return true #trkデータのスレッドのヘッダー書き出し #&trk_head($time) #$time time形式の時間（1970/1/1 00:00:00(UTC)からの秒数） sub trk_head{ my ($time) = @_; my @head_date = localtime($time); my ($year,$mon,$day) = ($head_date[5]+1900,$head_date[4]+1,$head_date[3]); print "\n"; print "H $year\/$mon\/$day\n"; print "\n"; print "H LATITUDE LONGITUDE DATE TIME ALT\n"; } 1; #return true #potデータの書き出し #&print_pot($ido,$kei,$address,$altitude,$col,$pda,$pti) sub print_pot{ my($ido,$kei,$address,$altitude,$col,$pda,$pti) = @_; print "ORD=POI/tex=/ido=$ido/kei=$kei/adr=$address/alt=$altitude/col=$col/pda=$pda/pti=$pti\n"; } 1; #return true #加速度と方向変化を求める #&get_acc_rot($lat0,$long0,$time0,$lat1,$long1,$time1,$lat,2$long2,$time2) #lat long ... degrees #time ... sec(time) sub get_acc_rot{ my ($lat0,$long0,$time0,$lat1,$long1,$time1,$lat2,$long2,$time2) = @_; #移動距離、方位角、東西、南北の移動量 #１点目と２点目 my ($dist0 ,$dir0 ,$dx0 , $dy0) = &hubeny($lat0/180 * $PI, $long0/180 * $PI, $lat1/180 * $PI, $long1/180 * $PI); my $sec0 = $time1 - $time0; #２点目と３点目 my ($dist1 ,$dir1 ,$dx1 , $dy1) = &hubeny($lat1/180 * $PI, $long1/180 * $PI, $lat2/180 * $PI, $long2/180 * $PI); my $sec1 = $time2 - $time1; my ($vx0,$vy0,$vx1,$vy1,$ax,$ay); #速度算出 if ($sec0 >0){ $vx0 = $dx0 / $sec0; $vy0 = $dy0 / $sec0; }else{ $vx0 = 0; $vy0 = 0; } if ($sec1 >0){ $vx1 = $dx1 / $sec1; $vy1 = $dy1 / $sec1; }else{ $vx1 = 0; $vy1 = 0; } my $v0 = sqrt($vx0**2 + $vy0**2); my $v1 = sqrt($vx1**2 + $vy1**2); #加速度算出 if ($sec1 >0){ $ax = ($vx1 - $vx0)/ $sec1; $ay = ($vy1 - $vy0)/ $sec1; }else{ $ax = 0; $ay = 0; } my $acc = sqrt($ax**2 + $ay**2); #回転角度 my $rote = $dir1 - $dir0; if($rote < -$PI){ $rote += 2*$PI; } if($rote > $PI){ $rote -= 2*$PI; } #print "$dist0,$dist1,$v0,$v1,$acc,$rote\n"; return ($acc,$rote,$dist0); } 1; #return true #&get_variables #potデータの１行を与えて経度、緯度、高度、時間などを返す sub get_variables{ my ($ORD , $tex , $ido , $kei , $adr , $alt , $col , $pda , $pti ) = (split /\// ,$_[0]); #緯度経度 split /=/ , $ido; my $ido_nanboku = $_[1]; split /=/ , $kei; my $keido = $_[1]; #スレッドNo.　−　もともと切れていたかつながっていたかを見る split /=/ , $adr; my $address = $_[1]; #高度　−　置き換えて使うのに便利だから。 split /=/ , $alt; my $altitude = $_[1]; #表示色　−　１２は赤。 split /=/ , $col; my $color = $_[1]; #日付、時刻 split /=/ , $pda; my $pdate = $_[1]; split /=/ , $pti; my $ptime = $_[1]; return ($ido_nanboku ,$keido ,$address ,$altitude ,$color ,$pdate ,$ptime); } 1; #return true #&hubeny ($lat1, $long1, $lat2, $long2) # lat long ... radian # # ヒュベニの距離計算式（旧日本測地系） Hubeny # http://www.kashmir3d.com/kash/manual/std_siki.htm #05/05/23 #05/10/04 方位角を追加 #05/10/06 関数名変更 #05/10/08 日付変更線をまたぐ動きに対応 sub hubeny{ my ($lat1, $long1, $lat2, $long2) = @_; my $dP = $lat2 - $lat1; my $dR = $long2- $long1; #日付変更線をまたぐとき if($dR < -$PI){ $dR += 2*$PI; } if($dR > $PI){ $dR -= 2*$PI; } my $P = ($lat1 + $lat2)/2; my $M = 6334834/sqrt( (1 - 0.006674 * sin($P) **2)**3); my $N = 6377397/sqrt( 1 - 0.006674 * sin($P) **2); my $y = $M*$dP; my $x = $N*cos($P)*$dR; my $dist = sqrt($x**2+$y**2); my $alpha = atan2($N*cos($P)*$dR , $M*$dP); my $d_alpha = sin($P) * $dR; my $direction = $alpha - $d_alpha / 2; if ($direction < 0){ $direction += $PI*2; } return ($dist ,$direction , $x , $y); } 1; #return true #&count_days ($yead, $mon, $day) # $year 1-9999 this year = 2001 # $mon 1-12 Jan=1, Feb=2, Mar=3, ... Dec=12 # $day 1-31 # 1/1/1 -> 1 #01/09/12 #01/09/22 バグ修正 sub count_days{ local ($year, $mon, $day) = @_; local (@m_day) = (0,31,59,90,120,151,181,212,243,273,304,334); local ($days) = ($year-1) * 365 + int(($year-1) / 4) - int(($year-1) / 100) + int(($year-1) / 400); $days += $m_day[$mon-1]; $days += $day; #うるう日を過ぎてる場合は１日追加 # if(($mon > 2) && (($year % 400) == 0) || ($year % 4) == 0 && ($year % 100) != 0 ) { $days++;} if((( (!($year % 4)) && ($year % 100)) || !($year %400)) && ($mon > 2)) { $days++;} return $days; } 1; #return true #dmspot2deg ($dmspot) # dmspot ... ddd'mm'ss'1/10s (eg. 034'40'50'5) #pot形式の度分秒を度に変換 sub dmspot2deg{ my ($D , $M , $S ,$dS); ($D , $M , $S , $dS) = (split /\'/,$_[0]) ; my $deg = (($S+$dS/10)/60+$M)/60+$D; return $deg; } 1; #return true #&pdapti2time($pda ,$pti); #pda=19-APR-03/pti=02:17:07 #を1970/1/1 00:00:00からの経過秒数に変換 sub pdapti2time{ my $YEAR_MAX = 70; #00~$YEARMAX -> 20xx ($YEAR_MAX-1)~99 -> 19xx my $DAY_1970_1_1 = &count_days (1970,1,1); my %MONTH = ('JAN',1,'FEB',2,'MAR',3,'APR', 4,'MAY', 5,'JUN', 6, 'JUL',7,'AUG',8,'SEP',9,'OCT',10,'NOV',11,'DEC',12); my ($pda , $pti) = @_; my ($pday , $pmon_name , $pyear) = (split /-/ , $pda); my ($phour , $pmin , $psec) = (split /:/ , $pti); my $pmon = $MONTH{$pmon_name}; if ($pyear > $YEAR_MAX){ $pyear+=1900; }else{ $pyear += 2000; } my $days = &count_days ($pyear , $pmon , $pday) - $DAY_1970_1_1; my $seconds = $days * 24 * 60 * 60 + $phour * 60 * 60 + $pmin * 60 + $psec; return $seconds; } 1; #return true