日本空中磁気異常ＣＤの巻

• 格子点データファイルのデータフォーマット
• この磁気図DBでは，観測高度が異なる探査データを高度リダクションせずに取り扱っており，その格子点データファイルでは，各格子点での磁力値とその観測高度の両方のデータを収録している。 　　　　　　
• 磁力値の格子点データは，ヘッダー情報を示す2行と，それ以降の格子点データの羅列で構成される。 　　　　　　
• 格子は，用いている座標展開図法の北向き座標軸，東向き座標軸に平行に固定されており，nT単位の格子点データが，全格子点を南西角からはじめて北方向に向う順序で，FORTRAN言語プログラムの
  　　　　　　 do 1 k=1,my
  　　　　　　 write(1,'((f7.1,9(1x,f7.1)))') (v(i,k),i=1,mx)
  　　　　　1 continue
  に相当する形式で列記されている。
  
• このデータ部分はFORTRANプログラムの
  　　　　　　 read(1,*) ((f(i,k),i=1,mx),k=1,my)
  で読み出すことができる。
  
• 観測高度のデータは，これと全く同じ形式（単位は nT を m に読み替える）で，磁力値データに引き続いて記録してある。
  すなわち，ファイル全体の構成は，『ヘッダー情報／磁力値格子点値(nT単位)／ヘッダー情報／高度格子点値(m単位)』となっている。

• 呼称	位置	桁数	形式	内　容
  (ヘッダー 1行目)
  area	0	8	A8	地域名略称等
  nc	8	4	I4	座標展開図法番号． 各地方の200mメッシュデータでは 「UTMゾーン番号＋200」，全国図の1kmメッシュデータでは「300」 （ランベルト正積方位図法を意味する）になっている．
  -	12	4	4X	（空白）
  iorg	16	8	I8	座標展開図法の原点情報，ispa＝ ispb＝0． 各地方図データでは，iorg＝0 で korg は中央経線の経度を分単位で表したもの，全国図データでは， iorg と korg が座標原点の緯度・経度を分単位で表 したものとなっている．
  korg	24	8	I8
  ispa	32	8	I8
  ispb	40	8	I8
  (ヘッダー 2行目)
  ixs	0	12	I12	格子の南西角の北向き(ixs)・ 東向き(isy)座標値（m単位）．UTM座標では，東向き座標値は中央経線で 500,000m である．
  iys	12	12	I12
  mszx	24	6	I6	北方向(mszx)と東方向( mszy)の格子間隔（m単位）．mszxとmszy を別の値に設定できる形になっているが，同じ値が埋められている．
  mszy	30	6	I6
  mx	36	6	I6	北方向(mx)と東方向(my)の格子点数 （両端を含む）
  my	42	6	I6
  vnul	48	8	F8.1	格子点に有効データが存在しないことを示す特別な値（9999.9)

• cm123ALL(area)	251(nc)		0(iorg)	7380(korg)	0(ispa)	0(ispb)
  2600000(ixs)	425000(iys)	200(mszx)	200(mszy)	3501(mx)	3376(my)	9999.9(vnull)
  9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9
  9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9 9999.9

１．上記のE.GRDというファイルを例にとればUTM51ゾーンで200m間隔のデータになっている。 まずはUTMの座標上の値でeast,north,zという形式に変換する。
オリジナルデータでは北向きがXになっているので注意。

      character area
      dimension ff(5000,5000)
C        --------------
         open(unit=10,file='D:\REG\REG\E.GRD',status='old')
         open(unit=20,file='D:\REG\REG\E.200m.txt')
C        --------------
      read(10,111) area,nc,iorg,korg,ispa,ispb
      read(10,112) ixs,iys,mszx,mszy,mx,my,vnul
  111 format(a8,i4,4x,i8,i8,i8,i8)
  112 format(i12,i12,i6,i6,i6,i6,f8.1)
      write(6,*) mszx,mszy,mx,my
      if(mx.gt.5000 .or. my.gt.5000) then
         write(6,*) 'too large!'
         go to 9
      end if
c
        read(10,*) ((ff(i,k),i=1,mx),k=1,my)
c
      write(20,*) nc
      do 71 k=1,my
         kest=iys+mszy*(k-1)
      do 72 i=1,mx
        zz=ff(i,k)
       if(zz .ge. vnul) go to 72
        north=ixs+mszx*(i-1)
        xcor=float(kest)
        ycor=float(north)
      write(20,*) xcor,ycor,zz
C     ----------------------------------------
   72 continue
   71 continue
    9 close(10)
      close(20)
      STOP
      END

データが大きいので読み込めない可能性が高い。 そのときはパワーのあるパソコンを持っている人にやってもらう。

2．経度、緯度、値に変換する
UTMから緯度経度への変換は駒澤さんのIKDOを使った（日本重力CD-ROM第２版にプログラムあり）。wgs84版。
作成したnansei.grd.xyzは｛経度、緯度、磁気異常｝の並びで200mメッシュ相当。
xylltm.grd.fのメインルーチン

      IMPLICIT REAL*8(A-H,O-Z)
C        --------------
         open(unit=10,file='E.200m.txt',status='old')
         open(unit=20,file='nansei.grd.xyz')
