QGISを用いた古代地図の復元(1) QGISの導入と基本地図作成

　この「古代史調査ナビ」の目的の1つは、古代史調査のためのさまざまなツールを紹介することである。「古墳GISデータの集め方」では、自治体等からGISデータを収集する方法について詳述した。構築した「GISシステム」にこれを取り込めば、古墳の分布状況が視覚化される。この「GISシステム」として「地理院地図」を想定していた。ある程度慣れたところで、フリーのGISソフトウエアである「QGIS」に進む計画だった。しかし、本格的なGISソフトである「QGIS」を導入してしまえば、「地理院地図」を経由する必要性はないのではないかと思うようになった。最初のハードルは高いが、いったん入ってしまえば「QGIS」が生み出す世界は「地理院地図」より遙かに広い。

　本来「古代史調査」において、GISのベースとなるのは「古代の地形を復元した地図」である。縄文・弥生・古墳時代の地形は現在のものとは相当異なっている。だから、「古代史」を捉える際には、その当時の地図が必要である。「香取海の地図」では、当時の状況を再現しようと、「地理院地図」を用いて古代の汀線を再現すべく「標高」を5ｍ、3ｍと変えてみた。今の海面より何メートル上げるかについてはとくに明確な根拠があるわけではなかった。だから、あくまでも「個人的推測」の域を出ない。もうひとつ、この方法である程度の古代の状況はつかめたとしても、海面変化だけでなく、河川による土砂の堆積や地盤沈下などによる地形の変化も起こっているから、どこでも同じ標高が使えるわけではない。加えて、人為的な埋め立て・河川改修操作もある。「古代地図の復元」を行おうとすると、これら複数の要素に基づく様々な変化を組み込む必要がある。こうなると、「地理院地図」ではとても無理である。やはり本格的なGISソフトが必要だと思う。

　このような観点から「QGIS」の導入を決めた。「QGIS」の機能を100％使って「復元古代地図」を作成する。それが、この一連のブログ「QGISを用いた古代地図の復元」のゴールである。

　分厚い本を一冊買った。喜多耕一氏の『業務で使うQGIS ver3 完全使いこなしガイド』（以下、「マニュアル本」と呼ぶ）である。それに加えて、Gemini Notebookに「QGISを用いた古代地図の復元」という新しいノートブックを作成した。

「古代地図の復元」はかなり大変な作業である。素人の私にはどこから手をつけていいか分からない。そこでAIの登場だと考えた。まず、AIにこれまで「それほどの根拠なしに単に今の海面より何メートル上げると仮定して地図を作ってきた」と告げ、「QGISを使ってもっと科学的にやりたいが、いい方法はないか？」と尋ねた。すると、「地質学、微地形学、考古学などのオープンデータを用いて、QGIS上で科学的根拠（エビデンス）に基づいて再現する手法が確立されている」という答えが返ってきた。さらに「科学的に正確な古代の海岸線（汀線）の復元は、現在の標高値（DEM）を一律にスライドさせるだけでなく、人間の手による『人工改変の剥ぎ取り』、氷期に刻まれた『埋没谷の復元』、そして『遺跡・貝塚データの分布に基づく制限』や『地域的な地盤隆起・沈下速度の補正』を有機的に組み合わせるQGISの多レイヤ分析手法が極めて有効であることが分かりました」といい、関係する論文のリストと検討方針まで作成してくれたのである。

　かなりのハルシネーションが含まれているかもしれない（私は疑り深くなっている！）が、やってみて損はない。しかし、AIの力を借りるだけで果たして素人の私が目的の「古代地図」を作れるのだろうか？　一抹の不安を覚えながら「QGIS」のインストールサイトを開く。さあ、同時進行ドキュメントの始まりである！

QGISのインストール

　Gemini Notebookにソースとして追加すべきものを訊いてみると、

QGIS公式ダウンロード・インストールガイド
QGIS 初心者向け基本操作・入門マニュアル
国土地理院データ（DEM・基盤地図情報）のQGIS読み込み手順
プラグイン（Qgis2threejs 等）の導入・操作ガイド

