産総研：研究者

序言

白嶺丸で得られた重磁力データ

データおよびデータ収録

調査船白嶺丸で得られたデータは、船上で処理され、1航海1ファイルに整理される。
航跡図はこちら。
このファイルはデータ交換によく用いられるMGD77フォーマットへの変換が容易なCR80フォーマット（岸本ほか、1984）に揃えられている。
具体的なフォーマットはこちらに載せた。
CR80フォーマットで残されている1974年から1996年までの日本周辺海域の40航海分のデータのうち、すべてに磁気異常が含まれているわけではない。
1970年代は広域調査が中心で、これらのデータの多くは100万分の1の地質図の作成に用いられている。そのため、測線間隔も粗く、また位置の精度も悪いので、今回の作業からはずし、 20万分の1海洋地質図用にデータが取得されたもののみを対象とした。
また、1997年以降の航海ではデータ取得システムに変更がありフォーマットがすこし違っていることと、まだ出版準備中のものも多いためここには含まれていない。

2008年までに出版されている海洋地質図を示した。

また、海洋地質図出版後に当該海域で取得されたデータがある場合は括弧のなかに示した。

作業当時の未解決点

白嶺丸の測位システムは、1990年にGPSが導入されたが、受信状態が安定したのは1996年ころである（上嶋・駒澤，1996）。そのため、90年代前半はまだ未成熟で、LORAN-C、NNSS、GPSによるハイブリッドシステムとなっていた。それ以前のデータは、LORAN-CとNNSSのハイブリッド、さらにはNNSSのみ、など精度が低くなっている。
そのため、測線ごとに磁気異常データを平滑化した上で、コンター図を描いても、交点誤差が大きく実用的なメッシュデータの作成は難しかった。
よって、磁気異常図（海洋地質図の付図）を作成した際には、空中磁気のデータ処理に用いられている交点コントロール手法を適用して、測線ごとのバイアス値を計算してやることによって交点誤差を小さくした。しかし、どうしても交点誤差が大きく残ってしまう部分もあった。
その原因として、測定精度が不十分であったことが挙げられるが、その他に、それぞれの時点で設定されていたIGRFモデルが地磁気の永年変化を正しく表現していない可能性もある（中塚、2001）。
そこでDGRF残差を計算する。さらに、位置を日本測地系から世界測地系へ移行したことにより、新しい緯度経度で、DGRFを計算する必要が生じた。
さらに、データ取得時には、日変化補正のための定点観測を行っていなかったので、海洋地質図出版に際しても日変化補正は行わなかった。夜間データが主で、このデータセットにおいて、日変化は交点誤差に比べて小さい、と判断された（森尻・山崎、1994）ことによる。今回も日変化補正は行わない。

データ処理

CR80フォーマットのファイルから、年、日、分、緯度、経度、船速、船首方位、補正水深、全磁力、IGRF残差、重力計の読み、絶対重力、フリーエアー異常、を抜き出す。
磁気異常測定中は、船はほぼ10ktで航走した。これは1分間に約300mというスピードになる。ちょうどプロトン磁力計のセンサーは船尾から後方約250mのところで曳航されているので、全磁力のデータを1分後のものにずらす
海洋情報部が公開しているプログラム（大門ほか、2000）を用いて緯度、経度を世界測地系に変換
新しい緯度経度でエトベス補正、正規重力値を再計算し、フリーエアー異常を再計算する。
DGRFを計算し、磁気異常をDGRF残差に入れ替える。DGRFの計算は係数IGRF10(Macmillan and Maus, 2005)を用いた。
ここから磁気異常について、主測線、交差測線のファイルに分け、メッシュデータを作成する。

データファイルの作成

STEP1

入力データ：GH90A.DATA
出力データ：gh90a.data0

プログラム:cr80selmag.f

Headerを削除し
緯度経度が入っていない部分（例：赤字のライン）を取り除き
磁気異常のデータが無い部分を取り除き
磁気異常の絶対値が2500nTを超えるものはプロトン磁力計の不調の可能性が大きいので取り除く
さらに数値はみな整数で入っているので正しい値にする。

GH90A.DATA

