PDS形式地図投影画像の取扱い

はじめに

Planetary Data System (PDS)とは、アメリカのNASAが運営している惑星科学分野のデータシステムである。主に月惑星探査機の観測データや、そこから派生して生成されたデータ、あるいは関連する実験室データをアーカイブし、公開している。

PDSとはまたこのデータシステム上で公開されているデータのフォーマットの名でもある。観測データ本体にメタデータ (PDS label) を付与する形式となっている。現在はメタデータをXMLで記載するPDS version 4 (PDS4) に移行しているが、少し前までは別の形式でメタデータを作成するversion 3 (PDS3) が主流であり、このバージョンのフォーマットで公開されているデータも多い。

PDSフォーマットは惑星探査データのデファクトスタンダードデータであり、日本の惑星探査ミッションでもデータアーカイブのフォーマットとして採用している。

ここでは月探査機「かぐや」（SELENE）に搭載されていた光学観測システムLISMの地図投影済みカメラデータを例に、PDS3形式の地図投影画像データを取り扱う方法を紹介する。

月面撮像・分光装置 (Lunar Imager/Spectrometer, LISM)

LISMは月探査機「かぐや」搭載されていた光学観測システムであり、

地形カメラ (Terrain Camera, TC)
マルチバンドイメージャ (Multiband Imager, MI)
スペクトルプロファイラ (Spectral Profiler, SP)

の三つの機器から構成されている。

データ入手

「かぐや」のデータはかぐやデータアーカイブから入手できる。

本稿で取り扱うのは「かぐや」一般公開プロダクト一覧にあるデータのうち、「マップデータ」と呼ばれる地図投影済みのデータである。例えば地形カメラのTCオルソマップなどが該当する。

フォーマット変換

gdalを使うのがよい

GDAL: https://gdal.org
- GDAL is a translator library for raster and vector geospatial data formats that is released under an X/MIT style Open Source License by the Open Source Geospatial Foundation.

Macの場合，MacPorts/homebrewどちらにもGDALパッケージがある．

macports: https://ports.macports.org/port/gdal/details/
homebrew: https://formulae.brew.sh/formula/gdal

gdalパッケージの内容とMacPorts/homebrewの違い

gdalは単体のコマンドと，Pythonスクリプト群から構成されている．
homebrewでは両方一度にインストールされるが，MacPortsでは単体のコマンドしかインストールされない．

PDS (.img) to FITS

基本的な書式は以下の通り．

# labelがimgファイル内に記述されている場合
$ gdal_translate -of FITS xxxxx.img xxxxx.fits

# lblファイルが別に存在する（detached）場合
$ gdal_translate -of FITS xxxxx.lbl xxxxx.fits

この方法で変換したFITS形式画像をds9などのFITSビューアで開く場合，上下反転は不要．

FITS以外の画像フォーマットにも変換することができる。GeoTIFFなどであれば地図投影情報も保存される。

SELENE/LISMデータの場合

SELENE/LISMデータの多くはimgファイルの画素値がunsigned形式になっている．この場合FITS変換後に正しい画素値を得るためには以下のようにする必要がある．

$ gdal_translate -unscale -of FITS TCO_MAP_02_S42E348S45E351SC.lbl TCO_MAP_02_S42E348S45E351SC.fits

これは，FITSは規格上signed形式のみしか存在せず，gdal_translateがこのことを正しくハンドリングできないためである．

なお，PDS labelに

  OFFSET = 0.000000
  SCALING_FACTOR = 0.018300

のようにオフセットやスケール値が定義されている場合は，

$ gdal_translate -a_scale 0.0183 -a_offset 599.6544 -of FITS TCO_MAP_02_S42E348S45E351SC.lbl TCO_MAP_02_S42E348S45E351SC.fits

とすればオフセットとスケール値が反映された画素値を持つFITS画像を生成できる．この時のオフセット値599.6544は以下のように算出する．

\begin{align}
a_{offset} &= 32768 \times a_{scale} \\
&= 32768 \times 0.0183 \\
&= 599.6544 \\
\end{align}

モザイク（複数画像の貼り合わせ）・地図投影法変換

これらの作業は地図投影法の情報が重要となるので，GeoTIFFまたは元のPDS形式で作業するとよい．FITS形式を経由すると問題が生じるかもしれない（FITS形式のデータが必要な場合は最後に変換するとよい）．

モザイクも，投影法変換もgdalwarpまたはgdal_merge.pyで行える．

モザイク（複数画像の貼り合わせ）

$ gdalwarp srcfile(s) dstfile

srcfile(s)で列挙した画像を貼り合わせて，dstfileとして出力する．

地図投影法変換

$ gdalwarp -t_srs srs_def srcfile(s) dstfile

画像srcfile（複数指定した場合は貼り合わせも行う）をsrs_defで指定した地図投影法でdstfileとして出力する．

srs_defの記述方法は少々ややこしい．
https://spatialreference.org/ に例が載っているので，それを元にすることができる．

IAUによる太陽系天体での定義例
月の投影法定義の例（上記のサブセット）

Well Known Text (WKT)形式（拡張子は.prj）のファイル（上記サイトでダウンロードできる）を読み込ませることができる．

Moon Equidistant Cylindrical形式で投影したい場合の例は，ここからダウンロードできるWKT形式のファイル30110.prjがカレントディレクトリに存在するとして，

$ gdalwarp -t_srs 30110.prj srcfile(s) dstfile

