1. はじめに

都市の姿は様々です。地形の制約を受けつつも何らかの思想や意図のもとで計画がなされ、技術の進歩に応じて再開発が繰り返されます。建物の入れ替わりや区画整理こそあれ、一度開拓された都市は大域的にはその構造を維持し続けることが多いでしょう。都市の様態には地理的・政治的・歴史的背景が詰まっています。

飛行機から見下ろすと、そのような都市の姿が明確に見て取れます。整然と格子状に区切られた京都の町・凱旋門を中心として放射状に広がるパリの通り・迷路のように入り組んだマラケシュのメディナ・山と海に挟まれた傾斜地……

ここで注目したいのが、こういった都市の特徴はいずれも道の方向という形で現れるという点です。ある地域で道路の方向がどのような性質を持っているかを可視化することで、このような都市の姿を明らかにしてみたいと思います。

なお、本記事の可視化手法はMapboxのブログに紹介されているものです。見せ方のアイデアはその記事に依りますが、ソースコードは特に参考にせず自前で実装しています。

2. 例

細かい実装の話は後回しにして、様々な地域で道路の方向を可視化した結果を見てみましょう ٩(๑❛ᴗ❛๑)۶

以下では地図上の青色の線が道路で、右上に道路の方向ヒストグラムが表示されています。

道路の「方向」は進行方向を考慮したものではなく180°周期で考えており、ヒストグラムも必ず180°対称になっています（「南 → 北」と「北 → 南」を区別しない）。

2.1. 京都駅周辺（碁盤の目）

京都駅周辺は道路が碁盤の目のように張り巡らされており、ヒストグラムも上下左右に突出しています。THE・碁盤の目という感じですね。

2.2. 長野県須坂市（迷路）

「迷路のような街」でググったところ、以下の記事で長野県須坂市が紹介されていました。

迷路の街で道に迷うと夢みたいな展開になる :: デイリーポータルZ

京都とは対照的に様々な方向の道路があることがわかります。

2.3. 彦根（城下町）

城下町は当然ながら城の影響を大きく受けており、彦根城と同方向の道路が多くなっています。

ズームアウトして見てみると、そもそも彦根城の方向自体が琵琶湖の存在に影響されているように推察されます。

2.4. 多摩川（川辺）

川の近くはやはり川の方向に大きな制約を受けていますね。かつてはこの辺りも田園地帯だったのかもしれません。

2.5. 神戸（海と山の間）

山と海に挟まれた神戸の街は南北に大きく傾斜しています。市民の多くが「山側」「海側」という形で方向を認識しており、線路は基本的に東西方向に存在します。

3. 実装

3.1. 基本方針

ブラウザ上に地図が表示されているとして、画面内地域の道路方向を可視化することが目的です。

明らかに必要なデータは地図上での道路に関するデータです。国土地理院が公開している道路中心線のベクターデータを利用します。

github.com

このデータを利用して、以下のような手順を踏めばよいでしょう。

画面上に表示されている地域の範囲を取得
範囲内における道路のベクターデータを取得
道路のベクターデータをまずは地図上にそのまま表示
道路（折れ線として座標が与えられている）の線分ごとにその方向を計算し、方向のヒストグラムを描画する

ヒストグラムは方向に関するものなので、円状にして方向に依存した色相を与えることでより直感的になります。方向を色相に対応づけるという考え方は方向と色相の周期性に着目したもので、複素関数の可視化でも用いた典型的なアイデアです。

sw1227.hatenablog.com

3.2. 今回の実装の特徴（フレームワーク等）

私の実装は、Mapboxなどで紹介されているものと比較すると以下のような特徴が存在します。

道路情報は国土地理院から取得しているため、日本国内で可視化ができる
地図ライブラリはMapboxではなくLeafletを使用
グラフの描画にはd3.jsを使用
JavaScriptのフレームワークとしてReactを使っている

3.3. 詳細

実装上重要な点をピックアップして詳しく解説します。

注意: 以下に記載するJavaScriptのコードはReactコンポーネントの一部だったりするため、そのままコピペしても動かない可能性があります。

地図・道路

地図の表示にはLeafletというJavaScriptのライブラリを使用しています。表示する地図の種類は任意のタイルを指定できるため、今回は国土地理院の地理院タイル（標準地図）を使用しました。

道路のベクターデータは、先述のように国土地理院のベクトルタイル提供実験で提供されているGeoJSONタイルを使用します。タイル座標でリクエストをすればGeoJSONが得られます。GeoJSONはLeafletのデフォルト機能で簡単に地図上に重ねることが可能です。

画面上の領域を取得