11GH90-A  90103FUNABASHI, JAPAN            90121NIIGATA, JAPAN              NBMG
12 IGRF85TOKYO GMT                      IAG1967  IGSN71  901039797894 42565   09
21  41313,1314,1414,1413,9999,9999,9999,9999,9999,9999,9999,9999,9999,9999,9999,
22  49999,9999,9999,9999,9999,9999,9999,9999,9999,9999,9999,9999,9999,9999,9999,
390103 240N99999999999999999V99999999B999999999999M9999999999G999999999999999999
390103 241N99999999999999999V99999999B999999999999M9999999999G999999999999999999

390106 517N 3711521 13789532V10091982B  705  69178M47381  -21G 436079798749 -390
390106 518N 3711204 13789581V10591984B  697  68378M47303  -97G 436199798750 -388
390106 519N 3710905 13789609V10221981B  643  63078M47294 -105G 436429798782 -353
390106 520N 3710617 13789511V 9702119B  612  60078M47321  -76G 436719798695 -438

gh90a.data0

90 106  517  37.11521  137.89532 10.09 198.2   691 47381   -21  4360.7  979874.9   -39.0
90 106  518  37.11204  137.89581 10.59 198.4   683 47303   -97  4361.9  979875.0   -38.8
90 106  519  37.10905  137.89609 10.22 198.1   630 47294  -105  4364.2  979878.2   -35.3
90 106  520  37.10617  137.89511  9.70 211.9   600 47321   -76  4367.1  979869.5   -43.8

iyr day min Lat Lon speed head depth mag-total mag-anomaly gravity-reading gravity-total free-air gravity

末尾に00000という行が入ってしまうので手動で削除しておく。

STEP2：測線ごとにプロトン磁力計のセンサーの位置補正を行う。

10kt前後で航行しているときは、位置を基準にして、1秒後の全磁力が正しい組み合わせ。
headingの変化の絶対値が15度以内のときは1秒後の全磁力値を採用する。
5分間隔のデータには適用しない。
入力データ：gh90a.data0
出力データ：gh90a.datam0
プログラム:cable.f

データを読み込み
連続した時刻のところは磁気の値を1分後のものに入れ替える。

gh90a.datam0

90 106  517  37.11521  137.89532 10.09 198.2   691 47303   -97  4360.7  979874.9   -39.0
90 106  518  37.11204  137.89581 10.59 198.4   683 47294  -105  4361.9  979875.0   -38.8
90 106  519  37.10905  137.89609 10.22 198.1   630 47321   -76  4364.2  979878.2   -35.3
90 106  520  37.10617  137.89511  9.70 211.9   600 47303   -94  4367.1  979869.5   -43.8

iyr day min Lat Lon speed head depth mag-total mag-anomaly gravity-reading gravity-total free-air gravity

STEP3:変針中の値を除く→データ末尾にflgを立てる。

iflg
データは通し番号を振る。変針中は番号が０．
1分間に10度以上headが変化した場合は変針中とみなす。iflg=1
時間の切れ目は測定の切れ目：時間が空いた初めのデータ＝2．ライン上のデータ＝0
入力データ：gh90a.datam0
出力データ：gh90a.datas01

プログラム：forlatlon0.f

gh90a.datas01

      1 90 106  517  37.11521 137.89532 10.09 198.2   691 47303   -97  4360.7  979874.9   -39.0 2
      2 90 106  518  37.11204 137.89581 10.59 198.4   683 47294  -105  4361.9  979875.0   -38.8 0
      3 90 106  519  37.10905 137.89609 10.22 198.1   630 47321   -76  4364.2  979878.2   -35.3 0
      0 90 106  520  37.10617 137.89511  9.70 211.9   600 47303   -94  4367.1  979869.5   -43.8 1
      4 90 106  521  37.10419 137.89276  9.74 206.7   533 47275  -122  4369.9  979876.7   -36.3 0
      5 90 106  522  37.10229 137.88928  9.80 208.3   484 47274  -122  4371.9  979877.1   -35.6 0

STEP4:変針中の値を除く

iflg
変針中のデータは除く。
入力データ：gh90a.datas01
出力データ：gh90a.datas1

プログラム：forlatlon1.f

