トップ地図アプリMap4 > 陸地の描画

陸地の描画

海岸線のレンダリング[2021.12.13]
OSM海岸線データ
LandParser.java

はじめに

陸地も OSMデータで描画する。しかし、広域森林や水域などとは描画方法が少し異なる。

OSMの海岸線Wayデータを繋ぎ合わせると、巨大なポリゴンとなる。 小さな島程度であれば、森林や湖と同じように描画できる。しかし、大陸は余りにも巨大であり、 日本地図に限っても、本州、北海道、九州などは巨大なポリゴンであるから、 森林や湖などと同じように、描画するのは無理である。

低ズームの場合は分割は要らないが、高ズームでは、 強大なポリゴンをプログラムで無数の小さいポリゴンに分割して、この小さくしたポリゴンを描画する。

陸地ポリゴンデータのダウンロード

サイト[1]からland-polygons-split-4326.zipおよびsimplified-land-polygons-complete-3857.zipをダウンロードした。

座標系は 4326系にするが、低ズーム用はこれが提供されていないため 3857系とした。いずれ、PostgreSQL/PostGIS で 4326系で CSV出力する。

QGISに読込む

低ズーム用ファイルを QGIS に読込み、異常がないことを確かめた。QGISで CSVファイルにエクスポートできるはずだが、 今回は、見送った。

PostgreSQL/PostGISに読込み、CSVファイルに出力する

データベース osm を作り、create extension postgis を実行して置く。

shp2pgsql -I -s 3857 d:/downloads/simplified-land-polygons-complete-3857/simplified_land_polygons.shp land_polygons > imp.sql
psql -h localhost -p 5432 -d osm -U postgres -f imp.sql

データベース osm に テーブル land_polygons が作られる。

COPY (select gid, ST_AsText(ST_Transform(geom,4326)) AS way from land_polygons) 
TO 'c:/map/simplified_land_polygons.tsv' WITH CSV DELIMITER E'\t';

高ズーム用は以下のようにする。

shp2pgsql -I -s 4326 d:/downloads/land-polygons-split-4326/land_polygons.shp land_polygons_split > imp2.sql
psql -h localhost -p 5432 -d osm -U postgres -f imp2.sql
COPY (select gid, ST_AsText(ST_Transform(geom,4326)) AS way from land_polygons_split) 
TO 'c:/map/land_polygons_split.tsv' WITH CSV DELIMITER E'\t';

TSVファイルをOSMバイナリレコードファイルに変換する

GISにおけるMultipolygonは一般に複数の outer polygon があり、それぞれの outer polygon が複数の inner polygon を持つものを言う。上で得られた TSVファイルは Multipolygon 形式であるが、以前は、inner polygon を持たない outer polygon が一つ だけのシンプルな Polygon であった。

しかし、2024年2月1日にダウンロードしたデータには inner polygon が含まれているようである。

標準OSM地図のレンダリングは osm2pgsql で PostgreSQL/PostGISデータベースに OSMデータを取り込むが、この場合の マルチポリゴンは outer polygon は一つである。複数の outer polygon がある場合、複数のレコードとする。 Map4 もこれに従っている。multipolygon の場合、先頭が outer polygon で、その後ろに複数の inner polygon が続く。 それぞれのノード数は先頭にまとめて置かれる。

陸地ポリゴンは実態に合わせて polygon データとする。{key,val} は {'natural','land'} とする。 数年前に作ったプログラムを参考にして、変換プログラムを作る。

TSVファイルの各行は次のようなものである。座標値は経度、緯度の順である。

173	MULTIPOLYGON(((17.6338182 59.0056128,17.6332013 59.0057863,17.6326863 59.0059161,17.632712 59.0061304)))

先頭の ID は当面、無視する。必要に応じて、osm_id と同じ扱いとする。

head(4),  [num_nds(4)], {lon,lat}*, {key,val}*  ... point/line/polygon 
head(4), {num_nds(4)}*, {lon,lat}*, {key,val}*  ... multipolygon
head:
  第0,1bit(0x03) 0: point、   1: line、    2: polygon、 3: multipolygon
  第2bit(0x04)   0: 差分座標、1: 絶対座標
  下位3バイト   レコード長

