ホーム < ゲームつくろー！ < IKD備忘録

何か色々徒然と
国土地理院の「地球地図」の行政界ベクターデータから都道府県の形を作る

(2019. 2. 21)


　もはやこの情報を必要とする人が何人いるんだという話ですか(^-^;、僕が必要になってしまったので備忘録として書き連ねます。

　勤め先のクラブ活動的なものに「ゲームハッカソン部」とういのがありまして、毎週何らかのゲームを作って発表し合おうという物なのですが、とある週のお題が「地図」になったんです。「地図を使ったゲームかぁ」とあれこれ考えて、まぁ地図を使ったジグソーパズルのような物でも作ろうかなと思い立った訳です。で、イメージしたのが都道府県のジグソーパズルかなと。そこで各都道府県の形を切り出す必要性が出てきたという事です。まぁ簡単に作るなら、どこぞやから白地図を持ってきて、そこから各都道府県の絵をスプライトで作って云々…とやれば良いのですが、せっかくなら正確な形で切り出したいなぁと(^-^;。んで色々と調べると、国土地理院が「地球地図」というデータを公開してくれていて、そこにどうやら各都道府県の境界線を表すベクターデータがある事が分かりました。それを使えばとても正確な都道府県の形が作れるはずだろうと(^-^)。いきさつはそういう感じです。



① 国土地理院の「地球地図」データ

　国土交通省 国土地理院は日本の地形を測定している行政組織です。今の日本のありとあらゆる地形情報をを測定し管理しています。国土地理院が発行する日本地図は、最新のデータを反映した日本の基準となる地図です。この国土地理院がまとめてくれている地形情報の一つが「地球地図日本」です。国土地理院が測定したデータを世界標準のフォーマットにしてくれています。ですから地球地図日本のデータを活用できるようになると、世界中の同じフォーマットの地図情報を扱えるようになる訳です。

　地球地図のデータは様々なカテゴリーに分類されています。地図と一括りに言っても、海岸線、行政区分、高度など様々な項目があります。今回扱いたいのは都道府県の形。それは「行政区分」、地球地図の言葉だと「行政界」というカテゴリーになります。

　地球地図の情報は一般公開されていて、ライセンスの範囲内で誰でも自由に扱う事ができます。地球地図のホームページからダウンロード可能です。ホームページにある利用規約をよく読んで使用する必要がありますが、掻い摘んで言えば、

以上が必須となります。

　今回欲しいのは都道府県の形の情報。それは地球地図日本の行政界カテゴリーにあります。先のリンク先の「第2.1版ベクタ」の「行政界」がそれに該当します。



② 地球地図行政界ベクターデータ

　地球地図の行政界ベクターデータをダウンロードして展開すると、次のようなファイル構成になっています：

　普段見慣れない拡張子ばかりです(^-^;;。これらは「地球地図国際運営委員会の地球地図仕様（Global Map Specifications）」で定められています。Global Map SpecificationsはGithubに情報をまとめてくれています。残念ながら日本語はありませんので、これを参考にするなら英語を一生懸命読む事になります、はい(^-^;;。ちょっと面倒です。

　コンテンツを読む前に、上のファイル名などをキーにGoogle先生に検索してもらうと、なるほど色々情報が出てきます。こちらのサイト様の情報から、.dbfが項目のデータベースでOpenOfficeなどで開く事が可能、.shpが形状情報、.shxがそのインデックス情報である事がわかりました。感謝です(>_<)/。で、.dbfファイルがOpenOfficeで開けるならExcelでも開けるんでしょ～ってそのまま開いてみたら、おーーちゃんと中を見られる！！polbnda_jpn.dbfの中身はこうなっておりました：

　JPNでHokkai DoでSapporo Shiと、なるほどpolbnda_jpn.dbfはどうやら日本のすべての市町村の情報のようです。popは人口みたいですね。ん？でもKushiro Shiが2つある？釧路市の人口は9千万人ではないのは明らかなので（笑）、多分これは釧路市に含まれる行政区分が二つあって、一番上のレコードに人口(180160人)があるので2つめ以降は「無効」という事かなと。

＜調査中…＞

　なるほど、市町村合併時、釧路市の左隣に白糠町さらにその左に音別町があったのですが、白糖町が合併拒否したのに対し音別町は合意したため飛び地になったとの事。つまり同じ名前になっているのは飛び地のためのようです。実際このリストには2914のレコードがありますが、日本の市町村の数は2018年10月時点で1741だそうです（Wikipedia：市町村）。

　都道府県の情報は無いのかなぁと、もう一つのpolbndl_jpn.dbfを開いてみるも…：

　…これは何だろう…。＜調査中＞
　はい、Global Map Specifications 2.1の20ページに意味がありました。これは「Political Boundary Line（行政境界ライン）」で「bst」というのはBoundary Status Type（境界状態タイプ）で、1が「境界確定」、2が「境界未確定」そして3が「論争中」なんだそうです。へぇぇぇ。確かに境界線が確定していない地域が日本にもいくつもありますものね。実に細かいデータがとられているんですねぇ。
　ただ、このデータは都道府県の形ではありませんので今は無視と。んー、この様子だとpolbanda_jpn.***から地図を作らないといけないのかなぁ…。でも市町村に分かれている…？嫌な予感しかしませんが、続けましょう。

　形状情報である.shpファイルはバッキバキのバイナリファイル。このフォーマット情報を知らなければ地図は作れません。これについてWikipediaに結構な情報が上がっていました（Wikipedia：シェープファイル）。中々のボリュームなので②以降でじっくり見ていきます。



② シェープファイル

　シェープファイルは地理情報システム（GIS）でデータを相互運用する時に使われるファイル形式のようです。.shp（シェープ規格）、.shx（シェープインデックス規格）そして.dbf（属性規格）の3つのファイルで地図上の幾何学的な位置とその属性情報を表します（やっと拡張子の情報がハッキリとわかった(T_T)）。シェープファイルというと.shpを指す事もあるようですが、実際は上の3つのファイルが無いと機能しないので注意が必要、との事。



③ シェープ規格（.shp）

　地理データを格納しているシェープファイルデータの本丸でバイナリファイルになっています。最初の100バイトに固定長のファイルヘッダー、その下に任意の数の可変長レコードが並びます。可変長レコードの先頭8バイトにはレコードヘッダーが付いていて、その下にレコード本体が続くようです。レコード本体の最初4バイトにシェープの種類、その下がそのデータとなっています。言葉だとさっぱりなのでww、表にすると以下の通り：

項目 バイト数 意味
ファイルヘッダー 100 シェープ全体の情報（ファイル長、シェイプ種別、地図のXY座標バウンディング領域、高さZの範囲、ユーザ情報Mの範囲）
レコード1 レコードヘッダー 8 レコード1番全体の情報（レコード番号、レコード長(バイト数ではなくて16bitWORD数のようです)）
レコード本体 Shape種別 4 (int32) 点(Point：1）、多角線（Polyline：3）、多角形（Polygon：5）など様々
Shape内容 n Shape種別により決まっている
レコードn

　より詳しくはWikipediaをご参照下さい。

　解析時に注意したいのは、多分歴史的な背景のせいかなと思うのですが、ヘッダー内の要素にビッグエンディアンとリトルエンディアンが混ざっている（！）という点と、サイズが16bitWORDの個数で表現されている所です。昨今のCPUはリトルエンディアンなので、ビッグエンディアンの数値はバイト列を逆に読む必要があります。またサイズを計算するには個数×2にしなければなりません。


④ シェープインデックス規格（.shx）

　シェープ規格は不定長にデータが並ぶので特定のレコードにジャンプするのが大変です。シェープインデックス規格は各レコードまでの絶対位置が固定長で記録されているため、素早くジャンプ出来ます：

項目 バイト数 意味
ファイルヘッダー 100 .shpと同じ内容
レコード1 レコード位置 4 .shpのレコード1の絶対位置。16bitWORD数で記述
レコード長 4 .shpのレコード1の長さ。16bitWORD数で記述
レコードn

ファイルヘッダー下から8バイトな構造体を作ればインデックスでおりゃっとジャンプ出来そうです。ただし、そのレコード（レコード位置及びレコード長）はいずれもビッグエンディアンです。



⑤ で、polbnda_jpn.shpはどうなっている？

　という事で各ファイルの中身について結構分かってきました。shapeの種別によってデータ形式が異なるので、そこはさらに調査が必要ですが、取り敢えずpolbnda_jpn.shpのヘッダー情報等を読み取るプログラムをざざっと作ってみます：

　※ 以下は中途版です。完全なshape.hは本章のサンプルプログラムに公開致します。

shape.h
#ifndef __IKD_SHAPEFILE_SHAPE_H__
#define __IKD_SHAPEFILE_SHAPE_H__

// シェープ規格（.shp）

#include <stdint.h>
#include <map>

namespace ShapeFile {

    enum ShapeType {
    ShapeType_Null = 0,         // 空データ
    ShapeType_Point = 1,        // 1点（XY）
    ShapeType_Polyline = 3,     // 多角線
    ShapeType_Polygon = 5,      // 多角形
    ShapeType_MultiPoint = 8,   // 複数点
    ShapeType_PointZ = 11,      // 3次元点（XYZ、オプションM）
    ShapeType_PolylineZ = 13,   // 3次元多角線
    ShapeType_PolygonZ = 15,    // 3次元多角形
    ShapeType_MultiPointZ = 18, // 3次元複数点
    ShapeType_PointM = 21,      // M付き1点
    ShapeType_PolilineM = 23,   // M付き多角線
    ShapeType_PoligonM = 25,    // M付き多角形
    ShapeType_MultiPointM = 28, // M付き複数点
    ShapeType_MultiPatch = 31,  // 複数パッチ
    ShapeType_Error = 999,      // 不明
};

class ShapeFileHelper {
public:
    // シェープタイプを取得
    static ShapeType getShapeType( int32_t shapeTypeNum ) {
        static std::map< int32_t, ShapeType > shapeTypeMap_g = {
            { (int)ShapeType_Null, ShapeType_Null },
            { (int)ShapeType_Point, ShapeType_Point },
            { (int)ShapeType_Polyline, ShapeType_Polyline },
            { (int)ShapeType_Polygon, ShapeType_Polygon },
            { (int)ShapeType_MultiPoint, ShapeType_MultiPoint },
            { (int)ShapeType_PointZ, ShapeType_PointZ },
            { (int)ShapeType_PolylineZ, ShapeType_PolylineZ },
            { (int)ShapeType_PolygonZ, ShapeType_PolygonZ },
            { (int)ShapeType_MultiPointZ, ShapeType_MultiPointZ },
            { (int)ShapeType_PointM, ShapeType_PointM },
            { (int)ShapeType_PolilineM, ShapeType_PolilineM },
            { (int)ShapeType_PoligonM, ShapeType_PoligonM },
            { (int)ShapeType_MultiPointM, ShapeType_MultiPointM },
            { (int)ShapeType_MultiPatch, ShapeType_MultiPatch },
        };
        if ( shapeTypeMap_g.find( shapeTypeNum ) == shapeTypeMap_g.end() )
            return ShapeType_Error;
        return shapeTypeMap_g[ shapeTypeNum ];
    }

    // エンディアン変換
    static int32_t toLittle( int32_t bigEndian ) {
        return
            ( ( bigEndian & 0xff ) << 24 ) |
            ( ( bigEndian & 0xff00 ) << 8 ) |
            ( ( bigEndian & 0xff0000 ) >> 8 ) |
            ( ( bigEndian & 0xff000000 ) >> 24 );
    }
};


#pragma pack(1)

// シェープレコード
struct ShapeFileRecord {
    int32_t shapeType_; // シェープタイプ
    int32_t data_; // データ先頭

    // データの先頭ポインタを取得
    void* getDataPtr() {
        return (void*)&data_;
    }

    // シェープタイプを取得
    ShapeType getShapeType() {
        return ShapeFileHelper::getShapeType( shapeType_ );
    }
};


// シェープレコードヘッダー
struct ShapeFileRecordHeader {
    int32_t index_big_; // レコードインデックス(1基底）
    int32_t len_big_; // レコード長（16bitWORD個数）
    ShapeFileRecord record_; // レコード先頭

    // インデックスを取得(1基底）
    int32_t getIndex() {
        return ShapeFileHelper::toLittle( index_big_ );
    }

    // レコードサイズを取得
    int32_t getSize() {
        return ShapeFileHelper::toLittle( len_big_ ) * 2;
    }
};

// シェープ規格ファイルヘッダー
struct ShapeFileHeader {
    int32_t id_big_;            // ファイル符号(0x0000270a) Bigendian
    int32_t reserved_big_[ 5 ]; // 未使用領域
    int32_t fileLen_big_;       // ファイル長（16bitWORD個数）
    int32_t version_;           // バージョン番号
    int32_t shapeType_;         // Shape種別
    double minX_;  // 最小境界矩形X座標最小値
    double minY_;  // 最小境界矩形Y座標最小値
    double maxX_;  // 最小境界矩形X座標最大値
    double maxY_;  // 最小境界矩形Y座標最大値
    double minZ_;  // 高さ情報最小値
    double maxZ_;  // 高さ情報最大値
    double minM_;  // ユーザ情報最小値
    double maxM_;  // ユーザ情報最大値
    char recordHeader_; // レコードヘッダー先頭

    // ファイルサイズを取得
    int32_t getFileSize() {
        return ShapeFileHelper::toLittle( fileLen_big_ ) * 2;
    }

    // シェープタイプを取得
    ShapeType getShapeType() {
        return ShapeFileHelper::getShapeType( shapeType_ );
    }

    // レコードヘッダーを取得
    ShapeFileRecordHeader* getRecordHeaderPtr() {
        return (ShapeFileRecordHeader*)&recordHeader_;
    }
};


// シェイプインデックス
struct ShapeIndex {
    int32_t recordPos_big_; // レコード位置
    int32_t recordLen_big_; // レコード長(16bitWORD個数)
};

// シェイプインデックスヘッダー
struct ShapeIndexHeader : ShapeFileHeader {
    // インデックスヘッダーを取得
    ShapeIndex* getShapeIndexPtr() {
        return (ShapeIndex*)&recordHeader_;
    }

    // レコード数を取得
    int32_t getRecordNum() {
        int32_t fileSize = getFileSize();
        return ( fileSize - 100 ) / sizeof( ShapeIndex );
    }
};

#pragma pack()

}  // end namespace

#endif

このshape.hを用いてファイル内の情報を調べます：

main.cpp
#include "stdafx.h"
#include "code/Shape.h"

using namespace ShapeFile;

char* getBuf( const char* fileName ) {
    FILE *f = 0;
    auto err = fopen_s( &f, fileName, "rb" );
    if ( f == 0 )
        return 0;

    fseek( f, 0, SEEK_END );
    auto fileSize = ftell( f );
    fseek( f, 0, SEEK_SET );

    char *buf = new char[ fileSize ];
    fread( buf, fileSize, 1, f );

    fclose( f );

    return buf;
}

int main( int argc, _TCHAR* argv[] )
{
    char *buf = getBuf( "polbnda_jpn.shp" );
    if ( buf == 0 )
        return -1;
    ShapeFileHeader *fileHeader = ( ShapeFileHeader* )buf;
    ShapeType shapeType = fileHeader->getShapeType();
    delete[] buf;

    buf = getBuf( "polbnda_jpn.shx" );
    if ( buf == 0 )
        return -1;
    ShapeIndexHeader *indexHeader = (ShapeIndexHeader*)buf;
    int32_t recordNum = indexHeader->getRecordNum();
    delete[] buf;

    return 0;
}

　上のクラスと実行コードでpolbnda_jpn.shpを調べた結果、ファイルヘッダーに格納されていたShape種別はShapeType_Polygon、つまり「多角形」でした。またレコード数は2914個である事がわかりました。これ、先ほどExcelで開いたpolbnda_jpn.bdf内のレコード数とピッタリ一致です(^-^)。よし。



⑥ 多角形フォーマット

　シェープファイルのShape種別が多角形と分かりました。これはどういうフォーマットであるか調べます。シェープファイルはESRIという団体のオープンフォーマットで、こちらに仕様書(PDF)がありました。これの8ページ目にありました「Polygon」。どうやら次のようになっているようです：

Polygonフォーマット（※リファレンス参照）
Polygon
{
    Double[4] Box       // Bounding Box
    Integer NumParts    // Number of Parts
    Integer NumPoints   // Total Number of Points
    Integer[NumParts] Parts // Index to First Point in Part
    Point[NumPoints] Points // Points for All Parts
}

Boxはポリゴン図形のバウンディングボックス (minX, minY) - (maxX, maxY)です。
NumPartsはポリゴンを形作る個々のパーツ（サブポリゴン）の数です。地図の形状は複雑で、穴あきなどあるため単純なポリゴンだけでは表現しきれない為パーツに分けられています。
NumPointsは全ポリゴンの頂点の数です。
Partsは各パーツの頂点インデックスの先頭番号です。それぞれのパーツはその先頭番号から次のパーツの一つ前の頂点までで形作られます。
Pointsはすべての頂点座標でdouble型×2です (x,y)

　パーツが複数ある場合、考えられる事としては「飛び地」と「穴あき」です。両者の違いはポリゴン面の法線方向の違いです。リファレンスを見ると、多分ですが右回りになっている場合が有効な地面、左回りで囲まれる部分は穴あき（湖など）と表現しているように見受けられます。つまり飛び地は右回り、穴あきは左回りです。概念としては分かりますが、実装となるとこの辺りは面倒な事になりそうです。ま、取り敢えず多角形のフォーマットがわかりました。この不定長な構造体がシェープ規格のレコード本体のShape内容の所にドカッとある訳です。ここから試しにポリゴンを描画してみる事にします。



⑦ polbnda_jpn.shpの最初のポリゴンを描画してみる

　作成したshape.hを用いてpolbnda_jpn.shpのレコード1番のポリゴン情報を抜き出し、描画してみました（実際のコードはサンプルをご覧下さい）：

左側が点を連結したマップです。レコード1番はpolbnda_jpn.dbfを見ると「札幌市」。Google Mapの札幌市（縮尺縦横比調整）と見比べると、おーーーーピッタリ！札幌市ってこういう形していたんですねぇ。札幌市のデータにはパーツが一つしかありませんでした。上のポリゴンのみという事です。

　そういえば、polbnda_jpn.dbfには釧路市のレコードが2つ（6番と7番）ありました。釧路市を描画してみます：

　なーるほど、確かに釧路市は2つに分かれています。飛び地はそれぞれが行政区分として分離しているんですね。

　という事で、すべての市町村とその飛び地を描画したのが以下の地図です：

この数だけの行政区分があるのだから凄い事です。もしこの単位でジグソーパズルにしたら…地獄だけどやってみたいかもw。製品化したい企業がありましたらデータ作成お手伝いしますよw



⑧ 都道府県を得るには？

　さて、何とか各市町村や飛び地、島などを配置し日本の形を作り出す事が出来ました。しかし今欲しいのは都道府県の形です。個々の行政区分がどの都道府県に属するかは分かっています。よって、例えば北海道に属する行政区分を全部くっつけてしまえば北海道の形が浮き彫りになります。しかしそれだと海岸線を線で描く事が出来ません。海岸線が描けないとポリゴン化できないのです。

　海岸線と県境のみを抽出するには、それ以外の内陸境界線を消してしまえば良いんです。つまり海岸線と内陸境界線を区別できれば良いと。両者の最大の違いは「共有者がいるかどうか」です。海岸線は一つの行政区分にしか属しませんが、隣の町との境界線は双方が共有しています。ポイントはこれです。

　地球地図日本の隣り合う町の辺、もっと言えば辺を構成する頂点は重なっているはずです。それが重なっていないと境界線がクロスしたり穴あきになるなどおかしな事が起こりますから。でも本当にそうなっているか一応テストしてみます。北海道旭川市と隣接する深川市の頂点に〇をつけて描画するとこうなります：

　予想通り双方の境界線の頂点は寸分たがわずぴったりと一致していました。頂点間を結ぶ辺も勿論一緒です。これを抽出して取り除けば両市は一つになります。

　共有境界線を取り除けたとして、次にしなければならない事は海岸線を繋げる作業です。上の両市の赤丸と青丸が海岸線だとして、双方はデータ上は接点で分離しています。これを繋げて一つの多角形にしないとポリゴン化出来ません。共有線の削除と海岸線の連結、それを各都道府県に所属する数十から週百の市町村で行わなければなりません。さて…どうしたものか…。



⑨ セル空間分割で共有境界線をチェック

　上の問題、結局は一つ一つの辺について調べるしか無さそうなのですが、問題はそのチェック数。例えば愛知県の場合、全部で70の行政区画があり、辺の数は5135本。ある辺が別の辺と共有であるかどうかを総当たりで調べると 5134 * 5135 / 2 = 1318万通り。まぁ昨今のPCなら力技で比較すればそれ程時間が掛かるものでは無いですが、それはあんまりにも無策かなと。そこで、最適解では無いかもしれませんが「セル空間分割」を使う事にしました。

　セル空間分割は対象とする空間をセル状に分割し、各セルに対象物を登録する事で、調べる対象物を絞り込む方法です。イメージは以下のような感じです：

　各セル内に所属する辺をそのセルに紐づいたリストに登録します。こうするとある辺と共有する辺を探す時に、自分が所属しているセルの中だけを見れば事足ります。これで検索する組み合わせ数が劇的に下がります。上の例は見やすくするため10×10のセルにしていますが、実際はもっと細かくセルを分割して効率を上げています。

　この方法で全ての辺についてそれと共有する辺をチェックし、共有されていない辺だけを残していくと…：

愛知県！綺麗に愛知県(^-^)。物の見事に内側にある辺だけが無くなり、アウトラインだけが残りました。ただ、見た目繋がっているように見えるこれらの辺はまだバラっバラの状態です。ここからこの愛知県をN角形のポリゴンにする必要があります。



⑩ セル空間分割で隣の辺を見つける

　見た目繋がっているように見える上の辺の集合。実際は共通境界線を消す作業でバラバラになってしまいます。これを再度つなぎ直すために、先程と同じセル空間分割を利用します。自分の隣にある辺を得るために自分が所属しているセルに登録されている辺だけをチェックすれば効率的です。

　線同士を繋げる方法は何でも良いのですが、僕は双方向リンクリストを作る事にしました。こんな感じの簡単な構造体を作ります：

struct LinkLine {
   Double2 s_, e_;    // 始点、終点
   LinkLine *next_;   // 次の辺
   LinkLine *pre_;    // 前の辺
   int32_t group_;    // グループ
}

自分の終点に繋がる辺を見つけたら、そのポインタをnext_に登録します。繋がる辺の方はpre_に自分を登録してもらいます。これで双方向リンクが完了です。これを黙々と続けていくと、いつか最初につなぎ始めた辺と連結されます。リング状になったという事です。リング状を確認出来たら適当なグループ番号をリングに所属するすべての辺に振って完了です。実際都道府県には陸続きになっていない島が沢山ありますので、一つの都道府県内のリンググループは複数になります。

　…と文章で書くと簡単なのですが、実際は結構メンドクサイコードを書く事になります。不具合と戦い試行錯誤を繰り返し、ついに出来た全47都道府県の白地図をご覧下さい！

　美しい～～。瀬戸内海や沖縄地方の小さな島々も全部閉じたN角形ポリゴンとして描画出来ているものです。全ての都道府県で頂点数は15万点くらいありますがw、世界地図日本の解像度で描くとその位の頂点数になるという事なんですねぇ。



⑪ 湖問題

　僕なりにこの段階で割と感無量だったりするのですが、実は上の白地図には大きな問題があります。それが「湖」。良く見ると、滋賀県に琵琶湖がありません。北海道に洞爺湖も支笏湖もありません。日本中の湖が無いんです。これはなぜか？

　上は滋賀県のMAPに今回の行政区画を重ねたものです。ご覧頂くと分かるように、実は行政区画の境界線は湖の上にもかかっています。これは考えてみると当然で、湖だって国内の土地の一つでして、自治体が管理する領域な訳です。つまり「地球地図日本の行政界」の情報からは湖の沿岸線を取る事が出来ないのです。なんてこったい orz。

　もう一つ、もし湖の沿岸線の情報が取れたとしても、N角形ポリゴンの中に穴を空けなければなりません。これはトポロジーやブーリアンの問題が絡んでくるため、ポリゴンとしてはとても面倒臭い部類のお話になってしまいます。残念ですが、今回はここまでが限界かなぁというのが僕の見解です。

　まぁ、でもここまでの知識と作ってきたコードがあれば、国土地理院が公開してくれている他のベクター形式のシェープファイルも同じプロセスで閲覧する事が出来ます。機会があればまたチャレンジしてみようと思います(^-^)