大学の授業（？）で Open Street Map (以下、OSMと記載) を触る機会がありました。
データの抽出方法の情報があまり見当たらなかったため、残しておきます。

今回は、GoとC++を使い、OpenStreetMapに含まれる道路情報を下の画像のようにJPG画像としてプロットしてみます。

東京都新宿駅周辺

コード及びDockerFile

ソースコードやOSM変換用のDockerファイルはGitHubにもアップしています。
合わせてご覧ください。

データを抽出するマップデータのURL取得

データがなければ何も始まらない為、初めにOSMのデータをダウンロードし、XML形式に変換します。 GeoFabrik より、.pbf形式のデータのURLを取得します。

Dockerを用いたXML形式への変換

今回は、Docker-Composeによるビルド時に、start.shに含まれるURLの地図データをダウンロードし、変換するようにしました。

続いて、バイナリ形式である.pbfデータを、XML形式である.osm形式に変換します。データ変換に関してはこのページに記載されています。

ただ、Windows版で変換しようとすると、文字化けしてしまい、XMLパースができませんでした。一方、Linux(64bit)版で変換すると、文字化けの問題もありませんでしたので、今回はLinux版を利用させていただきました。

私の環境はWindowsなので、Docker上に立ち上げたCentOS環境で変換を行いました。

start.sh

cd /osm
wget http://download.geofabrik.de/asia/japan/kanto-latest.osm.pbf
wget http://m.m.i24.cc/osmconvert64
./osmconvert64 kanto-latest.osm.pbf > kanto.osm

Dockerfile

FROM centos:centos8.1.1911
MAINTAINER Admin <admin@admin.com>

# "RUN" do at docker build
RUN dnf install -y wget
ADD ./start.sh /start.sh

CMD ./start.sh

docker-sompose.yml

version: '3'

services:
  centos8:
    build:
      context: .
      dockerfile: Dockerfile
    volumes:
        - ./osm:/osm

次のように、コンテナをビルドして実行すると、osmフォルダ内にXML形式である.osmファイルが生成されます。

docker-compose up -d

OSMのデータ構造

Node データ

<?xml version="1.0" encoding="UTF-8"?>
<osm version="0.6" generator="CGImap 0.7.5 (8774 thorn-02.openstreetmap.org)" copyright="OpenStreetMap and contributors" attribution="http://www.openstreetmap.org/copyright" license="http://opendatacommons.org/licenses/odbl/1-0/">
 <bounds minlat="35.6752400" minlon="139.7672300" maxlat="35.6790400" maxlon="139.7704900"/>
 <node id="31255230" visible="true" version="7" changeset="77771324" timestamp="2019-12-01T01:57:06Z" user="liuxinyu970226" uid="2008655" lat="35.6741521" lon="139.7714517">
  <tag k="highway" v="motorway_junction"/>
  <tag k="name" v="東銀座"/>
  <tag k="name:en" v="Higashi-ginza"/>
  <tag k="ref" v="5"/>
 </node>
 <node id="31255232" visible="true" version="4" changeset="74366018" timestamp="2019-09-11T17:30:45Z" user="penguin_the_dawn" uid="2929244" lat="35.6743695" lon="139.7709367"/>
 <node id="31255233" visible="true" version="3" changeset="74366018" timestamp="2019-09-11T17:30:45Z" user="penguin_the_dawn" uid="2929244" lat="35.6747901" lon="139.7698160"/>
 <node id="236621453" visible="true" version="16" changeset="34687216" timestamp="2015-10-17T04:55:48Z" user="muramototomoya" uid="610947" lat="35.6793882" lon="139.7716723">
  <tag k="highway" v="traffic_signals"/>
  <tag k="name" v="日本橋三丁目"/>
  <tag k="name:en" v="Nihonbashi 3"/>
  <tag k="traffic_signals" v="signal"/>
 </node>
 <node id="254366661" visible="true" version="6" changeset="15024546" timestamp="2013-02-14T01:01:25Z" user="Vlad" uid="24247" lat="35.6753419" lon="139.7716943">
  <tag k="highway" v="traffic_signals"/>
 </node>
 <!-- 後略 -->

