GMTで走向傾斜の記号を描く

最近はGPSロガーなどで位置情報を得ることが容易になってきました.ここではGMTを用いて,緯度経度が分かっている地点の走向傾斜記号を地図上に記入する方法について紹介します.

走向傾斜をGMTの方角表現に合わせる

一般に地層の走向傾斜はN30W60Sや330 60などのように書きます.しかし,GMTでは東が0°,北が90°,西が180°,南が270°であり,そのままではGMTで走向傾斜の記号をうまく描かせることができません.ここでは第一引数に走向を,第二引数に傾斜を与えて,GMTの方角に合わせます.走向傾斜には,N30Eといった四分式または360度式のどちらの形式でも与えることができます.

#!/usr/bin/perl -w
my $strike_in = "N5E";
my $dip_in = "30S";
my ($strike_gmt,$dip_gmt) = &strikedip2gmt($strike_in,$dip_in);
print "$strike_gmt $dip_gmt\n";

#----------------------------------------------------------
# subroutine for transformation of strike-dip to GMT format
# E: 0, N: 90, W: 180, S: 270
sub strikedip2gmt {
  my ($strike_qd, $dip_qd) = ($_[0], $_[1]);

  # strike
  if ( $strike_qd =~ /[NEWS]/ && !$dip_qd =~ /[NEWS]/) {
    die "ERROR: please confirm strike ($strike_qd) and dip ($dip_qd).\n";
  }

  if ($strike_qd =~ /^NS$/ && $dip_qd =~ /E$/) { $strike_qd = "N0W"; $dip_qd =~ s/E/N/; }
  if ($strike_qd =~ /^N\d+(\.\d+)?E$/ && $dip_qd =~ /E$/) { $dip_qd =~ s/E/S/; }
  if ($strike_qd =~ /^N\d+(\.\d+)?E$/ && $dip_qd =~ /W$/) { $dip_qd =~ s/W/N/; }
  if ($strike_qd =~ /^N\d+(\.\d+)?W$/ && $dip_qd =~ /E$/) { $dip_qd =~ s/E/N/; }
  if ($strike_qd =~ /^N\d+(\.\d+)?W$/ && $dip_qd =~ /W$/) { $dip_qd =~ s/W/S/; }
  if ($strike_qd =~ /^NS$/ && $dip_qd =~ /W$/) { $strike_qd = "N0(\.0+)?E"; $dip_qd =~ s/W/N/; }
  if ($strike_qd =~ /^EW$/ && $dip_qd =~ /S$/) { $strike_qd = "N90(\.0+)?E"; }
  if ($strike_qd =~ /^EW$/ && $dip_qd =~ /N$/) { $strike_qd = "N90(\.0+)?W"; }

  # quadrant
  if ( $strike_qd =~ /E$/ && $dip_qd =~ /S$/) {
    $strike_qd =~ s/[NEWS]//g; if ($strike_qd =~ /\d+(\.\d+)?/) { $strike_qd = 90 - $strike_qd; }
  } elsif ( $strike_qd =~ /W$/ && $dip_qd =~ /N$/ ) {
    $strike_qd =~ s/[NEWS]//g; if ($strike_qd =~ /\d+(\.\d+)?/) { $strike_qd = 90 + $strike_qd; }
  } elsif ( $strike_qd =~ /E$/ && $dip_qd =~ /N$/) {
    $strike_qd =~ s/[NEWS]//g; if ($strike_qd =~ /\d+(\.\d+)?/) { $strike_qd = 270 - $strike_qd; }
  } elsif ($strike_qd =~ /W$/ && $dip_qd =~ /S$/) {
    $strike_qd =~ s/[NEWS]//g; if ($strike_qd =~ /\d+(\.\d+)?/) { $strike_qd = 270 + $strike_qd; }

    # azimuth
  } elsif ($strike_qd =~ /^\d+(\.\d+)?$/ && $strike_qd >=0 && $strike_qd < 90 ) {
    $strike_qd = 90 - $strike_qd;
  } elsif ($strike_qd =~ /^\d+(\.\d+)?$/ && $strike_qd >= 90 ) {
    $strike_qd = 450 - $strike_qd;
  } else { print STDERR "Strike ($strike_qd) is not valid.\n"; }

  # dip
  $dip_qd =~ s/[NEWS]$//;

  return ($strike_qd,$dip_qd);
}

走向傾斜の記号を描く

上のサブルーチンを使って,GMTで地図上に走向傾斜の記号を描くためのシェルスクリプトを出力します.複数の走向傾斜データがある場合は,このサブルーチンを繰り返します.

#!/usr/bin/perl -w
my $output="out.eps";
my $id=1;
my $long="144.30";
my $lat="34.30";
my $strikeqd="N60W";
my $dipqd="30S";
my $wayup=1;
my $inverse=0;
my ($gmtline1,$gmtline2,$gmtline3,$gmtline4) = strike_symbol($id,$long,$lat,$strikeqd,$dipqd,$wayup,$inverse);