outer polygon は一つであるが、inner polygon を含むものがあるとして、プログラムを作る。

inner が二つの場合 MULTIPOLYGON(((outer),(inner 1),(inner 2))) となり、outer と inner、inner 同士は "),(" で区切られる。

今回のデータには低ズーム用で1つ(inner 1)、高ズーム用で3つのマルチポリゴン(inner 3, 1, 2)があることが分かった。

低ズームでのマルチポリゴンはカスピ海のようだ。高ズームではカスピ海は分割により、inner polygon ではなくなっていた。

c:\map>java -Dfile.encoding=UTF-8 -Xmx5g -classpath ./class OSMUtil -land high
215342: 4
449441: 2
744401: 3

c:\map>java -Dfile.encoding=UTF-8 -Xmx5g -classpath ./class OSMUtil -land low
63912: 2
LandPolygon.java

メッシュ分割

ポリゴンが大きいため、一般のOSMバイナリレコードよりも小さいzoom値で分割する。 取りあえず、低ズーム用は zoom 1、高ズーム用は zoom 5 と zoom 1 で分割する。

java -Dfile.encoding=UTF-8 -Xmx5g -classpath ./class OSMUtil -devide 1 1 land_polygons_low
java -Dfile.encoding=UTF-8 -Xmx5g -classpath ./class OSMUtil -devide 5 1 land_polygons_high

可変長バイトコード化

c:\map>java -Dfile.encoding=UTF-8 -Xmx5g -classpath ./class OSMUtil -bytecode lands low 1
29916KB/70981KB (42%)
Idx合計=4003KB
レコード数=255781 圧縮なし=284B 可変長=119B 管理=16B
実行時間: 0.05分

c:\map>java -Dfile.encoding=UTF-8 -Xmx5g -classpath ./class OSMUtil -bytecode lands split 1
2984KB/5637KB (52%)
Idx合計=81KB
レコード数=4163 圧縮なし=1386B 可変長=734B 管理=20B
実行時間: 0.01分

c:\map>java -Dfile.encoding=UTF-8 -Xmx5g -classpath ./class OSMUtil -bytecode lands split 5
1137363KB/2269387KB (50%)
Idx合計=48773KB
レコード数=3013466 圧縮なし=771B 可変長=386B 管理=16B
実行時間: 0.97分

zoom 5分割では、日本地図全体は 5/28/12(24.1MB)、5/28/11(4.62MB)、 5/27/12(17.4MB)、5/27/13(24.4MB) となる。 可変長バイトコードでは 5/28/12(13.1MB)、5/28/11(2.50MB)、 5/27/12(9.25MB)、5/27/13(13.0MB) となった。

陸地ポリゴンについては、単独でレンダリングして、結果をストレージに保存しておく。 以後は最初にストレージを調べて、あればこの画像データを使用する。

Androidタブレットではタイル画像を 256x256 とする場合と 384x384 とするケースがある。 画面サイズが小さいため、パソコンと同じ 256x256 では文字が小さすぎるため、 通常は 1.5倍した 384x384 で運用している。 これらのタイルは別のものとして扱う。

ディレクトリは lands256、lands384 とする。スマホでは lands768 となる。

Androidタブレットは主にデバッグ用であり、このために、プログラムを複雑にするのは避ける。

陸地ポリゴンタイルは陸地だけか海だけの単色になることが多い。 ストレージの無駄が問題になるようならば、データとして一括管理する。 zoom 12 で単色ならばその子供の zoom 13 の4タイルも同じ単色である。 zoom 14、zoom 15、... についても同じであるから、管理データをコンパクトにできる。

可変長バイトコードによる陸地ポリゴンの描画を始める

固定長(4バイト)コードによる描画時間

zoom 7 で中心タイルを 7_113_49 とした。タイルのレンダリング時間は約0.2秒、アプリのメモリ使用量は75MB である。

 render ixOsm=3 time=226mS
 tempText1 :72MB
 render ixOsm=24 time=320mS
 tempText1 :73MB
 render ixOsm=4 time=247mS
 render ixOsm=10 time=204mS
 render ixOsm=26 time=225mS
 tempText1 :75MB
 render ixOsm=7 time=206mS
 render ixOsm=15 time=214mS
 render ixOsm=2 time=194mS
 render ixOsm=2 time=229mS
 tempText1 :7_113_49