gh90a.datas1

      1 90 106  517  37.11521 137.89532 10.09 198.2   691 47303   -97  4360.7  979874.9   -39.0 2
      2 90 106  518  37.11204 137.89581 10.59 198.4   683 47294  -105  4361.9  979875.0   -38.8 0
      3 90 106  519  37.10905 137.89609 10.22 198.1   630 47321   -76  4364.2  979878.2   -35.3 0
      4 90 106  521  37.10419 137.89276  9.74 206.7   533 47275  -122  4369.9  979876.7   -36.3 2
      5 90 106  522  37.10229 137.88928  9.80 208.3   484 47274  -122  4371.9  979877.1   -35.6 0
      6 90 106  523  37.09972 137.88605 10.76 209.0   436 47256  -140  4372.9  979874.7   -37.8 0

STEP5：測地系を変換する（東京→世界）。水路部のプログラムを使うための前準備。
- プログラムが変ってからは試していない。以下は海洋情報部のサイトから
  - 一度に複数の点を変換する場合（とはいえ、そんなに大量のデータは無理な感じ）
    - エリアに順に緯度の度分秒、経度の度分秒を入力して下さい。
    - スペース区切りで一行に六つの数字を書いて下さい。
      例) 35 45 27.5 141 16 23.5
      36 00 00.0 135 00 00.0
      　　　・
      　　　・
      と入力して下さい。
      度および分の欄にも小数を入力可能になりました。
      <注意>
      分の欄に小数を入力する場合、秒の欄には0を入力して下さい。
      度の欄に小数を入力する場合、分および秒の欄には0を入力して下さい。
      0は省略しないで下さい。一行に六つの欄が必要です。
      重複して入力すると換算のうえ合算されます。
      下記の三例は全て同じ結果になります。
      35 45 30.0 141 15 24.0
      35 45.50 0 141 15.40 0
      35 45.25 15.0 141 14.75 39.0
  - 緯度経度のファイルを作る。
    書式にしたがって、緯度経度のみのファイルを作成する。
  - 入力データ：gh90a.datas1
  - 出力データ：gh90a.latlon
  - プログラム：fordatum3.f
    gh90a.latlon
```
 37  6 54.761 137 53 43.169
 37  6 43.349 137 53 44.927
 37  6 32.583 137 53 45.916
 37  6 22.214 137 53 42.400
```

STEP6：海洋情報部提供の測地系変換プログラムを使う

入力データ：gh90a.latlon　ウェブサイトの窓にコピー＆ペースト
出力データ：別窓に現れるデータをコピー＆ペーストしてテキストファイルを作成する　gh90a.wgs.latlon

結果：gh90a.wgs.latlon(2023.08.15)

「日本測地系 → 世界測地系(WGS84)」,「歪み補正を行う」 により変換しました。
  Version : 3.1.4H-I

N37ﾟ  6' 54.761" E137ﾟ 53' 43.169"  =>  N37ﾟ  7'  5.649" E137ﾟ 53' 31.991"  海上の点です。海洋情報部の式によって変換しました
N37ﾟ  6' 43.349" E137ﾟ 53' 44.927"  =>  N37ﾟ  6' 54.238" E137ﾟ 53' 33.749"  海上の点です。海洋情報部の式によって変換しました
N37ﾟ  6' 32.583" E137ﾟ 53' 45.916"  =>  N37ﾟ  6' 43.474" E137ﾟ 53' 34.739"  海上の点です。海洋情報部の式によって変換しました
N37ﾟ  6' 22.214" E137ﾟ 53' 42.400"  =>  N37ﾟ  6' 33.106" E137ﾟ 53' 31.223"  海上の点です。海洋情報部の式によって変換しました

「日本測地系 → 世界測地系(WGS84)」,「歪み補正を行う」により変換しました。
Version : 3.1.4H-I

この部分を手動で削除する。

STEP7：このデータをほしい数字にする。

入力データ：gh90a.wgs.latlon←headerを手動で削除すみ。
出力データ：gh90a-wgs.txt

datumからデータを出力予定フォルダへコピーしてくる。
gh90a-wgs.txt
```
   37.11824   137.89223
   37.11507   137.89272
   37.11208   137.89299
```

STEP8：位置情報・フリーエアー異常を入れ替える

入力データ：gh90a.datas1, gh90a-wgs.txt
出力データ：gh90a.bs1
プログラム：gnorm-gh.f