# output
pprint "$gmtline1 -K -O >> $output\n";
print "$gmtline2 -K -O >> $output\n";
print "$gmtline3 -K -O >> $output\n";
if ($wayup) {  print "$gmtline4 -K -O >> $output\n"; }

#-------------------------------------
# subroutine for strike-dip symbol
sub strike_symbol {
  my ($longitude,$latitude) = ($_[1],$_[2]);
  my ($strike_in,$dip_in) = ($_[3],$_[4]);
  my ($wayup,$inverse) = ($_[5],$_[6]);

  # defaults
  my $strike_bar_length = "0.5";
  my $dip_bar_length = $strike_bar_length*0.3;
  my $bar_width = $strike_bar_length*0.08;
  my $color = "black";
  if ($inverse) { $color="blue"; }
  my $label_font_size = 7;
  my $wayup_symbol_size = "3p";

  # strike symbol
  my ($strike_gmt,$dip_symbol_label) = strikedip2gmt($strike_in,$dip_in);
  my $strike_symbol = "echo $longitude $latitude $strike_gmt $strike_bar_length |  psxy -J -R -Svb$bar_width/0/0 -G$color";
  my $dip_gmt = $strike_gmt - 90;

  # dip symbol
  my $dip_symbol = "echo $longitude $latitude $dip_gmt $dip_bar_length | psxy -J -R -Svt$bar_width/0/0 -G$color";

  # wayup symbol
  my $wayup_symbol;
  my $pi = atan2(1, 1)*4;
  if ($wayup) {
    my $x_shift_wayup = -1*$strike_bar_length/2*cos($strike_gmt*$pi/180); $x_shift_wayup .= "c";
    my $y_shift_wayup = -1*$strike_bar_length/2*sin($strike_gmt*$pi/180); $y_shift_wayup .= "c";
    $wayup_symbol = "echo $longitude $latitude | psxy -J -R -Sc$wayup_symbol_size -G$color -Xa$x_shift_wayup -Ya$y_shift_wayup";
  }

  # dip label
  my $shift = $label_font_size*0.035;	# 1pt ~ 0.035c
  my $x_shift_label = $shift*sin($strike_gmt*$pi/180)-$shift/6; $x_shift_label .= "c";
  my $y_shift_label = $shift*-1*cos($strike_gmt*$pi/180)-$shift/2; $y_shift_label .= "c";
  my $dip_label = "echo $longitude $latitude $label_font_size 0 0 1 $dip_symbol_label| pstext -J -R -G$color -Xa$x_shift_label -Ya$y_shift_label";

  # return
  return ($strike_symbol,$dip_symbol,$dip_label,$wayup_symbol);
}

上のサブルーチンをまとめて実行すると,以下のようなgmtスクリプトが出力されます.


echo 144.30 34.30 330 0.5 |  psxy -J -R -Svb0.04/0/0 -Gblack -K -O >> out.eps
echo 144.30 34.30 240 0.15 | psxy -J -R -Svt0.04/0/0 -Gblack -K -O >> out.eps
echo 144.30 34.30 7 0 0 1 30| pstext -J -R -Gblack -Xa-0.163333333333333c -Ya-0.334676223927187c -K -O >> out.eps
echo 144.30 34.30 | psxy -J -R -Sc3p -Gblack -Xa-0.21650635094611c -Ya0.125c -K -O >> out.eps

面構造を示す記号を描くには,dip_symbolのところを次のようにします.

my $dip_symbol = "echo $longitude $latitude $dip_gmt $dip_bar_length | psxy -J -R -Svt$bar_width/$dip_bar_length/$dip_bar_length -G$color";

出力例です.左はASTER GDEM,中央はSRTM3_v2.1,右は国土地理院の基盤地図情報(数値標高モデル 10mメッシュ:JPGIS2.0(GML)形式)の地形データを用いて描きました.

test4a.png test4c.png test4b.png

度分秒から度に変換する

引数にコロン区切りの度分秒を与えると,度となって返ってきます.

#!/usr/bin/perl
my $lat_dms="44:36:30"; # 度分秒
my $lat_deg=dms2deg($lat_dms); # 度
print "$lat_deg\n";

# Transform from dd:mm:ss to degree
sub dms2deg {
  my $sep = ":";
  my ($dms) = shift;
  my ($deg) = 0;
  my @dms = split(/$sep/, $dms);
  for (my $i = 0; $i <=$#dms; ++$i) { $deg += $dms[$i] / (60**$i); }
  return ($deg);
}

度分秒から度に変換する

引数に度を与えると,コロン区切りの度分秒となって返ってきます.

#!/usr/bin/perl -w
my $lat_deg="44.50";
my $lat_dms=deg2dms($lat_deg);
print "$lat_dms\n";

# Transform from degree to dd:mm:ss
sub deg2dms {
  my $sep = ":";
  my $dd = shift;
  my $deg = int($dd);
  my $min = int(($dd-$deg)*60);
  my $sec = ((($dd-$deg)*60)-$min)*60;
  return "$deg$sep$min$sep$sec";
}

Back