　を挙げてきたので、AIに任せてこれらをソースに追加した。続いて、 QGISのダウンロードとインストールである。「推奨バージョン: LTR版（Long Term Release：長期サポート版）のダウンロードをおすすめします」とあったので、それを選んでダウンロード。「QGIS-OSGeo4W-3.44.14-1.msi」をダブルクリックしてインストール。「無事にQGISを起動できましたら教えてください」とGemini Notebookから指示があった。

　無事インストール完了。デスクトップに「QGIS 3.44.14」というフォルダができていた。この中の「QGIS Desktop 3.44.14」をダブルクリックする。

QGISの起動と地理院タイルの表示

　初期画面は図1。ここでGemini Notebookに起動できたことを伝える。マニュアル本に書かれていた画面の描画速度を速くする「レンダリング設定」はすでに「レンダリングに使うCPUコア最大数」が「すべて」になっていたので設定不要だった。

　Gemini Notebookから操作画面上に背景地図となる国土地理院の「地理院タイル」を表示させてみましょうと指示がある。これはマニュアル本の「第5章インターネット地図を使う」にあたる。インターネット地図だとパソコンがネットに繋がっていないと使用できないが、フィールドで使うわけではないので十分だろう。

　画面左側のブラウザパネル内にある［XYZ Tiles］という項目を探して右クリックし、［新規接続…］を選択。設定画面に、名前「地理院地図（淡色）」、URLにhttps://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.pngと入力。（標準地図や白地図を表示したければ、それぞれに対応するURLを入力すればよい。）

　すると、ブラウザパネルの「XYZ Tiles」の中に「地理院地図（淡色）」が追加され、ここをダブルクリックするとレイヤのパネルに「地理院地図（淡色）」が現れた。これを右クリックして「レイヤの領域にズーム」を選ぶと、右側に地図が表示される。と、ここまでとても順調。

ここで注意事項がある。

地理院地図を利用する場合には、出典を明示する必要がある
その地図が「基本測量成果」である場合には測量法の申請が必要。標準地図や淡色地図などはこれに該当する。
ただし、「基本測量成果」のところにつぎのような記載がある。「地理院タイルをウェブサイトやソフトウェア、アプリケーション上でリアルタイムに読み込んで利用する場合、地理院タイルは出典の明示のみで申請不要でご利用いただけます。出典は、『国土地理院』または『地理院タイル』等と記載していただき、地理院タイル一覧ページ（https://maps.gsi.go.jp/development/ichiran.html）へのリンクを付けてください。」

　QGISでインターネット地図の「地理院タイル」を使う場合はこの3に該当するため、出典の記載のみでよく、申請は不要だ。なお、QGISのマップビューに地図の出典を表示させることができる。メニュー→ビュー→地図装飾→著作権ラベルで設定する。

　播磨地方を表示させたのが図2。なお、この図では左のレイヤパネルに「地理院地図（淡色）」の他のレイヤが見えている（非表示にしてある）が、これはあとから追加したもので、この段階では「淡色地図」だけが表示されていた。

標高データ（DEM）の導入

「古代の海岸線や入り江（溺れ谷）を可視化するには、土地の高さを示す『標高データ（DEM: 数値標高モデル）』が必要となります。国土地理院が無料公開しているDEMデータ（XML形式）をQGISで扱えるGeoTIFF形式に変換するため、専用プラグイン『QuickDEM4JP』をインストールしましょう」と指示が出る。

　DEMとは「Digital Elevation Model」の略で、地表面をメッシュで区切って、その位置の標高を数値化（ｍ）した標高のメッシュデータである。画像データ（ラスタデータ）として提供される。国土地理院の基盤地図情報から入手できるものがXML形式なので、QGISが読めるGeoTIFF形式に変換する必要があるのだ。この変換にここでは専用プラグインの「QuickDEM4JP」（株式会社MIERUNE）を使うが、マニュアル本では株式会社エコリスの「基盤地図情報標高DEMデータ変換ツール」を紹介している。

　メニューから［プラグイン］＞［プラグインの管理とインストール…］をクリックすると、プラグインの一覧が表示される。検索窓に「QuickDEM4JP」と入れ、表示されたプラグインをインストール。すると、ツールバー上に「XML」と書かれたアイコンが出現！

