この件に関しては、和歌山県オープンデータカタログサイトの担当者の方にお伝えし、現在は修正済みです。
はじめに
和歌山県オープンデータカタログサイトからダウンロードできる津波浸水想定図の Shapefile データが正しい位置で表示されない。という質問を目にしました。
質問の方は既に回答済みですが、私の方でも読み込んでみると再現され、配布されているデータの方が壊れていることが確認できました。
津波浸水想定図のデータだけでなく保安林など、ダウンロードできる他の Shapefile でも（一部を除き）同様の現象がみられます。
prj ファイル
より正確にいうと、使用している座標系の情報を格納している prj ファイルの記述がおかしなことになっていました。
まず、正解として EPSG:3857 WGS 84 / Pseudo-Mercator の標準的な定義のうち、 WKT1 によるものをみてみましょう。
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3857"]]
次に問題の prj ファイルをみてみます。
PROJCS["WGS 84 Web メルカトル",
GEOGCS["WGS 84",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree minute second hemisphere",0.017453292519943295,
AUTHORITY["EPSG","9108"]],
AXIS["Long",EAST],
AXIS["Lat",NORTH],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator"],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["E",EAST],
AXIS["N",NORTH],
AUTHORITY["EPSG","3857"]]
S["WGS 84 Web メルカトル",
GEOGCS["WGS 84",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree minute second hemisphere",0.017453292519943295,
AUTHORITY["EPSG","9108"]],
AXIS["Long",EAST],
AXIS["Lat",NORTH],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator"],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["E",EAST],
AXIS["N",NORTH],
AUTHORITY["EPSG","38
ともに、わかりやすく改行とインデントを加えております。
問題その1
津波浸水想定（南海トラフ）.prj の前半部分で、AUTHORITY["EPSG","3857"]] と称する座標系の情報を記述していますが、そのあと、間違ってコピペをしてしまったかのような PROJCS の S からはじまり、 AUTHORITY["EPSG","3857"]] の AUTHORITY["EPSG","38 で終わる謎の記述が入り込んでいます。
他の GIS ソフトはわかりませんが、 QGIS においては、特にエラー等は起こさず前半部分の読めるところまで読んで解釈したようです。
問題その2
後半部分に変な記述が入り込んでいること自体おかしなことですが、それはそれとして、前半部分の記述も不適切であり、それゆえ位置がずれて表示されています。
標準的な座標系定義とくらべて細かな差異がありますが、一番の問題は下記の記述が不足していることです。
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],
疑似メルカトル（ウェブメルカトル）は、計算を容易にするため +a=6378137 +b=6378137 にあるように赤道半径と極半径を同一の大きさ、すなわち真球とみなしたメルカトル図法で投影しています。（メルカトル図法に投影する元となる緯度経度を定める際には WGS84 楕円体を用いています。）
座標系の WKT 表記にも複数種類があり、本文中で引用しているのは WKT1 の記述法となります。 WKT2 や ESRIWKT 形式では EXTENSION["PROJ4", ...] の記述はありませんが、かわりに METHOD["Popular Visualisation Pseudo Mercator", ID["EPSG",1024]] や PROJECTION["Mercator_Auxiliary_Sphere"] の記述がなされ、それにより真球としてメルカトル投影を行う座標系であることを示しています。
津波浸水想定（南海トラフ）.prj にはそれに関する記述がないため、 WGS 84 Web メルカトル というラベルと AUTHORITY["EPSG","3857"] としながらも、真球ではない普通の（楕円体の）メルカトル図法の座標系定義になり、実質的に EPSG:3395 WGS 84 / World Mercator として読み込まれてしまいます。
その他、細かな疑問点
UNIT["degree minute second hemisphere",0.017453292519943295,
AUTHORITY["EPSG","9108"]],
確かに、単位として EPSG:9108 degree minute second hemisphere は登録されていますが、普通は使われないです。
座標値が degree ではなく dms であるのならば、使う分には別によいと思いますが、その単位を使っている GEOGCS の定義を既存の AUTHORITY["EPSG","4326"] としてよいのか。。。
原因の推察
津波浸水想定図だけでなく保安林、道路規制情報、埋蔵文化財包蔵地などの複数のデータでみられました。
ただし https://data.bodik.jp/dataset/... からダウンロードする大雨による主要水系河川の洪水浸水想定区域図は発生しませんでした。
津波浸水想定図などの問題になっているデータは、ウェブ地図上では正しく表示されています。
これより推察となります。実際に製品のシステムにより生成されたものであるかは不明です。
Shapefile のダウンロード先は
https://wakayamaken.geocloud.jp/webgis/OpenData/Download/Shp/Polygon/1/0 となっていることから、ウェブ GIS である GeoCloud のダウンロード機能を使っているのではないかと思われます。
おそらくデータ作成時は正しく作成され、ウェブ地図にインポートする際も問題なく行われたにも関わらず、 GeoCloud のエクスポート機能がポンコツな prj ファイルを作成しているのではないかと推察できます。
この Shapefile が、システムにより動的に生成されているものか、システムにより過去に生成したものなのか、他のシステムや業者が生成したものをダウンロードできるようにしているものなのかは不明です。
もしウェブ GIS のシステムが生成したものだとしたら、オープンデータとしては珍しく EPSG:3857 なデータとなっているのは、ウェブ GIS のエクスポート機能によるものだったからでしょうか？ （一般的には地理座標系か平面直角座標系）
別の自治体について
ほかに GeoCloud のダウンロード機能でオープンデータを配布している自治体を探したところ、甲府市をみつけました。座標系として EPSG:3857 もどきではなく、 EPSG:4326 でした。ただし単位の定義は degree ではなく dms でした。
GEOGCS["WGS 84",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree minute second hemisphere",0.017453292519943295,
AUTHORITY["EPSG","9108"]],
AXIS["Long",EAST],
AXIS["Lat",NORTH],
AUTHORITY["EPSG","4326"]]
甲府市のデータでは、一応ちゃんとした座標系定義になっているようです。ただ、記述内容のクセが同じであるため、和歌山県と同じシステムが生成した prj ファイルと思われます。
更に他の自治体
深谷市の例
GEOGCS["WGS 84",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree minute second hemisphere",0.017453292519943295,
AUTHORITY["EPSG","9108"]],
AXIS["Long",EAST],
AXIS["Lat",NORTH],
AUTHORITY["EPSG","4326"]]
^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@^@
^@ は NUL 文字 \x00 。
もしほんとにシステム側で生成しているとしたら、ちょっと杜撰すぎやしませんかね。。。
実験
（作成されたデータに対し）適用される座標系が不適切な場合の位置のズレについて、確かめてみようと思います。
作業の流れ
私自身も少し混乱したため、まず最初に整理しておきます。
EPSG:3857 Pseudo Mercator は赤道半径＝極半径として計算したメルカトルで、 WGS84 におけるちゃんとしたメルカトル投影になっていないので、疑似メルカトルと呼ばれています。
EPSG:3395 World Mercator は、 WGS84 楕円体を使用した楕円体考慮のメルカトル投影です。 GIS 的にはあまり使われませんが、地図学的には通常のメルカトル図法といえると思います。
作業の流れとしては
- 北緯$34^\circ$における地点を EPSG:3857 Pseudo Mercator で投影したときの座標を計算（データの準備）
- EPSG:3857 における座標値を、座標値の値はそのままで EPSG:3395 World Mercator であると誤認識
- EPSG:3395 World Mercator として認識してしまった場合、地球上のどの位置を示しているか確認するため緯度に逆変換
- 計算した緯度は、当初意図していた北緯$34^\circ$とは別の地点。その地点が EPSG:3857 上でどこに位置するか求めてみる
| ステップ | 緯度 | Pseudo Mercator | World Mercator |
|---|---|---|---|
| 1 | $34.000000000^\circ$ | $4,028,802.026$ m | |
| 2 | $4,028,802.026$ m | ||
| 3 | $34.178573659^\circ$ | $4,028,802.026$ m | |
| 4 | $34.178573659^\circ$ | $4,052,805.384$ m |
ここで得られるのは、あくまで緯度座標間あるいは EPSG:3857 Pseudo Mercator の座標値間の差異なので、楕円体上の距離も求めます。
QGIS
QGIS で差を確認するには
- EPSG:3857 で点レイヤを作成
- もし緯度経度で指定した点で計測したい場合は EPSG:4326 で点レイヤを作成
- 点を追加し、編集用の頂点ツールより頂点エディタで任意座標を指定
- 「ベクタレイヤを再投影」で EPSG:3857 に再投影
- 作成したレイヤを複製
- 複製したレイヤのプロパティを開き EPSG:3395 に設定
- 元のレイヤと複製したレイヤ間の距離を、計測ツールを利用して測る
- スナップを有効にするとよい
和歌山県付近の北緯34度だと EPSG:3857 の地図上座標（デカルト座標）で24,003.358 m 、地球上（楕円体）で19,808.105 m ほどの差が生じます。
Proj
コマンドだと Proj の cs2cs を使って
$ # 東経135度 北緯34度の地点を EPSG:3857 疑似メルカトルに再投影したときの座標値
$ echo 135 34 | cs2cs -f %.6f +init=EPSG:4326 +to +init=EPSG:3857
15028131.257092 4028802.026134 0.000000
$ # 東経135度 北緯34度の地点を EPSG:3857 Pseudo Mercator （疑似メルカトル）に再投影したときの座標値を求め
$ # その値を EPSG:3395 World Mercator （楕円体のメルカトル）とみなし、緯度経度に再投影
$ echo 135 34 | \
cs2cs -f %.6f +init=EPSG:4326 +to +init=EPSG:3857 | \
cs2cs -f %.12f +init=EPSG:3395 +to +init=EPSG:4326
135.000000000001 34.178573658865 0.000000000000
$ # 東経135度 北緯34度の地点を EPSG:3857 Pseudo Mercator （疑似メルカトル）に再投影したときの座標値を求め
$ # その値を EPSG:3395 World Mercator （楕円体のメルカトル）とみなし、緯度経度に再投影し、
$ # 得られた緯度経度の座標を、 EPSG:3857 の地図上で示すため再投影
$ echo 135 34 | \
cs2cs -f %.6f +init=EPSG:4326 +to +init=EPSG:3857 | \
cs2cs -f %.12f +init=EPSG:3395 +to +init=EPSG:4326 | \
cs2cs -f %.6f +init=EPSG:4326 +to +init=EPSG:3857
15028131.257092 4052805.383709 0.000000
最後の例は、2つめと3つめの cs2cs は EPSG:4326 をかまさずに直接 +init=EPSG:3395 +to +init=EPSG:3857 してもよいのですが、わかりやすく。
緯度にして0.1786度（10.7分）ほどの差。
$ # 入力の緯度経度座標と計算結果の緯度経度座標を COORD.dat に保存する
$ echo 135 34 | tee COORD.dat | \
cs2cs -f %.6f +init=EPSG:4326 +to +init=EPSG:3857 | \
cs2cs -f %.12f +init=EPSG:3395 +to +init=EPSG:4326 | \
tee -a COORD.dat
135.000000000001 34.178573658865 0.000000000000
$ # 与える始終点の座標の形式を「緯度1 経度1 緯度2 経度2」にして
$ # invgeod (geod -I) で方位と測地線距離を計算
$ cat COORD.dat | tr "\n" " " | \
awk '{printf("%s %s %s %s", $2, $1, $4, $3)}' | \
geod -I -F %.6f +ellps=WGS84
0d -180d 19808.104553
$ rm COORD.dat
invgeod (geod -I) で計算した楕円体上の測地線距離は19,808.105 m とわかります。
数式による計算
メルカトル図法の投影式
緯度 $\varphi$ の地点を、メルカトル図法で投影したときの $y$ 座標は
y = a \left[ \tanh^{-1} \sin \varphi - e \tanh^{-1} \left( e \sin \varphi \right) \right]
になります。
ただし、 $a$ は楕円体の赤道半径、 $e$ は楕円体の第一離心率です。もし、赤道半径＝極半径とする真球モデルの場合、離心率 $e = 0$ なので第2項が $0$ となります。
緯度より疑似メルカトル
緯度が北緯34度の地点の場合、 EPSG:3857 WGS84 / Pseudo Mercator では、
\begin{align}
a &= 6,378,137 \ \mathrm{m} \\
e &= 0 \\
\varphi &= 34.0^{\circ} = 0.593411945678 \\
\\
y &= a \left[ \tanh^{-1} \sin \varphi - e \tanh^{-1} \left( e \sin \varphi \right) \right] \\
&= a \tanh^{-1} \sin \varphi \\
&= 4,028,802.026 \ \mathrm{m}
\end{align}
と計算されます。今回の場合、（仮に北緯34度に点があった場合）座標値として Shapefile ファイルに格納される値は 4028802.026 となります。
楕円体メルカトルから緯度
この $4,028,802.026$ m という値が、 prj ファイルの不備などで EPSG:3395 WGS84 / World Mercator として認識された場合、 WGS84 の緯度になおすとどのような位置になるか計算すると、 $34.178573659^{\circ}$ となります。
数式は前述の式の逆変換ですが、楕円体の場合、ニュートン法などで繰り返し計算を行うか、下記式のような展開した式で近似計算を行います。式は EPSG Guidance Note #7-2 より。
\begin{align}
a &= 6,378,137 \ \mathrm{m} \\
e &= 0.0818191908426
\end{align}
\begin{align}
\varphi \simeq \chi + \left( \frac{1}{2} e^2 + \frac{5}{24} e^4 + \frac{1}{12} e^6 + \frac{13}{360} e^8 \right) \sin \left( 2 \chi \right) \\
+ \left( \frac{7}{48} e^4 + \frac{29}{240} e^6 + \frac{811}{11520} e^8 \right) \sin \left( 4 \chi \right) \\
+ \left( \frac{7}{120} e^6 + \frac{81}{1120} e^8 \right) \sin \left( 6 \chi \right) \\
+ \left( \frac{4279}{161280} e^8 \right) \sin \left( 8 \chi \right) \\
\end{align}
ただし $\chi$ と $t$ は
\begin{align}
\chi &= \frac{\pi}{2} - 2 \tan^{-1} t \\
t &= \exp \left( - \frac{y}{a} \right)
\end{align}
緯度より疑似メルカトル
EPSG:3395 WGS84 / World Mercator として読み取ると、 $34.178573659^{\circ}$ に位置すると（誤って）解釈されたことになります。この解釈した結果を、 EPSG:3857 WGS84 / Pseudo Mercator の地図上でどこに位置するか計算するには、再度先程の式を使い
\begin{align}
a &= 6,378,137 \ \mathrm{m} \\
e &= 0 \\
\varphi &= 34.178573659^{\circ} = 0.596528643985 \\
\\
y &= a \left[ \tanh^{-1} \sin \varphi - e \tanh^{-1} \left( e \sin \varphi \right) \right] \\
&= a \tanh^{-1} \sin \varphi \\
&= 4,052,805.384 \ \mathrm{m}
\end{align}
となります。
距離
経度が異なる場合、経緯度を用いた2地点間の測地線長、方位角を求める計算にあるとおり、非常に複雑な計算が必要となります。しかし2点の経度が同一であるため、赤道から緯度 $\varphi_1$ までの子午線弧長 $S_{\varphi_1}$ と赤道から緯度 $\varphi_2$ までの $S_{\varphi_2}$ の差を計算することで（若干楽に）求まります。
赤道から緯度 $\varphi$ までの子午線弧長 $S$ は
\begin{align}
S &= \int^\varphi_0 \frac{a (1 - e^2)}{(1 - e^2 \sin^2 \theta)^\frac{3}{2}} \mathrm{d} \theta \\
&= \frac{a}{1 + n} \sum^\infty_{j=0} \left( \prod^j_{k=1} \varepsilon_k \right)^2 \left\{ \varphi + \sum^{2j}_{l=1} \left( \frac{1}{l} - 4l \right) \sin 2l \varphi \prod^l_{m=1} \varepsilon^{(-1)^m}_{j + (-1)^m \lfloor m/2 \rfloor} \right\} \\
& \simeq \frac{a}{1+n} \left\{ \left( 1 + \frac{n^2}{4} + \frac{n^4}{64} \right) \varphi - \frac{3}{2} \left( n - \frac{n^3}{8} \right) \sin 2 \varphi + \frac{15}{16} \left( n^2 - \frac{n^4}{4} \right) \sin 4 \varphi - \frac{35}{48} n^3 \sin 6 \varphi + \frac{315}{512} n^4 \sin 8 \varphi \right\}
\end{align}
1行目は楕円積分を含んだ一般的な式で、2行目は国土地理院の河瀬の式。3行目は河瀬の式を $j=2$ までで打ち切ったものと一致する Helmert の近似式。 $n$ は楕円体の第三離心率です。詳しくは緯度を与えて赤道からの子午線弧長を求める一般的な計算式を参照ください。
$\varphi_1 = 34.0^\circ$ 、 $\varphi_2 = 34.178573659^\circ$ としたとき、赤道から指定緯度までの子午線弧長 $S_{\varphi_1} = 3,763,661.442$ m 、 $S_{\varphi_2} = 3,783,496.547$ m となり、2点間の距離 $| S_{\varphi_1} - S_{\varphi_2} | = 19,808.105$ m となります。
オマケ
なお属性データを格納する dbf ファイルの文字エンコーディングは UTF-8 になっています。
もともと Shift_JIS (CP932) を使うのは日本語 Windows における慣習であり、 Mac や Linux では UTF-8 で作られたりしますので、個人的には UTF-8 にすること自体は特に気になりません。
ですが、実際問題として属性名は10バイト制限であるため、日本語1字が3バイトの UTF-8 では日本語3文字しか使うことができません。一方 Shift_JIS は1字が2バイトなので5文字まで可能です。40%減ですね。
ASCII 文字であれば UTF-8 でも1バイトなので、10文字使えます。ただ国際化の流れのなかで、 UTF-8 で作成するとかえって日本語文字が使いにくくなる現象が生じるのは皮肉的ですね。