トップ＞ OpenStreetMap ＞ OSM陸地ポリゴンシェイプファイルをPostgreSQLにインポートする

OSM陸地ポリゴンシェイプファイルをPostgreSQLにインポートする

１．はじめに

海岸線データ(natural=coastline)は osm.pbf に含まれているが、 osm2pgsql ではデータベースにインポートできない。

シェープファイルが入手できる。これをデータベースにインポートすれば tsv ファイルとして出力できる。

このデータをみると、これまでの osm tsv ファイルとは異なり、 MULTIPOLYGON形式であった。

プログラムで簡単に調べた結果では、実体は outer だけのシンプルな POLYGON であった。inner はない。

海岸線データといっても、本州、九州といった大きな島が一つのポリゴンになっているわけではなかった。適当な大きさのブロックに分割されている。海岸線が全くない内陸部にも矩形データがある。要するに、海岸線データというより、陸地描画用ポリゴンデータと言える。つまり、誤りを自動補正しただけではなく、レンダリングに適したように加工されている。

OSMでは禁止タグであるが、仮想的には natural=land と考えた方が分かりやすい。

陸地ポリゴンとは双対的な水域ポリゴン(water polygon)も提供されている。こちらは、水域(海)を outer polygon として、島々を inner polygon とする OSM ではマルチポリゴンと呼ぶデータとなっている。

２．海岸線データをインポート

データベース osm のテーブル land_polygons に海岸線データをインポートするには、次のようにする。

shp2pgsql -I -s 3857 d:/osmdata/world_boundaries/land_polygons land_polygons > imp.sql
psql -h localhost -p 5432 -d osm -U postgres -f imp.sql

パイプで繋いで１コマンドで実行することもできる。

shp2pgsql -I -s 3857 d:/osmdata/world_boundaries/land_polygons land_polygons | ^
psql -h localhost -p 5432 -d osm -U postgres

psql osm で postgres に入って \d コマンドでテーブル内容を確認すると以下のようになる。

                                      Table "public.land_polygons"
 Column |            Type             | Collation | Nullable |                  Default                   
--------+-----------------------------+-----------+----------+--------------------------------------------
 gid    | integer                     |           | not null | nextval('land_polygons_gid_seq'::regclass)
 fid    | double precision            |           |          | 
 geom   | geometry(MultiPolygon,3857) |           |          | 
Indexes:
    "land_polygons_pkey" PRIMARY KEY, btree (gid)
    "land_polygons_geom_idx" gist (geom)

世界地図データであるが、昨年(2020年)のデータでは約70万行であり、次のようにして TSVファイルに出力すると、約1GBになった。

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

このデータは以下のように、outer polygon が２つ以上あってもよい MULTIPOLYGON形式であるが、分割データの場合、実際は outer がひとつで、 inner はない単純な polygonデータである。

173	172	MULTIPOLYGON(((17.6338182 59.0056128,17.6332013 59.0057863,17.6326863 59.0059161,
17.632712 59.0061304,17.6332914 59.0061481,17.6340467 59.0058299,17.6338182 59.0056128)))

[2022.5.22にdownloadしたファイルに対して実行]

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

shp2pgsql -I -s 3857 c:/osm/land-polygons-split-3857/land_polygons land_polygons > imp.sql
psql -h localhost -p 5432 -d wrd -U postgres -f imp.sql

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

c:\gisa>java -Dfile.encoding=UTF-8 -Xmx5000m Parser world
landParserWorld start
Parser実行時間 = 1.9分

world_polygons.dat は552MB となった。

２．海岸線データのブロック分割

現在、zoom 13以上のレンダリング用のosmデータは zoom 11 で分割している。世界地図対応の海岸線データを zoom 11 で分割するとファイルサイズが膨大になる。このため、海岸線データは zoom 3, 8 で分割している。 zoom 8 での分割ファイルは元のデータのままであるが、 zoom 3 では、許容誤差 500ｍで、間引きを行っている。 zoom 3, 8 のファイルサイズは 16MB, 1.04GB である。

レコードフォーマット

海岸線データのレコード例を下に示す。最初の２項目はレンダリングには使っていない。後続のデータは整数化した経度、緯度を16進数で表したものである。 OSMデータのように差分を可変長コードで表せば、ファイルサイズは 1/3 程度に縮小するであろう。

50695	50694	b920ce13 28f197bb,b920ce13 27fad03a,b4efc298 27fad03a,b4efc298 28f197bb,b920ce13 28f197bb

描画プログラム

現在は、下の海岸線(陸地)描画専用のプログラムを使っている。 (int)(val + 0.5) は val が負の値のときは必ずしも妥当ではないが、特に精度が要求されるのは、中・高ズームであり、日本地図が対象となるため、val(緯度、経度)は正の値であるため、四捨五入に当たる。

現実のOSMデータは仕様上の精度(小数点以下７桁)はないので、さほど気にすべきことではない。単に (int)v でも問題はない。

共通化を試したが、いい結果は得られなかったので、別プログラムのままとする。