Edge データ

 <!-- 前略 -->
 <node id="7040633094" visible="true" version="1" changeset="78091429" timestamp="2019-12-07T16:51:17Z" user="Diego Sanguinetti" uid="309382" lat="35.6747193" lon="139.7733165"/>
 <node id="7040633095" visible="true" version="1" changeset="78091429" timestamp="2019-12-07T16:51:17Z" user="Diego Sanguinetti" uid="309382" lat="35.6747667" lon="139.7733395"/>
 <way id="4849049" visible="true" version="11" changeset="70272861" timestamp="2019-05-15T12:15:41Z" user="Francesco_Citoni" uid="9396136">
  <nd ref="499185617"/>
  <nd ref="6114407385"/>
  <nd ref="499185467"/>
  <tag k="highway" v="unclassified"/>
  <tag k="maxspeed" v="30"/>
  <tag k="oneway" v="yes"/>
  <tag k="source" v="Bing 2010"/>
  <tag k="surface" v="asphalt"/>
 </way>
 <way id="23770583" visible="true" version="33" changeset="66736808" timestamp="2019-01-29T12:29:38Z" user="Diego Sanguinetti" uid="309382">
  <nd ref="1893305687"/>
  <nd ref="254366661"/>
  <nd ref="3790001996"/>
  <nd ref="6244220935"/>
  <nd ref="1105125752"/>
  <nd ref="6244220931"/>
  <nd ref="1984364215"/>
  <nd ref="6240643728"/>
  <nd ref="689881127"/>
  <nd ref="1130812185"/>
  <nd ref="6240643719"/>
  <nd ref="499185708"/>
  <nd ref="6227268576"/>
  <nd ref="2135828685"/>
  <nd ref="4114751664"/>
  <nd ref="499185715"/>
  <tag k="highway" v="tertiary"/>
  <tag k="lanes" v="3"/>
  <tag k="maxspeed" v="50"/>
  <tag k="name" v="鍛冶橋通り"/>
  <tag k="name:en" v="Kajibashi-dori"/>
  <tag k="name:es" v="Calle Kajibashi"/>
  <tag k="oneway" v="yes"/>
  <tag k="source" v="Bing 2010"/>
  <tag k="surface" v="asphalt"/>
 </way>
 <way id="24402716" visible="true" version="13" changeset="69094805" timestamp="2019-04-10T19:24:00Z" user="Diego Sanguinetti" uid="309382">
  <nd ref="264873291"/>
  <!-- 後略 -->

XMLをパースし、CSV形式に変換する

GoでXMLを処理するなどの記事を参考に、.osmファイルをパースしていきました。

Node データをCSVに書き出すコード

package main

import (
    "encoding/csv"
    "encoding/xml"
    "fmt"
    "os"
    "log"
    "io/ioutil"
)

type OSM struct {
    XMLName xml.Name  `xml:"osm"`
    Xml     string    `xml:",innerxml"`
    Nodes   []Node `xml:"node"`
}

type Node struct {
    Id  string `xml:"id,attr"`
    Lat string `xml:"lat,attr"`
    Lon string `xml:"lon,attr"`
}

func main() {
    xmlFile, err := os.Open("kanto.osm")
    if err != nil {
        log.Fatal(err)
        return
    }
    defer xmlFile.Close()
    xmlData, err := ioutil.ReadAll(xmlFile)
    if err != nil {
        log.Fatal(err)
        return
    }
    var data OSM
    if err := xml.Unmarshal(xmlData, &data); err != nil {
        fmt.Println("XML Unmarshal error:", err)
        return
    }

    file2,err2 := os.OpenFile("node.csv", os.O_WRONLY|os.O_CREATE, 0600)
    defer file2.Close()

    if err2 != nil {
        log.Fatal(err)
        return
    }
    writer := csv.NewWriter(file2)
    for itm := 0;itm<len(data.Nodes);itm++ {
        writer.Write([]string{data.Nodes[itm].Id, data.Nodes[itm].Lat,data.Nodes[itm].Lon})
    }
    writer.Flush()
}

このコードにより、次のようなCSVファイルを生成することができます。

