トップMy OpenStreetMap > 陸地ポリゴン
前ページ
OSMデータから陸地ポリゴン構築
OSM海岸線パーサー(by SAX)

陸地ポリゴン

はじめに

OpenStreetMapの陸地のレンダリングは一般にはOSMデータから作り出されたシェイプファイル[3]を使う。

自作OSM地図では、現在は Mapnik を使わないため、シェイプファイルのままでは使えない。 これまでは、このシェイプファイルをPostgreSQL/PostGISデータベースに一旦インポートして、 OpenGIS の WKT(Well Known Text)形式でエクスポートしたものを使ってきた。

手間はかかるが、プログラム的には楽である。

今回は、より深く OpenStreetMap を理解するために、 OSMデータから直接陸地ポリゴンデータを取り出すことを試みた。

当面、世界地図は陸地の描画だけであるが、いずれ、国境、国名、大都市名なども描画したい。 Overpass APIには制約があるため、OSMデータから抽出する方がよい。 このためにも、OSMデータの扱いに慣れておきたい。

europe-latest.osm.pbf などから海岸線データを得る

最初に、nodeセクションをスキップして wayセクションから 海岸線データを得る。europe と asia にまたがる海岸線データは同じものが両者にある。

europe の海岸線データだけでは、asia にまたがる場合は単独で繋ぎ合わせても閉ループとはならない。

繋ぐ合わせは行わず、wayオブジェクト単位の海岸線データをファイル出力しておく。 xml形式であっても、世界全体で 2GB弱の大きさであるから、CSVファイルでいいだろう。 行先頭を way id として、その後を node id 列とすればよい。

osmconvert で解凍した xml(テキスト)をパイプで受ける。 全体のパフォーマンスは osmconvert の解凍時間に支配される。日本地図データでは 10分弱である。

nodeセクションは空読みでいいが、解凍に時間がかかる。

世界全体では合計3~4時間程度かかった。同時に2、3の並列実行とした。

osmconvert d:/downloads/asia-latest.osm.pbf | java -Dfile.encoding=UTF-8 -Xmx5000m LandPolygon asia
    // xmlファイルは標準入力(パイプ)から 
    static void parseWay(String name) {
        char quote = '"';
        String pathout = "c:/osm/" + name + "_coastline.csv";
        List<Long> listNodes = new ArrayList<>();
        long idWay = 0;

        try (BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
                FileWriter fw = new FileWriter(new File(pathout)) ) {
            String line;
            boolean fCoastline = false;
            boolean fWaySection = false;

            while ((line = br.readLine()) != null) {
                int ix = line.indexOf('<');
                if (ix < 0) continue;
                String s = line.substring(ix);  // 行先頭のスペースを削除する
                if (fWaySection && s.startsWith("<nd ")) {
                    int bgn = 9;
                    int end = s.indexOf("\"", bgn);
                    String ref = s.substring(bgn, end);
                    long idRef = Long.parseLong(ref);
                    listNodes.add(idRef);
                } else if (fWaySection && s.startsWith("<tag ")) {
                    int kbgn = 8;
                    int kend = s.indexOf(quote, kbgn);
                    String key = s.substring(kbgn, kend);
                    int vbgn = kend + 5;
                    int vend = s.indexOf(quote, vbgn);
                    String val = s.substring(vbgn, vend);
                    if (key.equals("natural") && val.equals("coastline")) {
                        fCoastline = true;
                    }
                } else if (s.startsWith("<way ")) { // wayオブジェクト始まり
                    if (!fWaySection) fWaySection = true;
                    listNodes.clear();
                    fCoastline = false;
                    int bgnId = 9;
                    int endId = s.indexOf(quote, bgnId);
                    String strId = s.substring(bgnId, endId);
                    idWay = Long.parseLong(strId);
                } else if (fWaySection && s.startsWith("</way")) {  // wayオブジェクト終わり
                    if (fCoastline) {
                        fw.write("" + idWay);
                        for (int i = 0; i < listNodes.size(); i++) {
                            fw.write("," + listNodes.get(i));
                        }
                        fw.write("\n");
                    }
                }
            }
        } catch (Exception e) {
            e.printStackTrace();
        }
    }