可変長バイトコードに変更する

変更前(固定4バイト整数)をMap416、変更後(可変長バイトコード)を Map417 とする。

ファイルは Map4/byte 下に置いた。フォルダは将来変更するかも知れない。

Block.get() の Map.DIR を Map.DIRBYTE に変更する。

block.vals = readBytes(block.path); を readAll(block, block.path); に変更する。

block.idxs の生成はなくなる。パソコンで作成したものを readAll で読み込んでいる。

idxs の中身が変更されているため、空間検索 Block.getOSMs も修正した。 こちらは修正が多いので、バグが入り込んだ。 これを好機として、プログラムを見直す。 可変長バイトだけでなく、インタフェースを合わせた固定長整数読み込みメソッドを追加する。 プログラムが少し分かりやすくなる。

type 0 では width, height を省略した。しかし、陸地ポリゴンでは出てこないはず。 まず、管理データのスキャンに問題がある。

先頭の10レコードをチェックした。type = 0 となるのは、データの誤りかも知れない。 LandParserでは type は 2 か 3 になっている。

ByteCoder の次の行が間違い。これでは type はシフトされず、record_length に加算されてしまう。

int head4 = (type | (fShort ? 0 : 0x10)<<24) | record_length;

次のように修正して、データを作り直す。

int head4 = ((type | (fShort ? 0 : 0x10)) << 24) | record_length;

当面シングルスレッドでテストする。

稀にタグエラーがある。

空間検索でレコードが全く抽出されない。低ズームでは精度の変更に未対応である。 高ズームでは東京を中心にすると次のようになった。先頭の固定整数の 経度、緯度からおかしい。タグでもエラーがある。

 render ixOsm=1 time=12mS
 41194: off=165892, head=12012649, type=2, length=75337, x0=1389993858
 -871843827, -1811564029
 parseTag: error ikey=-20178
 -1.5864358E7, NaN

このファイルのバイナリレコード部の先頭を表示してみた。

最初のレコードでは lon、lat は 50CC022D、15284EAC である。これは妥当な値である。 10進で -871843827, -1811564029 となるような値ではない。

しかし、空間検索で取り出したのは先頭レコードではない。

head(1), length,             [num_nds],  {lon,lat}*, {key,val}*
2 18 4  50 CC 2 2D  15 28 4E AC  51 90 7 E0 8 9B 5 8D 8 F3 1 0 B 1 84
2 1F 6 50 C5 34 57 15 2A 5D CC AE 3 B1 2 E9 
2 CB 2 9B 4 4A 24 DE 2 B4 3 D6 1 0 B 1 84 

idxs のオフセットと vals のオフセットがある。ここにミスがあった。 修正で、一応、陸地ポリゴンが描画されるようになった。しかし、まだ、バグが残っている。

低ズームも部分的は描画されるようになった。

以前も似たようなエラーを経験しているが、原因は記憶にない。

int[] vals と byte[] vals のインデックスがごっちゃになっていた。

低ズームで日本地図は描画されるようになったが、大陸がほとんど描画されない。

低ズームで描画されないのは概ねユーラシア大陸のみのようである。 type = 3 で num_polys = 0 となっているのがエラー原因。 ByteCoder.java での設定漏れであった。

これで、ユーラシア大陸も描画されたが、inner polygonの描画が少しおかしい。

マルチポリゴンの復号化にバグがあった。この修正で、世界地図が正しく描画された。

かなりの変更であったが、これで概ね完了した。低ズームの奄美大島だけ描画が欠けている。 タイル境界線で切れているので、多分、描画関連に問題が残っているのであろう。

リファレンス

[1] Data Derived from OpenStreetMap for Download
[2] [C# GIS]SharpMapでシェープファイルの情報を読み取ってみる
[3] QGIS を使ってシェープファイル(Shapefile)を CSV に変換する方法
[4] Daylight Coastlines