とする．投影中心などパラメータを持つ投影法を用いる場合は，WKT形式のファイルを編集して使うことができる．

Moon Transverse Mercator（横メルカトル図法）で投影中心経度を東経350度にする場合は，WKT形式のCentral_Meridianを350とする．

PROJCS["Moon_Transverse_Mercator_E350",
    GEOGCS["Moon 2000",
        DATUM["D_Moon_2000",
            SPHEROID["Moon_2000_IAU_IAG",1737400.0,0.0]],
        PRIMEM["Greenwich",0],
        UNIT["Decimal_Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",0],
    PARAMETER["False_Northing",0],
    PARAMETER["Central_Meridian",350],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0],
    UNIT["Meter",1]]

これを例えば30164_E350.prjとして，

$ gdalwarp -t_srs 30164_E350.prj srcfile(s) dstfile

とすると投影中心経度東経350度の横メルカトル図法にすることができる．

画像座標・地理座標変換

地図投影された画像中のある地点について画像座標がわかっている時，その地点の緯度経度（地理座標）は，地図座標を経由して以下のようにして得る．

画像座標→地図座標

地図座標は地図平面上の座標である．地図座標の原点は，地図投影の基準点であり，単位は画素のサイズである．地図投影済み画像の投影情報はgdalinfoコマンドで調べることができる．

$ gdalinfo test_tm_E350.tiff
Driver: GTiff/GeoTIFF
Files: test_tm_E350.tiff
Size is 20912, 14580
Coordinate System is:
PROJCRS["Moon_Transverse_Mercator_E350",
    BASEGEOGCRS["Moon 2000",
        DATUM["D_Moon_2000",
            ELLIPSOID["Moon_2000_IAU_IAG",1737400,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Decimal_Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["Decimal_Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",350,
            ANGLEUNIT["Decimal_Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (-112645.896858343825443,-1273067.585431055165827)
Pixel Size = (6.463906771405839,-6.463906771405839)
Metadata:
  AREA_OR_POINT=Area
  DATA_SET_ID="SLN-L-TC-5-ORTHO-MAP-V2.0"
  INSTRUMENT_ID="TC"
  INSTRUMENT_NAME="TERRAIN CAMERA"
  MISSION_NAME="SELENE"
  PRODUCT_CREATION_TIME=*
  PRODUCT_ID=*
  START_TIME=UNK
  STOP_TIME=UNK
  TARGET_NAME="MOON"
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -112645.897,-1273067.585) ( 14d59'29.83"W, 41d53'29.91"S)
Lower Left  ( -112645.897,-1367311.346) ( 15d15'16.57"W, 44d59'18.38"S)
Upper Right (   22527.322,-1273067.585) (  9d 0' 0.00"W, 41d59'43.94"S)
Lower Right (   22527.322,-1367311.346) (  8d56'49.34"W, 45d 6'15.28"S)
Center      (  -45059.288,-1320189.466) ( 12d 3' 2.04"W, 43d32'10.15"S)
Band 1 Block=20912x1 Type=UInt16, ColorInterp=Gray
  NoData Value=0

この中の

Origin = (-112645.896858343825443,-1273067.585431055165827)

は，画像座標原点の座標を地図座標（原点地図座標）で表現したものである（地図座標の原点ではないことに注意）．また，Corner Coordinatesとして画像の4隅の地図座標値も表示される．

地図座標系の単位である画素サイズは

Pixel Size = (6.463906771405839,-6.463906771405839)

である．Y軸方向の画素サイズが負なのは，通常地図座標は北が上（つまり+Y方向が上）なのに対し，画像座標では+Y方向が下であることによる．

地図画像上で画像座標がわかっている時，その場所の地図座標は

$$
\begin{align}
地図座標 = ([Origin] + [画像座標]\times [Pixel Size])
\end{align}
$$

として求めることができる（Pixel Sizeの正負はそのまま使うこと）．

地図座標→地理座標

gdaltransformで地図座標から地理座標を得ることができる．実際にはこのコマンドは2種の地図投影法の間の地図座標の変換を行うので，変換先として地理座標系を指定すればよい．

月の場合Moon 2000を用いる．ここからダウンロードできるWKT形式のファイル30100.prjをカレントディレクトリに置いておく．

先に見た投影中心経度東経350度の横メルカトル図法（prjファイル30164_E350.prj）で画像左下（地図座標Lower Left ( -112645.897,-1367311.346)）の地理座標は，

$ gdaltransform -s_srs 30164_E350.prj -t_srs 30100.prj
Enter X Y [Z [T]] values separated by space, and press Return.
-112645.897 -1367311.346
-15.2546021973711 -44.9884375919914 0

として，西経14.99度，南緯44.99度とわかる．gdalinfoの出力

Lower Left  ( -112645.897,-1367311.346) ( 15d15'16.57"W, 44d59'18.38"S)

とも一致している．

gdaltransformは対話形式で動作し，地図座標を入力するたびに変換後の座標が表示される．コマンドラインで非対話的に結果を得たい場合は，

$ echo "-112645.897 -1367311.346" | gdaltransform -s_srs 30164_E350.prj -t_srs 30100.prj
-15.2546021973711 -44.9884375919914 0

のようにする．地図座標が1行に一つ記載されたテキストファイルがあれば，

$ gdaltransform -s_srs 30164_E350.prj -t_srs 30100.prj < input.txt

で一括変換することもできる．出力もリダイレクトでファイルに保存できる．

$ gdaltransform -s_srs 30164_E350.prj -t_srs 30100.prj < input.txt > output.txt