ここで行っているのは
1. 新測地系による緯度を使って正規重力を再計算
2. 旧測地系の緯度からCR80で与えたエトベス補正値を計算
3. 新測地系の緯度からエトベス補正値を計算
4. フリーエアー異常を再計算
5. データを入れ替える。
  - CR80ファイルでgravity anomalyとなっているのはフリーエアー異常(ga)である。
    具体的には、以下の式で計算されていた。
    etv=(7.5*cos(Lat)*sin(head)+0.0042*speed)*speed
    ga=(gt+etv)-gnorm+0.87
    測地基準系1967に基づく正規重力式
    gnorm = 978031.85*(1+0.005278895*sin²(Lat) +0.000023462*sin⁴(Lat))
  - 測地基準系1980に基づく正規重力式を使う。
    gnorm = 978032.67715*(1+0.0052790414*sin²(Lat) +0.0000232718*sin⁴(Lat))
    gnorm = 978032.67715 + 5163.07*sin²(Lat) +22.76*sin⁴(Lat)
  - エトベス補正も緯度が変わるので計算しなおす。
    物体が赤道に沿って地球に対して速度Vで東向きに移動すると自転速度ωR（465m/sec）から（ωR+V）に増える。
    したがって遠心力の差は
    2ωV＝((ωR+V)²/R)-((ωR)²/R)
    ωは地球の自転角速度、Vは船速度、Latは緯度、αが方位角で、エトベス補正値は
    δgE=2ωVcos(Lat)sin(α)+(V²/R)
    GRS80(WGS84)では、地球の半径R = 6378137.000ｍ、扁平率F = 1.0d0/298.257223563d0
    速度Vが１ノット（1.8km/h=0.5m/sec）とすると、
    2ωV＝2*(0.5m/sec)*((465m/sec)/6378137m)= 2*0.5*0.000729 m/sec2= 0.00729 gal= 7.29mGal
    となる。
  - 具体的には、Vの単位はノットなので、エトベス補正値は
    etvn=（7.3*cos(Lat)*sin(head)+ 0.00392*V）*V
  - 絶対重力値(wgo)＝元の値(wgobs)―元のエトベス(etv)＋新しいエトベス(etvn)
  - よって、新しいフリーエアー異常は
    ganew=wgo-wgnorm+0.87
  - gh90a.bs1
```
      1 90 106  517  37.11824 137.89223 10.09 198.2   691 47303   -97  4360.7  979875.4   -39.7 2
      2 90 106  518  37.11507 137.89272 10.59 198.4   683 47294  -105  4361.9  979875.5   -39.3 0
      3 90 106  519  37.11208 137.89299 10.22 198.1   630 47321   -76  4364.2  979878.7   -35.8 0
      4 90 106  521  37.10722 137.88966  9.74 206.7   533 47275  -122  4369.9  979877.4   -36.7 2
      5 90 106  522  37.10532 137.88618  9.80 208.3   484 47274  -122  4371.9  979877.8   -36.1 0
      6 90 106  523  37.10275 137.88295 10.76 209.0   436 47256  -140  4372.9  979875.5   -38.2 0
```
    最終行に変な値が入っていたら削除しておく。

STEP9：DGRF残差を計算する

入力データ：gh90a.bs1
出力データ：gh90a.bs2

プログラム：igrf10-02.f
igrf10coeffs.dataが必要

ここで行っているのは

新測地系による位置を使ってDGRFを計算

磁気異常を再計算

データを入れ替える。

IGRF-10のプログラムはここから取った（↓）
http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
これをすこし書き直した。

gh90a.bs2

     1 90 106  517  37.11824  137.89223  10.09 198.2   691  47303.0   -94.8 4360.7  979875.4   -39.7 2
     2 90 106  518  37.11507  137.89272  10.59 198.4   683  47294.0  -101.9 4361.9  979875.5   -39.3 0
     3 90 106  519  37.11208  137.89299  10.22 198.1   630  47321.0   -73.1 4364.2  979878.7   -35.8 0
     4 90 106  521  37.10722  137.88966   9.74 206.7   533  47275.0  -117.4 4369.9  979877.4   -36.7 2
     5 90 106  522  37.10532  137.88618   9.80 208.3   484  47274.0  -118.4 4371.9  979877.8   -36.1 0

データをチェックして、この形式のファイルを保存しておくと便利

メッシュデータを作れなかった言い訳