　つぎに「国土地理院」からDEM（標高データ）をダウンロードするため、ブラウザで「基盤地図情報ダウンロードサービス」（https://service.gsi.go.jp/kiban/）にアクセス。「数値標高モデル」を選択。対象地域を選択して、［まとめてダウンロード］（ZIP形式）を行う。なお、ダウンロードには無料の利用登録・ログインが必要となる。1m、5m、10mメッシュがあるが、5ｍでもデータ容量が大きすぎるので10mメッシュを選択。兵庫県と岡山県をカバーする範囲だけでも、容量合計が430.86MBもあった。

　つぎに「QuickDEM4JP」にダウンロードしたファイルを読み込ませる。［XML］アイコンをクリックして開いたボックスの「入力設定」にダウンロードした ZIPファイルを指定する。出力設定で「ファイルに保存」を選び、保存先フォルダとファイル名を設定。ここで実行を押すと、XMLがGeoTIFFに変換されるはずなのだが、すぐにエラーが発生した。どうやらダウンロードしたままのZIPファイルではダメなようだ。ZIPファイルを解凍してみると、中にもZIPファイルがたくさん並んでいる。そのひとつを試しに解凍してみるとXMLファイルが1つ現れた。どうやら、XMLファイルを1つずつZIP化し、その全体をまたZIP化するという手の込んだ構成になっているようだ。この解凍だが、「7-Zip」を使えば簡単にできる。まず、ダウンロードしたZIPファイルを普通に解凍して、フォルダの中に複数のZIPファイルが含まれている状態にする。「7-Zip」をインストール後、すべてのZIPファイルを選択して右クリック、メニューの一番下にある「その他のオプションを表示」をクリックすると「7-Zip」が現れるので、これをクリックすると全てのZIPファイルを一発で解凍できた。

　解凍して得られたXMLファイルを一つのフォルダに入れて、そのフォルダを「QuickDEM4JP」の「入力設定」で指定しようとしたがうまくいかない。よく見ると、「XMLを含むZIP/XML」と記載がある。つまりXMLが入った単なるフォルダは使えないのだ。XMLを一つ一つ指定していては日が暮れるので、ここはZIP化しかない。なんとせっかく解凍したXMLファイルをまた圧縮するという手間のかかることになってしまった。もうちょっとなんとかならないものか！　

　さて、このZIPファイルを読み込ませ実行させるとまたエラー。大多数のファイルは頭に「FG-GML」がついているのだが、他にも「fmdid」というファイルがあった。これも一緒にZIP化したのだが、どうもそれがいけなかったらしい。これらは国土地理院のメタデータ、スキーマ定義ファイル、基本項目データなどで、本来のDEMデータではないのでエラーになったようだ。これらを削除して、再度ZIP化して実行。なんとか最後まで読み込んでくれた。画面は図3のように変化！　このように真っ黒だが、つぎで標高別に色分けする作業を行う。

段彩図を作る

　標高DEMのレイヤをダブルクリックするとウィンドウが開く。サイドメニューからシンボロジを選んで、図4のように設定してみた。

一番上の「レンダリングタイプ」が「単バンドグレー」になっていたので、「単バンド疑似カラー」に切り替え
最大値/最小値設定の補間方法のところは、「線形」「離散」「厳密」と切り替えられるが、「地理院地図」の「自分で作る色別標高図」の設定と同じになるように、範囲内が同じ色となる「離散」にした
色分けは7段階。左の標高の数値をダブルクリックすると数字が変えられるようになり、色をダブルクリックすると配色が変えられるようになる。配色は図5のようにした。ここでは標高4ｍまでを海の色と同じにしている
ここで左下の「スタイル」を押して、メニューからスタイル保存を選び、保存しておく
最後に左のサイドメニューから透明度を選び、70％に設定

　適用→OKと押して実行させると図6のように、海・平地・山が色別されて表示された。「段彩図」の完成である。

陰影図の追加

「段彩図」だけではもう一つ地形の状況がハッキリしないので、「陰影図（hillshade）」を追加してみよう。地形の影（起伏）を生成して重ねることで、台地の崖や谷の立体感をリアルに強調できる。その手順をAIに訊くと、以下の様にすれば良いと回答があった。