C        --------------
      read(10,*) nc
      write(6,*) nc
c
      olat=0.0
c     zone 55=14700.0; 54=14100.0; 53=13500.0; 52=12900.0; 51=12300.0
      if(nc.eq.255) olon=14700.0
      if(nc.eq.254) olon=14100.0
      if(nc.eq.253) olon=13500.0
      if(nc.eq.252) olon=12900.0
      if(nc.eq.251) olon=12300.0
      QLAT = DEG(OLAT)
      QLON = DEG(OLON)
c
    1 continue
        read(10,*,end=71) xcor,ycor,zz
        xcor=xcor*1000
        ycor=ycor*1000
C     ----------------------------------------
      CALL  IKDO(YCOR,XCOR,PLAT,PLON,QLAT,QLON)
C     -----------------------------------------
      WRITE(20,600) PLon,PLat,zz
  600 format(f12.6,1x,f11.6,1x,f8.1)
      go to 1
   71 continue
      close(10)
      close(20)
      STOP
      END

3．全国１ｋｍメッシュデータも変換（例：wj1.grd→北海道）
ファイル形式は200mメッシュと同じなので同様の作業をする。

図法（ランベルト正積図法）で1000m間隔のデータになっている。
座標原点の緯度(I8)35degN;35*60=2100
座標原点の経度(I8)135degE;135*60=8100
まずはランベルトの座標上の値でeast,north,zという形式に変換する。
オリジナルデータでは北向きがXになっているので注意。
データが大きいので読み込めない可能性が高い。
そのときはパワーのあるパソコンを持っている人にやってもらう。

４．次にこれを、経度、緯度、値に変換する。

• ランベルト等積方位図法から緯度経度への変換は 中塚さんのライブラリーを使った。
• 作成したwjm0.1km.xyzは｛経度、緯度、磁気異常｝の並びで１ｋｍメッシュ相当。
• 同様に高度データ部分も呼び出し、｛経度、緯度、海面高度｝の並びで１ｋｍメッシュ相当のwjh0.1km.xyzを作成。

  ５．全国１ｋｍメッシュデータも変換（例：wj1.grd→北海道）
  ファイル形式は200mメッシュと同じなので同様の作業をする。
  
  

  hkdALL(area)	300(nc)		2100(iorg)	8100(korg)	0(ispa)	0(ispb)
  750000(ixs)	300000(iys)	1000(mszx)	1000(mszy)	451(mx)	601(my)	9999.9(vnull)
  183.0 216.7 240.1 239.8 269.8 255.0 211.5 182.9 155.6 135.6
  142.6 156.9 150.1 110.3 73.0 50.6 27.8 5.2 -16.1 -41.5

１．基本的に重力異常とおなじ。
　　　　

２．データがないところも多いので通常は再グリッドして用いる。（6秒メッシュ）
　　　　

３．200mメッシュのデータファイルは大きいので、グリッドファイルの作成は必要なところだけその都度やるほうが楽。
　　　　

４．一気にやるにはパソコンのパワーが足りなかった。

blockmean D:\hkd.grd.xyz -R142:27:00/142:49:00/42:05:00/42:26:40 -I6c -V > airmag.blk

surface airmag.blk -R142:27:00/142:49:00/42:05:00/42:26:40 -I6c -T0.2 -Gairmag.grd -V

grdmask airmag.blk -R142:27:00/142:49:00/42:05:00/42:26:40 -I6c -NNaN/1/1 -S4k -Gairmag.mask.grd -V

grdmath airmag.grd airmag.mask.grd OR = airmag.masked.grd

CD-ROMの解説文にあるように、この磁気図DBでは，観測高度が異なる探査データを高度リダクションせずに取り扱っている。

• 北海道・東北（A.GRD)
• 東北・関東（B.GRD)
• 中部・近畿（C.GRD)
• 九州（D.GRD)
• 沖縄（E.GRD)
• 観測高度１ｋｍメッシュ

５．自分のためのメモ
１ｋｍメッシュのデータ作成→airmag.masked.grd

makecpt -Chaxby -T-800/800/10 -Z > ghmag.cpt

blockmean I:\airmag-db\wjm0.1km.xyz -R120:00/150:00/20:00/46:00 -I30c -V > airmag.blk

surface airmag.blk -R120:00/150:00/20:00/46:00 -I30c -T0.2 -Gairmag.grd -V

grdmask airmag.blk -R120:00/150:00/20:00/46:00 -I30c -NNaN/1/1 -S4k -Gishimag.mask.grd -V

grdmath airmag.grd ishimag.mask.grd OR = airmag.masked.grd

grdimage airmag.masked.grd -R120:00/150:00/20:00/46:00 -JM15 -Cghmag.cpt -P -K -V > airmag.eps

psscale -D5/-1/8/0.2h -Ba250f50g50:MagneticAnomaly(nT): -Cghmag.cpt -K -O -V >>airmag.eps

pscoast -R120:00/150:00/20:00/46:00 -JM15 -Df -W2 -K -O -V >> airmag.eps

psbasemap -R120:00/150:00/20:00/46:00 -JM15 -Ba5f1g5WSne -O -V >> airmag.eps