諸般の事情がありまして、中断しました
石原丈実（2021）GSJ研究資料集No714のデータをご利用ください。広域にわたってレベル補正済みです。
解説はこちら　⇒　石原丈実（2021）白嶺丸重磁力データの整備・公開．地質調査研究報告 Vol.72 No.5, 421-445
gh**.mag.dataを主（横）測線のheadingが280-300°と100-120°、交差（縦）測線のheadingが20-40°と200-220°として横測線と縦測線にawkを使って分ける。
そのあとエクセルで番号を振る。縦は1000番から。
極端なスパイクノイズは手作業ではずす

交点を抜き出し誤差をチェックする→ プログラムはこちら

交差測線(tlat,tlon)に対して主測線(yokoa,yokob)の緯度経度が、以下の条件以下のものを交点とみなす
- dlat=abs(tlat-yokoa(kk))
- dlon=abs(tlon-yokob(kk))
- if(dlat.le.0.008 .and. dlon.le.0.008) then
- dd=sqrt(dlat*dlat+dlon*dlon)
- if(dd.le.0.005) then：30秒が約１ｋｍ

gh90-res.txtを作成。
横側線、縦側線、縦緯度、縦経度、縦磁気異常、横緯度、横経度、横磁気異常

   11 1004 37.99943 138.16904  -140.7 37.99806 138.16782  -129.2
   11 1004 37.99943 138.16904  -140.7 37.99913 138.16435  -144.8
   11 1004 38.00234 138.16997  -143.1 37.99806 138.16782  -129.2
   30 1006 38.11720 138.72757  -165.2 38.12177 138.72604  -152.0
   30 1006 38.11720 138.72757  -165.2 38.12083 138.72919  -162.6
   21 1007 38.00813 138.71748  -252.3 38.00684 138.72217  -243.3
   21 1007 38.00813 138.71748  -252.3 38.00772 138.71906  -253.7
   21 1007 38.00813 138.71748  -252.3 38.00861 138.71594  -282.1
   21 1007 38.00813 138.71748  -252.3 38.00957 138.71288  -312.5
    4 1007 37.33921 138.20020   135.4 37.34419 138.19974   104.2
    4 1007 37.33921 138.20020   135.4 37.34314 138.20319   104.8
    5 1007 37.29670 138.17868   127.9 37.29358 138.17926   121.4
    5 1007 37.29670 138.17868   127.9 37.29435 138.17615   119.1
    5 1007 37.29382 138.17700   117.0 37.29358 138.17926   121.4
    5 1007 37.29382 138.17700   117.0 37.29435 138.17615   119.1
    5 1007 37.29382 138.17700   117.0 37.29507 138.17291   123.7
    5 1007 37.29102 138.17523   120.1 37.29358 138.17926   121.4
    5 1007 37.29102 138.17523   120.1 37.29435 138.17615   119.1
    5 1007 37.29102 138.17523   120.1 37.29507 138.17291   123.7

これを整理する。
必ず2測線以上と交わる必要があるので、これで言うと1004と11は削除。
エクセルで読み込んでとりあえず手作業。
5と1007のような複数の点で同じ線同士が接近する例では縦のデータも単純平均を取り、横の磁気値も単純平均を取る。これを交点磁力値差とする。

横側線、縦側線、縦緯度平均、縦経度平均、縦磁気異常平均、横磁気異常平均
21	1007	138.7175	38.00813	-252.3	-272.9
4	1007	138.2002	37.33921	135.4	104.5
5	1007	138.1787	37.2967	127.9	120.25
5	1007	138.177	37.29382	117	121.4
5	1007	138.1752	37.29102	120.1	121.4
15	1007	138.1494	37.24348	2	18.95
15	1007	138.1482	37.24052	2.3	18.95
9	1009	137.947	38.2365	92.8	73.93333
10	1009	137.9474	38.1197	52.3	72.65
28	1009	137.9516	37.88848	-30.65	-25.05
4	1010	137.9459	37.42342	34.3	21.4
46	1010	138.6865	38.65596	-9.6	23.5
9	1012	137.9044	38.25132	79.35	63.55
8	1012	137.93	38.32343	59.48889	51.72222
43	1012	137.9329	38.33104	52.35	40.475
7	1012	137.944	38.36149	11.23333	-5.8
6	1016	137.9618	38.41098	-40.6	-62.8333
45	1012	137.994	38.52003	-113.28	-117.1
47	1013	138.1876	38.98818	-33.8	-9.6
46	1013	138.0607	38.84655	-11.6125	7.55
45	1013	138.0201	38.80275	-48.5	-39.0667
30	1014	138.5235	38.18173	-88.34	-78.06
21	1014	138.8838	38.77664	-91.3	-108.767
6	1016	138.0112	38.39714	-57.3	-56.35
33	1018	137.9243	38.59323	-90.4	-86.7333
43	1018	137.9564	38.62909	-45.16	-50.16
44	1018	137.9667	38.64054	-39.4	-40.45
48	1018	138.0073	38.6858	-33.7	-32.8
24	1018	138.0431	38.72579	20.5	19.26667
47	1018	138.0942	38.7828	-58.9	-59.2714
45	1018	138.0349	38.80332	-79.5	-66.6667