NodeId , Node Lat , Node Lon

Go によるXMLのパースはかなりメモリを使うらしく、5.03GBほどのXMLデータのパースを始めると、ぐんぐんとメモリ使用量が増加しました。

ここまでくると、もう圧巻です。

Edge データをCSVに書き出すコード

package main

import (
"encoding/csv"
    "encoding/xml"
	"fmt"
	"os"
	"log"
	"io/ioutil"
)

type OSM struct {
XMLName     xml.Name  `xml:"osm"`
Xml         string    `xml:",innerxml"`
	Ways	[]Way `xml:"way"`
}

type Way struct {
  Id string `xml:"id,attr"`
  Nds []nd `xml:"nd"`
}

type nd struct {
    Ref string `xml:"ref,attr"`
}

func main() {
    xmlFile, err := os.Open("kanto.osm")
    if err != nil {
        log.Fatal(err)
        return
    }
    defer xmlFile.Close()
    xmlData, err := ioutil.ReadAll(xmlFile)
    if err != nil {
        log.Fatal(err)
        return
    }
//fmt.Println(xmlData);
    var data OSM
    if err := xml.Unmarshal(xmlData, &data); err != nil {
        fmt.Println("XML Unmarshal error:", err)
        return
    }

    file2,err2 := os.OpenFile("edge.csv", os.O_WRONLY|os.O_CREATE, 0600)
//    failOnError(err2)
    defer file2.Close()

//    err2 = file2.Truncate(0) // ファイルを空っぽにする(実行2回目以降用)
//    failOnError(err2)
    if err2 != nil {
        log.Fatal(err)
        return
    }
    writer := csv.NewWriter(file2)
    for itm := 0;itm<len(data.Ways);itm++ {
        for itmnd := 0;itmnd<len(data.Ways[itm].Nds)-1;itmnd++{
            writer.Write([]string{data.Ways[itm].Id, data.Ways[itm].Nds[itmnd].Ref, data.Ways[itm].Nds[itmnd+1].Ref})
        }
    }
    writer.Flush()

//fmt.Println(data.XMLName)
//    fmt.Println(data.Name)
//    fmt.Println(data.Nodes)
//    fmt.Println(data.Nodes[0].Lat)
 //   fmt.Println(data.osmdata.Nodes[1].lon)
}

このコードにより、次のようなCSVファイルを生成することができます。

WayId , Node1 , Node2

Tagデータのタグ情報の抽出

さらに、道路情報のみを抜き出すために、wayに含まれるタグ情報を取得していきます。

package main

import (
    "encoding/csv"
    "encoding/xml"
    "fmt"
    "os"
    "log"
    "io/ioutil"
)

type OSM struct {
    XMLName xml.Name  `xml:"osm"`
    Xml     string    `xml:",innerxml"`
    Ways    []Way     `xml:"way"`
}

type Way struct {
  Id    string  `xml:"id,attr"`
  Tags   []tag    `xml:"tag"`
}

type tag struct {
    K string `xml:"k,attr"`
    V string `xml:"v,attr"`
}

func main() {
    xmlFile, err := os.Open("kanto.osm")
    if err != nil {
        log.Fatal(err)
        return
    }
    defer xmlFile.Close()
    xmlData, err := ioutil.ReadAll(xmlFile)
    if err != nil {
        log.Fatal(err)
        return
    }
    var data OSM
    if err := xml.Unmarshal(xmlData, &data); err != nil {
        fmt.Println("XML Unmarshal error:", err)
        return
    }

    file2,err2 := os.OpenFile("tag.csv", os.O_WRONLY|os.O_CREATE, 0600)
    defer file2.Close()

    if err2 != nil {
        log.Fatal(err)
        return
    }
    writer := csv.NewWriter(file2)
    for itm := 0;itm<len(data.Ways);itm++ {
        for itmnd := 0;itmnd<len(data.Ways[itm].Tags);itmnd++{
            writer.Write([]string{data.Ways[itm].Id, data.Ways[itm].Tags[itmnd].K, data.Ways[itm].Tags[itmnd].V})
        }
    }
    writer.Flush()
}

