OSMデータに海岸線データ(natural=coastline)は含まれるが、適当なノード数の折れ線データが分散して存在する だけである[1]。中には、place=island relation としてポリゴンを形成するものもあるが、そうでないケースもある。
陸地ポリゴンあるいは水域ポリゴンとしてシェイプファイルが提供されているが[2]、テキストファイル形式にするには、例えば、PostgreSQL/PostGISデータベースにインポートして、それをテキスト形式でエクスポートするといった手間がかかる。
そこで、OSMデータから海岸線データを取り出して、それらを繋ぎ合わせて、陸地ポリゴンを構築する。 wayデータだけであれば、OSM Overpass APIで全世界のデータを入手できるが、node id列データであり、 座標値データではない。 狭い範囲であれば、node 座標値データを含めて、海岸線データを取り出せるが、 世界規模では到底無理である。
海岸線データに使われているnode座標値データを後から少しずつ Overpass APIで入手する場合、 OSMデータは時々刻々変動していることへの対応が必要となる。
そこで、今回は、Geofabrikのサイト[3]から europe.osm.pbf、asia.osm.pbf など８つの分割ファイルで 世界全体の OSMデータを入手する。元々は一つのデータを分割したものであるから、８つのファイルは同一時刻の データである。
ダウンロードにも時間がかかり、そこから海岸線データを取り出すのにも時間はかかるが、 将来的には海岸線データだけでなく、国境、国名など、別のデータも取り出す予定であるから、無駄にはならない。
OSM（Google Map方式）では 南緯85.06より南、および北緯85.06より北をカットする[5]。
ポリゴンの面積など、極座標で扱う方が分かりやすいものもあるが、最終的にメルカトル図法で地図を平面に描画する。 ポリゴンの分割はその描画時間短縮のためのものであるから、XY平面座標で考える必要がある。
zoom 0 で言えば、世界地図が１枚のタイルで表される。地図上の左上(西北端)が原点(0, 0)で、右下(東南端)が (1,1) である。 マイOSM地図では、原点を(0, 0)、右下を (230, 230) としている。
南極大陸の海岸線データには、南緯85.06より南の値もあるが、Y座標値としては 230 とみなす。
座標変換は次のメソッドを用いる。
static double Lon2X(double lon, int zoom) {
return (lon + 180.0) / 360.0 * (1<<zoom);
}
static double Lat2Y(double lat, int zoom) {
return (1 - Math.log(Math.tan(lat*Math.PI/180.0) +
1 / Math.cos(lat*Math.PI/180.0)) / Math.PI)/2.0 * (1<<zoom);
}
変換例を下に示す。
System.out.println(Lat2Y(85.05112878, 0)); System.out.println(Lat2Y(-85.051, 0)); System.out.println(Lat2Y(-89.999, 0)); System.out.println(Lat2Y(-90, 0)); -6.227685034332353E-12 0.9999958533607622 2.354016547796194 Infinity
Geofabrik から得たデータでは緯度の最小値は -85.0511287 であった。
海岸線データには wayオブジェクトの id があるが、現時点では使用しないため、海岸線オブジェクト毎に、 node id列を、CSV形式で出力する。最も大きな europe_coastline.csv でも 228MB に過ぎない。 バイナリファイルにすればファイルサイズはかなり小さくなるが、csvファイルの方がチェックしやすいため、 CSVファイルとした。
合計ノード数は約 7,000万である。
*.osm.pbfファイルは巨大である。例えば africa.osm.pbf は 5.6GB であるが、それを解凍すると 117GB となる。 世界全体では約１TB となる。そのため、下に示すように、解凍した xml形式のテキストデータをパイプで受ける。
メモリ使用量は１～２GBであり、３、４プロセスを同時実行できる。実際の処理は wayセクションのみであり、 最初の nodeセクションは空読みに過ぎないが、解凍に時間がかかる。
osmconvert d:/downloads/europe-latest.osm.pbf | java -Dfile.encoding=UTF-8 -Xmx5000m LandPolygon -way europe
62194466,775829176,775829177,775829178,775829176 62202434,775893901,775893911,775893910,775893909,775893908,775893907,775893906,775893905,775893904,775893903,775893902,775893901 62202854,775897086,775897087,775897088,775897089,775897090,775897091,775897086 …
約7,000万ノードの座標値データを得る必要がある。
nide idを座標値に変換するには、HashMap が便利である。 しかし、Java の HashMap はキー、値ともオブジェクト型でなければならない。 int、long、doubleなどプリミティブ型はそのラッパー Ingeger、Long、Double としなければならない。 プログラム上は簡単であるが、メモリ使用量が大きいため、 7,000万もの登録は５GBのメモリでは難しい。
ハッシュを使わない方法も考えられるが、ここでは、自前のハッシュ表を使った。
プリミティブ型の long、int を扱うため、HashMap に比較するとオーバヘッドは微々たるもの である。
ハッシュ関数をシンプルなものとして、ハッシュ表に十分な余裕を持たせた。
static int MAXSIZE = 130000000;
static long[] nodes = new long[MAXSIZE]; // node id
static int num_nodes;
static int[] lons;
static int[] lats;
static int hash(long id) {
return (int)(Math.abs(id * 137) % MAXSIZE);
}
static int put(long id, int lon, int lat) {
int n, h = hash(id);
for (n = 0; n < MAXSIZE; n++) {
int ix = (h + n) % MAXSIZE;
if (nodes[ix] == 0) {
nodes[ix] = id;
lons[ix] = lon;
lats[ix] = lat;
num_nodes++;
if (num_nodes % 1000000 == 0) {
System.out.printf("%d万\r", num_nodes / 10000);
}
return ix; // 新規登録
} else if (nodes[ix] == id) {
return ix; // 登録済み
}
}
System.out.println("Error: Hash Table Full");
return -1;
}
static int get(long id) {
int n, h = hash(id);
for (n = 0; n < MAXSIZE; n++) {
int ix = (h + n) % MAXSIZE;
if (nodes[ix] == 0) {
return -1; // 登録なし
} else if (nodes[ix] == id) {
return ix; // 登録済み
}
}
return -1; // 存在しない
}
このステップでの出力はバイナリで node id（８バイト)、lon（４バイト）、lat（４バイト） の順とした。node id の昇順のデータとなる。
最も大きい europe_idlonlat.dat のサイズは 325MB となった。
この場合は node セクションだけの処理のため、第１ステップより実行時間は短い。
osmconvert d:/downloads/europe-latest.osm.pbf | java -Dfile.encoding=UTF-8 -Xmx5000m LandPolygon -node europe
static void parseNode(String name) {
lons = new int[MAXSIZE];
lats = new int[MAXSIZE];
String file = "c:/osm/" + name + "_coastline.csv";
try (BufferedReader br = new BufferedReader(new FileReader(file))) {
String line;
while ((line = br.readLine()) != null) {
String[] nids = line.split(",");
for (int n = 1; n < nids.length; n++) {
put(Long.parseLong(nids[n]), 0, 0);
}
}
} catch (IOException e) {
e.printStackTrace();
return;
}
int cntFound = 0;
String pathOut = "c:/osm/" + name + "_idlonlat.dat";
// 標準入力
try (BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
DataOutputStream dos = new DataOutputStream(new FileOutputStream(pathOut))) {
String line;
while ((line = br.readLine()) != null) {
int ix = line.indexOf('<');
if (ix < 0) continue;
String s = line.substring(ix); // 行先頭のスペースを削除する
if (s.startsWith("<node ")) { // nodeオブジェクト
String sId = getAttr(" id=", line);
long idNode = Long.parseLong(sId);
if (get(idNode) > 0) { // 海岸線ノード
cntFound++;
String sLon = getAttr(" lon=", line);
String sLat = getAttr(" lat=", line);
int lon = (int)(Double.parseDouble(sLon)*E7);
int lat = (int)(Double.parseDouble(sLat)*E7);
dos.writeLong(idNode);
dos.writeInt(lon);
dos.writeInt(lat);
}
} else if (s.startsWith("<way ")) { // wayオブジェクト
break;
}
}
} catch (Exception e) {
e.printStackTrace();
}
}
第３ステップ（最終処理）の処理時間は短いが、ここがこのプログラムの主題である。 とりわけ、ポリゴン分割は数年来の願望であったが、ようやく実現にこぎつけた。
wayデータの連結は relation処理で経験済みのため、問題はなかった。
海岸線wayデータは先頭node の id をキー、wayデータを値とする HashMap で管理する。 この場合のエントリ数は少ないため、問題はない。
新しい海岸線 way データについては、HashMapのエントリーの前後につながるものがないかを調べる。つながるものがあれば、繋いだ上で、 エントリーを更新する。前と後ろにつながるものがあれば、エントリーは一つ減少する。また、way が閉ループになれば、それは陸地ポリゴンリストに移す。
最終的に閉ループにならないものがあれば、OSMデータそのものにエラーがあることになり、そのときは、 陸地ポリゴンの更新をやめ、以前のデータを使う。自分で誤りを正せる場合にはOSMデータを修正し、そうでないときは誰かかが修正してくれるのを待つ。
今回ダウンロードしたデータでは問題となるエラーはなかった。
自分として初めての経験はポリゴンの分割である。海岸線は入り組んでいるため、均等に都合よくは分割できない。 水平線とか垂直線を引いてみて、東西または南北で海岸線と一度しか交わらない線をポリゴンの分割線とする。
最初に水平分割を終えてから、次に垂直分割を行う。または、その逆でもよい。若干、比較的大きなポリゴンが残るかも知れないが、 それはよしとする。
最終処理の流れは次のようになる。
最初に node の座標値データを読み込み、自前のハッシュ表(配列 nodes, lons, lats)にセットする。
次に、海岸線データ(node id列)を順に読み込み、海岸線の繋ぎ合わせを行う。小さな島などは、繋ぎ合わせは必要なく、 そのまま、ポリゴンとなっている。
小さなポリゴンはそのまま一つのポリゴンレコードとして出力する。今回は、記事[4]のフォーマットで出力する。
繋ぎ合わせによってポリゴンが完成した場合も同じように、小さなポリゴンはそのまま、 大きなポリゴンは分割して、出力する。
閉ループが完成していない途中段階の海岸線データは先頭ノードの id と末尾ノードの id をキーとしたものが HashMap に置かれる。 前に繋がるケースと後ろに繋がるケースがあるため、二つの登録が必要である。 最後に閉ループが完成するときは、前と後ろ同時に繋がる。その前の段階では、前のデータは削除して、繋ぎ合わせた新しいデータを登録する。
// 最終処理
static void generateLandPolygon() throws Exception {
dos = new DataOutputStream(new BufferedOutputStream(new FileOutputStream(pathLands())));
num_nodes = 0;
lons = new int[MAXSIZE];
lats = new int[MAXSIZE];
// ノード座標値データを nodes, lons, lats にセット
for (String area : areas) {
loadIdLonLat(area);
}
System.out.printf("ノード座標値読み込み完了 num_node=%d\n", num_nodes);
for (String area : areas) {
procCoastline(area);
}
unclosedCoastlines();
}
[プログラム・ソース]南極大陸やオーストラリアなどは細かく分割できるが、ユーラシア大陸などは海岸線が入り組んでいるため、分割しにくい。
ただし、世界地図は低ズームでの描画となるため、分割ニーズは日本地図である。 描画時間は計測の上で、必要に応じて、分割をより細かくする。
以前からの分割プログラムを使っており、LandPolygon.java に含まれている。 現時点では特に問題はない。不具合が見つかったときに見直しを行う。