横―縦を磁気値差とする。

通し番号を振る。縦、横、元縦、元横、横―縦

1	1	1007	4	-30.9
3	1	1010	4	-12.9
1	2	1007	5	-0.65
7	3	1016	6	-10.6417
4	4	1012	7	-17.0333
4	5	1012	8	-7.76667
2	6	1009	9	-18.8667
4	6	1012	9	-15.8
2	7	1009	10	20.35
1	8	1007	15	16.65
1	9	1007	21	-20.6
6	9	1014	21	-17.4667
8	10	1018	24	-1.23333
2	11	1009	28	5.6
6	12	1014	30	10.28
8	13	1018	33	3.666667
4	14	1012	43	-11.875
8	14	1018	43	-5
8	15	1018	44	-1.05
4	16	1012	45	-3.82
5	16	1013	45	9.433333
8	16	1018	45	12.83333
3	17	1010	46	33.1
5	17	1013	46	19.1625
5	18	1013	47	24.2
8	18	1018	47	-0.37143
8	19	1018	48	0.9

ここからプログラムclc.fの入力データを作る。（reform.f）→ プログラムはこちら
最小二乗法で測線同士のバイアスを計算する→ プログラムはこちら

出力　memoc.txt

これは縦が先に来るから、gh90は縦測線が8本だったので、１－８までは縦の番号、それ以降は横の番号に手入力で戻す。
```
1007	5.8
1009	-12.6
1010	23.8
1012	-9.5
1013	9.8
1014	8.9
1016	-13.4
1018	-1.1
4	36.7
5	6.4
```
このバイアス値を全てのデータに足す。→ プログラムはこちら

図化を行う

測線図

gh90.bs2.dataより
awk '{print $6, $5}' gh90.bs2.data > gh90.pos
として以下のファイルを作る。
```
   137.89223 37.11824
   137.89272 37.11507
   137.89299 37.11208
  
```

GMTを起動して500mメッシュの地形データと重ねがきをする。

   grdcontour bathy500m.grd -R135:00/142:00/37:00/44:00 -JM10 -C500 -P -K -V > D:\Myhome\yoreyore\gh90.eps
   psxy D:\GHdata2008\gh90.pos -R135:00/142:00/37:00/44:00 -JM10 -Sc0.01 -G255/0/0 -K -O -V >> D:\Myhome\yoreyore\gh90.eps
   pscoast -R135:00/142:00/37:00/44:00 -JM10 -Df -W2 -K -O -V >> D:\Myhome\yoreyore\gh90.eps
   psbasemap -R135:00/142:00/37:00/44:00 -JM10 -Ba5f1g1WSne -O -V >> D:\Myhome\yoreyore\gh90.eps

sample:

水深図

佐渡周辺
awk '{print $6, $5, -$9}' gh90.bs2.data > sado.dep.data
awk '{print $6, $5, -$9}' gh91.bs2.data >> sado.dep.data
として以下のファイルを作り、Nodataを手動で削除する。
GMTを起動して500mメッシュの地形データと重ねがきをする。
sample:

フリーエアー異常図（交点補正なし）

佐渡周辺
awk '{print $6, $5, $14}' gh90.bs2.data > sado.fag.data
awk '{print $6, $5, $14}' gh91.bs2.data >> sado.fag.data
として以下のファイルを作り、Nodataを手動で削除する。