Leafletではmap.getZoom(), map.getBounds()で地図上に表示されている領域の情報が得られるため、それをタイル座標に変換し、必要な道路中心線GeoJSONタイルを全て取得することになります。

以下に、

Leafletで地図上に表示されている地域を取得するメソッド
あるズームレベルのもとで、特定の緯度経度を含むタイルのタイル座標を計算する関数

のコードを記載します。これらを組み合わせて画面内の道路データが取得できます。

    getCurrentBounds = () => {
        // Get bounding box of tile coordinates in the current view
        const zoom = this.state.map.getZoom();
        const bounds = this.state.map.getBounds();
        return {
            "northEast": getTileCoords(bounds._northEast.lat, bounds._northEast.lng, zoom),
            "southWest": getTileCoords(bounds._southWest.lat, bounds._southWest.lng, zoom)
        };
    }

// Get tile coordinates from latitude, longitude and zoom
function getTileCoords(lat, lon, zoom) {
    const xTile = parseInt(Math.floor((lon + 180) / 360 * (1 << zoom)));
    const yTile = parseInt(Math.floor((1 - Math.log(Math.tan(lat * Math.PI / 180) + 1 / Math.cos(lat * Math.PI / 180)) / Math.PI) / 2 * (1 << zoom)));
    return { "z": zoom, "x": xTile, "y": yTile }
}

道路方向の計算

道路のGeoJSONは曲線ではなく折れ線なので、各線分の両端の座標から $arctan$ を計算すればOKです。

単純に緯度をy方向・経度をx方向として扱いました。

// Direction of vector p2-p1
function direction(p1, p2) {
    return Math.atan2(p2[1] - p1[1], p2[0] - p1[0]);
}

本来は線分の長さに応じて重み付けをした上でヒストグラムを作るべきなのかもしれませんが、今回は妥協して（高速化も考慮して）全ての線分を対等に扱うことにしました。

ヒストグラムの描画

d3.jsで頑張って描画します。色は方向に対して虹色のカラースケールを適用することで設定できます。

関連するメソッドは以下のようになります。

    drawHist = () => {
        const margin = { top: 20, right: 20, bottom: 20, left: 20 };
        const width = svgSize - margin.left - margin.right;
        const height = svgSize - margin.top - margin.bottom;

        const svg = d3.select(this.refs.svg)
          .append("g")
            .attr("transform", `translate(${margin.left + width / 2}, ${margin.top + height / 2})`);

        svg.append("circle")
            .attr("cx", 0).attr("cy", 0)
            .attr("r", this.maxRadius)
            .attr("style", "stroke: #aaa; stroke-width: 1; fill: none;");

        this.updateHist([]);
    }

    updateHist = data => {
        const histogram = d3.histogram()
            .value(d => d)
            .domain([-Math.PI, Math.PI])
            .thresholds(d3.range(-Math.PI, Math.PI, Math.PI / this.histSplit));
        const bins = histogram(data)

        const rScale = d3.scaleLinear()
            .domain([0, d3.max(bins.map(x => x.length))])
            .range([0, this.maxRadius]);

        const arc = d3.arc()
            .innerRadius(0).outerRadius(d => rScale(d.length))
            .startAngle(d => Math.PI / 2 - d.x0)
            .endAngle(d => Math.PI / 2 - d.x1);

        const path = d3.select(this.refs.svg).select("g")
            .selectAll("path")
            .data(bins);
        path.enter().append("path")
            .merge(path)
            .attr("d", arc)
            .attr("fill", (d, i) => d3.interpolateRainbow(i / this.histSplit));
    }

インタラクション

地図上を自由に移動し、気になる地点でボタンを押すと可視化が実行されるようにしました。タイルごとに非同期で道路データを取得・可視化するため、ボタンを押すとヒストグラムがにょきにょきと更新されていくはずです（下図）。

f:id:sw1227:20181226151251g:plain:w400 — 非同期で更新（GIF）

今回は実装しませんでしたが、地図を動かすたびにイベント駆動でヒストグラムを更新するという方法も考えられるでしょう。

4. 今後やりたいこと

タイル座標ごとに道路方向を要約した量（最頻方向・方向のばらつきなど）を計算し、その値に応じて色分けしたグリッドを地図に重ねて表示する。地域ごとの違いがより大域的に分かるはず。
ヒストグラムを描画する際、折れ線の線分数ではなく線分長で重み付けされた量を使用する。定性的にはあまり変わらないような…？
ベクトルタイルの仕様の関係でズームレベルが16のときのみヒストグラム描画を実行できるようにしているが、別のズームレベルにも対応する（ズームレベル間の座標変換を実装すれば良い）。