NodeデータとEdgeデータの統合

続いて、NodeデータとEdgeデータの統合を行います。
ここまでGoで変換してきたので、本来はGoを利用するべきかとも思うのですが、慣れているC++の方に逃げました…。

# include <string>
# include <vector>
# include <iostream>
# include <cstdio>
# include <fstream>
# include <cstring>

using namespace std;

class mycsv {
private:
    FILE* fp;
    int bufsize;
    char* line;
    char comma;
    void init() {
        fp = NULL;
        line = new char[bufsize];
        comma = ',';
    }
public:
    // set comma charcter
    void setcomma(char c) {
        comma = c;
    }
    // open csv file
    bool open(string filename) {
        if (fp != NULL) fclose(fp);
        fp = fopen(filename.c_str(), "r");
        if (fp == NULL) return false;
        return true;
    }
    // close csv file
    void close() {
        if (fp != NULL) fclose(fp);
        fp = NULL;
    }
    vector<string> loadnextline() {
        memset(line, bufsize, sizeof(char) * bufsize);
        fgets(line, bufsize, fp);

        size_t len = strlen(line);
        line[len - 1] = comma;
        line[len] = '\0';

        vector<string> cells;
        string l;

        for (int i = 0; i < len; i++) {
            if (line[i] == '\r') continue;
            if (line[i] == comma) {
                cells.push_back(l);
                l.clear();
                continue;
            }
            l.push_back(line[i]);
        }
        return cells;
    }
    bool eof() {
        if (fp == NULL) return true;
        return feof(fp);
    }
    vector<vector<string>> loadall() {
        vector<vector<string>> ans;
        while (!eof())
        {
            ans.push_back(loadnextline());
        }
        return ans;
    }
    mycsv(string filename, int linemax = 4096) {
        bufsize = linemax;
        init();
        open(filename);
    }
    mycsv(int linemax) {
        bufsize = linemax;
        init();
    }
    ~mycsv() {
        delete[] line;
        if (fp == NULL) return;
        fclose(fp);
        fp = NULL;
    }
};
# include <unordered_map>
# include <unordered_set>

typedef pair<string, string> node;

int main() {

    unordered_set <string> g_street_wayids;

    // ノード情報の読み取り
    mycsv csv("node.csv");
    unordered_map<string, node> nodes;
    while (!csv.eof())
    {
        vector<string> cells = csv.loadnextline();
        if (cells.size() < 3)continue;
        nodes[cells[0]] = node(cells[1], cells[2]);
    }

    // タグ情報の読み取り
    csv.open("tag.csv");
    while (!csv.eof())
    {
        vector<string> cells = csv.loadnextline();
        if (cells.size() < 3)continue;
        if (cells[1] == "highway") g_street_wayids.insert(cells[0]);
        else if (cells[1] == "highway") g_street_wayids.insert(cells[0]);
    }

    // way情報の読み取り
    ofstream ofs("combind.csv");
    csv.open("edge.csv");
    while (!csv.eof())
    {
        vector<string> cells = csv.loadnextline();
        if (cells.size() < 3)continue;
        if (nodes.find(cells[1]) == nodes.end()) continue;
        if (nodes.find(cells[2]) == nodes.end()) continue;
        if (g_street_wayids.find(cells[0]) == g_street_wayids.end()) continue;
        ofs << cells[0] << ',' << cells[1] << ',' << nodes[cells[1]].first << ',' << nodes[cells[1]].second << ',' << cells[2] << ',' << nodes[cells[2]].first << ',' << nodes[cells[2]].second << endl;
    }
}

Node を画像にプロットする

道路データが取得できたので、nodeデータを座標平面上にプロットしてみました。画像をC++で扱うのは、環境構築が面倒なので、再びGoに回帰します。

package main

import (
    "os"
//  "math"
    "image"
    "image/jpeg"
    "image/color"
    "strconv"
    "encoding/csv"
//  "fmt"
)

const scale = 100000
const pie = 3.14159265359

const area_minlat,area_minlon,area_maxlat,area_maxlon=35.55,139.6,35.8,139.9

var minlat,minlon,maxlat,maxlon int
var datamap,datamap2 [][]int
var coun,num_thread int