GMTを起動して500mメッシュの地形データと重ねがきをする。

  makecpt -Chaxby -T-100/100/10 -Z > ghfag.cpt
  blockmean D:\GHdata2008\sado.fag.data -R137:50/140:10/37:00/39:00 -I30c -V > ghfag.blk
  surface ghfag.blk -R137:50/140:10/37:00/39:00 -I30c  -T0.2 -Gghfag.grd -V
  grdmask ghfag.blk -R137:50/140:10/37:00/39:00 -I30c -NNaN/1/1 -S10k -Gghfag.mask.grd -V
  grdmath ghfag.grd ghfag.mask.grd OR = ghfag.masked.grd

  grdimage ghfag.masked.grd -R137:50/140:10/37:00/39:00 -JM15 -Cghfag.cpt -P -K -V > ghfag.eps
  grdcontour bathy500m.grd -R137:50/140:10/37:00/39:00 -JM15 -C500 -K -O -V >> ghfag.eps
  psxy D:\GHdata2008\sado.fag.data -R137:50/140:10/37:00/39:00 -JM15 -Sc0.02 -G255/0/0 -K -O -V >> ghfag.eps
  psscale -D5/-1/8/0.2h -Ba50f10g5:Free-airAnomaly(mgal): -Cghfag.cpt -K -O -V >> ghfag.eps
  pscoast -R137:50/140:10/37:00/39:00 -JM15 -Df -W2 -K -O -V >> ghfag.eps
  psbasemap -R137:50/140:10/37:00/39:00 -JM15 -Ba0.5f0.25g0.25WSne -O -V >> ghfag.eps

sample:

DGRF異常図（交点補正なし）

佐渡周辺
awk '{print $6, $5, $11}' gh90.bs2.data > sado.mag.data
awk '{print $6, $5, $11}' gh91.bs2.data >> sado.mag.data
として以下のファイルを作り、Nodataを手動で削除する。
GMTを起動して500mメッシュの地形データと重ねがきをする。

sample:

参考文献

Susan Macmillan and Stefan Maus (2005)
International Geomagnetic Reference Field-the tenth generation. Earth Planets Space, Vol. 57 (No. 12), pp. 1135-1140, 2005
岸本清行・石原丈実・玉木賢策（1984）
地質調査所における海洋地球物理データ処理の現状。財団法人日本水路協会、シンポジウム資料、最近の海底調査―その４－、39-45
大門　肇・加藤　剛・片山真人・久保良雄（2002）
新・測地系変換ソフトウエアMGC2000A，水路部技報，vol.20, pp.58-70.
岸本清行（1999）
海陸を合わせた日本周辺のメッシュ地形データの作成：Japan250m.grd　地質調査所研究資料集（GSJ Open file Report）No.353
飛田幹男（2002）
世界測地系移行のための座標変換ソフトウエア"TKY2JGD",国土地理院時報,vol.97, pp.31-51．
上嶋正人・駒澤正夫（1996）
積丹半島の地磁気・重力・三成分地磁気異常測定。平成7年度研究概要報告―北海道南西沖海域、13-29、地質調査所
Wessel, P. and Smith, W.H.F. (1998)
New, improved version of the Generic Mapping Tools released, EOS Trans. AGU, 79, 579
地質調査総合センター(2005)
日本空中磁気データベース。数値地質図　P-6、地質調査総合センター
地質調査総合センター(2004)
日本重力 CD-ROM 第2版。数値地質図　P-2、地質調査総合センター
地質調査総合センター(2013)
日本重力データベース DVD版。数値地質図　P-2、地質調査総合センター
地質調査所・東、東南アジア沿岸地球科学計画調整委員会（CCOP）(1996)
400万分の1東アジア磁気異常図。数値地質図　P-1、地質調査所
中塚　正 (1989)
空中磁気探査のシステム化について(II) ―データ処理ソフトウェアシステム―. 地調月報, vol.40, 99-111.
中塚　正（1984）
空中磁気探査のシステム化について(I) ―ハードウェアシステム―. 地調月報, vol.35, 341-364.
中塚　正（2001）
日本周辺空中磁気異常のデータベース構築について。地質調査所月報、vol.52, No.2/3, pp.125-132

白嶺丸CR80形式のファイルから、の巻