海岸線wayから参照される node の id, lon, lat を出力する

node の id をキー、lon, lat を値として HashMap で管理するのはノートパソコンでは メモリが足りない。記事[7] に示したハッシュ法を採用した。

LandPolygon は nodeセクションが終わると、実行を終了するため、 前段の osmconvert は write error で終了するが、問題はない。

c:\gisa>osmconvert  d:/downloads/antarctica-latest.osm.pbf | java -Dfile.encoding=UTF-8 -Xmx5000m LandPolygon -node antarctica
実行時間=0.6分
osmconvert Error: write error.
    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;
        }

        char quote = '"';
        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 && ixNode < num_nodes) {
                int ix = line.indexOf('<');
                if (ix < 0) continue;
                String s = line.substring(ix);  // 行先頭のスペースを削除する
                if (s.startsWith("<node ")) {   // nodeオブジェクト
                    int bgnId = 10;  // <node id="31252942" lat="35.6528795" lon="139.7573991" ... >
                    int endId = s.indexOf(quote, bgnId);
                    String strId = s.substring(bgnId, endId);
                    long idNode = Long.parseLong(strId);

                    if (get(idNode) > 0) {  // 海岸線ノード
                        cntFound++;
                        int bgnLon = s.indexOf("lon=", endId);
                        int endLon = s.indexOf(quote, bgnLon+5);
                        String strLon = s.substring(bgnLon+5, endLon);
                        int bgnLat = s.indexOf("lat=", endId);
                        int endLat = s.indexOf(quote, bgnLat+5);
                        String strLat = s.substring(bgnLat+5, endLat);
                        int lon = (int)(Double.parseDouble(strLon)*E7);
                        int lat = (int)(Double.parseDouble(strLat)*E7);
                        dos.writeLong(idNode);
                        dos.writeInt(lon);
                        dos.writeInt(lat);
                    }
                } else if (s.startsWith("<way ")) {   // wayオブジェクト
                    break;
                }
            }
        } catch (Exception e) {
            e.printStackTrace();
        }
    }

同時に2、3の並列実行で約2時間かかった。

陸地ポリゴン構築

海岸線データ(wayデータ)を繋ぎ合わせ、 ノード番号を座標値に変えて、ポリゴンデータとする。

経度 -180度は 180度と同じであることを考慮して繋ぎ合わせている。

プログラムソースコード: LandPolygon.java

ノード座標値が得られないものが1点あった。除外しても問題はないことを JOSM で確認した。

閉ループにならない海岸線がひとつだけあったが、 ノードは3点だけで、-180度の位置であるから、無視して問題はない。

全体で約73万ポリゴンにもなるのにエラーがわずかこの二つに過ぎないのが驚きである。

procCoastLine antarctica
mapCoastLines エントリ数=0
procCoastLine africa
mapCoastLines エントリ数=66
procCoastLine europe
mapCoastLines エントリ数=74
procCoastLine asia
mapCoastLines エントリ数=4
procCoastLine north-america
Error 座標値が得られない Id=7670000000
mapCoastLines エントリ数=8
procCoastLine central-america
mapCoastLines エントリ数=6
procCoastLine south-america
mapCoastLines エントリ数=4
procCoastLine australia-oceania
mapCoastLines エントリ数=2
3 (-1800000000 675522267) (-1800000000 669296325)

ポリゴン数=729300
実行時間=1.53分

おわりに

結果には2つ問題がある。

一部の描画がおかしくなった。ポリゴンの外側が塗りつぶされているようだ。

第二に、ポリゴンを適当に分割しないと大陸ポリゴンが大きすぎる。 世界地図は間引きをした低ズーム用のみとして、日本地図は間引き無しとする案があるが、 パフォーマンス上は分割が望ましい。

なんとか分割方法を考えたい。


時計回りでないものがあった。これを反転するのを忘れていた。この追加修正で描画エラーはなくなった。

リファレンス

[1] OSM陸地ポリゴンをOSMデータから構築する
[2] OSM海岸線パーサー(by SAX)
[3] Data Derived from OpenStreetMap for Download
[4] 空間オブジェクトの表現: OpenGIS WKT形式
[5] OSM Overpass API
[6] OpenStreetMapのデータ構造
[7] ハッシュ法