using System;
using System.Drawing;
using System.Drawing.Imaging;
using System.Drawing.Drawing2D;
using System.Globalization;
using System.IO;
using System.Linq;
using System.Text;
using System.Collections.Generic;
using System.Collections.Specialized;
using System.Text.RegularExpressions;
using System.Threading;
using System.Diagnostics;


class Coastline {

    Brush brLand  = new SolidBrush(Color.FromArgb(255, 0xf3, 0xf0, 0xea));

    OrderedDictionary caches; 

    public Coastline() {
        caches = new OrderedDictionary();
    }

    public void Render(Graphics g, int zoom, int x, int y) {
        OSM[] osms0 = Read(zoom, x, y);
        OSM[] osms1 = Select(zoom, x, y, osms0);
        foreach (OSM osm in osms1) {
            FillPolygon(g, brLand, osm.points); // 陸地
        }
    }

    OSM[] Read(int zoom, int x, int y) {
        var listOsm = new List<OSM>();
        string path;
        if (zoom < 8) {
            path = "d:/osmdata/w3/" + (x>>(zoom-3)) + "/" + (y>>(zoom-3)) + ".dat";
        } else {
            path = "d:/osmdata/w8/" + (x>>(zoom-8)) + "/" + (y>>(zoom-8)) + ".dat";
        }
        if (caches.Contains(path)) {
            return (OSM[])caches[path];
        }

        if (!File.Exists(path)) {
            return listOsm.ToArray();
        }

        var hex = System.Globalization.NumberStyles.HexNumber;
        IEnumerable<string> lines1 = File.ReadLines(path);
        foreach (string line in lines1) {
            string[] v = line.Split('\t');
            if (v.Length < 3) continue; 
            var rec = new OSM();
            int xmin=int.MaxValue, ymin=int.MaxValue;
            int xmax=int.MinValue, ymax=int.MinValue;
            var list = new List<Pnt>();
            var way = v[2];
            for (int n = 0; n < way.Length; n += 18) {
                var lon = way.Substring(n, 8);
                var lat = way.Substring(n+9, 8);
                int xf = int.Parse(lon, hex);
                int yf = int.Parse(lat, hex);
                list.Add(new Pnt(xf,yf));
                if (xf < xmin) xmin = xf;
                if (xf > xmax) xmax = xf;
                if (yf < ymin) ymin = yf;
                if (yf > ymax) ymax = yf;
            }
            rec.points = list.ToArray();
            rec.min = new Pnt(xmin, ymin);
            rec.max = new Pnt(xmax, ymax);
            listOsm.Add(rec);
        }
        OSM[] osms = listOsm.ToArray();
        if (caches.Count >= 32) {
            caches.RemoveAt(0);     // 最初に登録したものを削除する
            GC.Collect();
        }
        caches[path] = osms;
        return osms;
    }

    OSM[] Select(int zoom, int X, int Y, OSM[] osms) {
        int bx1 = (int)(Zone.X2Lon(X, zoom)*10000000 + 0.5);
        int by1 = (int)(Zone.Y2Lat((Y+1), zoom)*10000000 + 0.5);
        int bx2 = (int)(Zone.X2Lon((X+1), zoom)*10000000 + 0.5);
        int by2 = (int)(Zone.Y2Lat(Y, zoom)*10000000 + 0.5);
        var list = new List<OSM>();
        foreach (var osm in osms) {
            int ex1 = osm.min.x, ey1 = osm.min.y;
            int ex2 = osm.max.x, ey2 = osm.max.y;
            if (bx1 <= ex2 && ex1 <= bx2 && by1 <= ey2 && ey1 <= by2) {
                list.Add(osm);
            }
        } 
        return list.ToArray();
    }

    void FillPolygon(Graphics g, Brush brush, Pnt[] points) {
        Stopwatch sw = new Stopwatch();
        sw.Start();
        if (points.Length == 0) return;
        PointF[] pts = new PointF[points.Length-1];
        for (int n = 0; n < pts.Length; n++) {
            Pnt pnt = points[n];
            pts[n] = new PointF(pnt.xf, pnt.yf);
        }
        sw.Stop();
        //Console.WriteLine("conv: {0}ticks", sw.Elapsed.Ticks);
        sw.Start();
        g.FillPolygon(brush, pts, FillMode.Alternate);
        //Console.WriteLine("draw: {0}ticks {1}ms", 
        //      sw.Elapsed.Ticks, sw.Elapsed.Milliseconds);
    }

}

３．おわりに

陸地描画は森林や農地などポリゴン描画と何ら変わらない。しかし、低ズーム用データをメモリ常駐とする必要がないことと世界地図を対象にすると、zoom 11分割ではファイル数が膨大になることから、日本地図用osmデータと一体化しない方がいい。

ただし、プログラム自体は共通化できる部分は残っている。シンプルな方法で共通化できるならばそうした方がよい。 zoom 10以上では陸地描画は描画時間、メモリ使用量とも全体に占める割合は極めて小さいので、パフォーマンスに与える影響は僅少である。

Ａ．リファレンス

[1] Data Derived from OpenStreetMap for Download

Ｂ．来歴およびノート

2021.4.17