QGISの上部メニューから［ラスタ］＞［解析］＞［陰影図（Hillshade）…］をクリック
開いた画面で設定を確認：
- 入力レイヤ: 先ほど作成した標高データ（GeoTIFF）が選択されていることを確認
- Z係数：1～2（大きくすると高さが強調される）
- その他の設定はそのままでよい
右下の ［実行］ ボタンをクリック
処理が終わったら［閉じる］をクリック
画面左側の「レイヤ」パネルに 陰影起伏 という新しいレイヤが追加され、画面にはモノクロ（白黒）の立体的でリアルな山肌や谷の影が表示される

　この通り実施してみると、確かに「陰影図レイヤ」が追加された。レイヤパネルでの順番を上から、標高DEM、陰影起伏、淡色地図の順にする。さらに、

陰影起伏レイヤの透明度を50％に変更
陰影起伏レイヤをダブルクリックしてシンボロジを開き、レイヤレンダリングを「乗算」に変更

　と指示の通りに進んだ結果が図7だ。起伏の表示が滑らかではなく、多数の線が現れており、予想していたものとはかなり違う。

　実はこの問題の解決に結構時間を費やした。AI の指示通りにいろいろやったのだが、一向に問題が解決しない。マニュアル本を調べてみると、「傾斜区分図」の「ひとことメモ」に「作成した傾斜図にスジのような模様が入ってしまう場合があります。これは標高DEMの座標参照系がWGS 84/Pseudo-Mercator（EPSG:3857）になっている場合に起こります。標高DEMをレイヤのエクスポートで座標系を平面直角座標系に変更してから傾斜図を作りましょう」と書かれていた。関係しそうだが、これはいったい何だろうか？　Gemini Notebook に質問を投げかける。

「hillshadeに線がたくさん現れていてスムーズではない。座標系の問題だと本にありましたが、どうしたらいいですか？」

「本に書いてあった通り、まさに『座標系（CRS）』の違いが原因です！　国土地理院からダウンロードした直後のDEMデータは、位置を『緯度・経度』で表す『地理座標系（単位：度 °）』（例: JGD2011 / EPSG:6668）になっています。しかし、標高（高さ）の単位は『メートル（m）』です。横方向が『度』で縦方向が『メートル』のまま陰影計算（Hillshade）を行うと、QGISが距離の計算を正しく処理できず、画面上に不自然な縦横の格子線（縞模様やノイズ線）が現れてしまいます。これを解決するには、DEMの単位をメートルに揃える『ワープ（再投影）』を行います。手順は以下の通り非常に簡単です！」という答え。

バカ！　もっと早く言ってくれ！

座標参照系

　座標参照系に関する詳しい説明は以下を見てもらいたい。マニュアル本では第4章でこの問題を扱っている。

8. 座標参照系 — QGIS Documentation ドキュメント

QGIS 3.44 documentation: 8. 座標参照系

　大きく分けて座標系には緯度・経度で位置を表す「地理座標系」と、球体の地形を平面に投影し原点からの距離をメートルで表す「平面座標系（投影座標系）」がある。地理院地図・Google Mapsなどは「地理座標系」を用いている。QGISでは「地理座標系」「平面座標系」ともに使用できるが、両方が混在すると時には正しく処理できない場合が生じるというわけだ。この地球上の曲面（球体や楕円体）の位置を、地図やコンピュータ上で正確に表すための基準と仕組みを「座標参照系（CRS: Coordinate Reference System）」と呼ぶ。

　各座標参照系にはそれぞれ名前があるが、「EPSG」という座標参照系コード番号を使用すると便利だ。QGISの画面の右下にはプロジェクトの座標参照系が表示されている。また、レイヤパネルで各レイヤを右クリックし「CRSを表示」で座標参照系を表示できる。プロジェクトの座標参照系が背景にWeb 地図をつかったので EPSG:3857（平面座標系）、また標高DEMレイヤはGeoTiff変換時のデフォルトのEPSG:4326（地理座標系）となっていた（図8）。AIが言うように「地理座標系」と「平面座標系」が混在したのでノイズが発生したようだ。