func main() {

    minlat,minlon,maxlat,maxlon = area_minlat*scale,area_minlon*scale,area_maxlat*scale,area_maxlon*scale

    // メモリ確保
    datamap = make([][]int,int(maxlat-minlat))
    for i := range datamap{
        datamap[i]=make([]int,int(maxlon-minlon))
    }
    datamap2 = make([][]int,int(maxlat-minlat))
    for i := range datamap2{
        datamap2[i]=make([]int,int(maxlon-minlon))
    }

    loaddatafile()

    // 画像処理用
    x := 0
    y := 0
    width := maxlon-minlon
    height := maxlat-minlat

    // RectからRGBAを作る(ゼロ値なので黒なはず)
    img := image.NewRGBA(image.Rect(x, y, width, height))

    imagefile, _ := os.Create("map.jpg")
    defer imagefile.Close()

    for i:=0;i<height;i++{
        for j:=0;j<width;j++{
            lonx:=j
            laty:=height-i
            col := datamap2[i][j]
//          col := math.Log10(float64(datamap2[i][j]+1))*100

            if col > 0{
                col = 255
            }
            img.Set(lonx, laty, color.RGBA{uint8(col), uint8(col), uint8(col), 1})
        }
    }

    // JPEGで出力(100%品質)
    if err := jpeg.Encode(imagefile, img, &jpeg.Options{100}); err != nil {
        panic(err)
    }
}

func numthreadfunc(){
    num_thread--
}

func loaddatafile(){
    file, err := os.Open("combind.csv")
    if err != nil {
        panic(err)
    }
    defer file.Close()

    reader := csv.NewReader(file)
    var line []string

    for {
        line, err = reader.Read()
        if err != nil {
            break
        }
        var aflat,aflon,bflat,bflon float64
        aflat, _ = strconv.ParseFloat(line[2],64)
        aflon, _ = strconv.ParseFloat(line[3],64)
        bflat, _ = strconv.ParseFloat(line[5],64)
        bflon, _ = strconv.ParseFloat(line[6],64)

        if aflat < area_minlat || area_maxlat < aflat || aflon < area_minlon || area_maxlon < aflon{
            continue
        }
        if bflat < area_minlat || area_maxlat < bflat || bflon < area_minlon || area_maxlon < bflon{
            continue
        }

        aintlat := int(aflat*scale)
        aintlon := int(aflon*scale)
        bintlat := int(bflat*scale)
        bintlon := int(bflon*scale)

        if aintlat >= maxlat || aintlon >= maxlon{
            continue
        }
        if aintlat < minlat || aintlon < minlon{
            continue
        }
        if bintlat >= maxlat || bintlon >= maxlon{
            continue
        }
        if bintlat < minlat || bintlon < minlon{
            continue
        }

        dlat,dlon := aintlat-bintlat,aintlon-bintlon;
        if dlat < 0 {
            dlat = -dlat
        }
        if dlon < 0 {
            dlon = -dlon
        }
        ddd:=dlat
        if dlon > dlat{
            ddd = dlon
        }
        if ddd==0{
            continue
        }
        for ilat:=0;ilat<=ddd;ilat++{
            ilon:=ilat
            tlat := int(ilat*aintlat/ddd+(ddd-ilat)*bintlat/ddd)
            tlon := int(ilon*aintlon/ddd+(ddd-ilon)*bintlon/ddd)

/*          fmt.Println("aa")
            fmt.Println(aintlat)
            fmt.Println(bintlat)
            fmt.Println("bb")*/

            if tlat-minlat < 0 || tlon-minlon < 0 {
                continue
            }
            if maxlat-minlat <= tlat-minlat || maxlat-minlat <= tlon-minlon{
                continue
            }

            datamap[tlat-minlat][tlon-minlon]++
            datamap2[tlat-minlat][tlon-minlon]++
        }

//      datamap[aintlat-minlat][aintlon-minlon]++
//      datamap2[aintlat-minlat][aintlon-minlon]++
    }
    num_thread--
}

プロットした結果、下の画像のように、道路を線としてあらわすことができました。

拡大図