　ここで日本でよく使われる座標参照系の「EPSGコード」を整理しておこう。

地理座標系

EPSG:4326：WGS84　世界測地系、GPSで記録されるデータなどに使用
EPSG:6668：JGD2011　新日本測地系、2011年の東日本大震災の地殻変動を加味した座標参照系。WGS84とほぼ同じ

平面座標系

EPSG:3857：WGS84/Pseudo Mercator Webメルカトルと呼ばれ、Google Maps や地理院地図などのインターネット地図に使われている。投影座標系に分類されるが、あとの平面直角座標系とは異なり、高緯度・中緯度地域では距離や面積にかなりの拡大歪みが生じるので、面積計算など注意が必要
EPSG:6669～6687：日本測地系2011（JGD2011）平面直角座標系。精度が要求される場合に向いているが、一つの座標系でカバーされる地域は狭くなる。各コードに対応する地域については例えばLL.meを参照のこと。関心のある岡山県・兵庫県は EPSG:6673
EPSG:6688～6692：JGD2011/UTM zone 51N-55Nユニバーサル横メルカトル。世界を東西に経度 6度の幅を持つ 60等分ゾーンに分割している。投影法は平面直角座標系と同じ横メルカトル図法だが、ゾーンの幅が平面直角座標系より広いため、ゾーン端では歪みが大きくなる。そのため精度最優先の国内測量には平面直角座標系が適しており、UTMは広域データとの相互運用性を重視する場面に向いている。各コードに対応する地域については例えばLL.meを参照のこと。関心のある岡山県・兵庫県はEPSG:6690 （UTM zone 53N）

座標参照系変換そして再度陰影図の作成

　座標参照系についての知識を頭に入れたところで、「標高DEM」の座標参照系を「平面直角座標系」に変換しよう。このステップはつぎの通り。

ステップ1：DEMの座標系を「投影座標系」に変換する

上部メニューから［ラスタ］＞［投影法］＞［再投影（warp）］を開く
入力レイヤ: 作成したDEMデータ（GeoTIFF）を選択
ターゲットCRS（重要）: 右側の地球儀アイコンをクリックし、対象地域に合わせた平面直角座標系（メートル単位）を選択。6690で検索して、EPSG:6690 （UTM zone 53N）を設定した
リサンプリング方法（重要）:初期設定の「最近接法」から、「バイリニア」 に変更した。ピクセルのつなぎ目が綺麗に補間され、線や階段状のドット感が消えて滑らかになるらしい
出力レイヤ: ［…］をクリックし、ファイル名を指定。
［実行］。

ステップ2：変換したDEMから再度「陰影図（Hillshade）」を作成する

新しく作成された GeoTiff ファイルを使って、再度［ラスタ］＞［解析］＞［陰影図（Hillshade）…］を開く。
Zファクター: 1.0 〜 2.0 程度に設定。
［実行］をクリック。

　できた陰影図を比較してみよう（図9）。変更前（左）は0～100ｍの領域に線が現れていたが、変更後（右）はそれが消えてグレーとなり、起伏が見やすくなっている。

　二つのレイヤの座標参照系がEPSG:6690 （UTM zone 53N）となったが、プロジェクトの座標系はEPSG:3857（Webメルカトル）のままである。QGISのプロジェクトではいくつもレイヤを組み合わせて一つの地図を作る。プロジェクトの座標参照系とは、画面に読み込んだ地図やデータを、地球上のどの位置・どの形状で表示・計算するかを決めるための土台となるもので、レイヤの座標系が違っていても自動的にプロジェクトの座標参照系を基準として位置を合わせてくれる。問題はこのEPSG:3857（Webメルカトル）だ。これは高緯度になるほど面積が極端に拡大する歪みがでる。従って、測量に使う平面直角座標系とは精度・用途ともに別物であることを理解しておく必要がある。問題は既に記したように一つの座標系でカバーされる地域が狭いこと。ここではカバー範囲がより広いユニバーサル横メルカトル（UTM）をレイヤの座標参照系に選んだので、プロジェクトの方も合わせておこう。

　このことをGemini Notebookに訊いてみると、つぎのような回答があった。