へっぽこ実験ブログ

国土数値情報の鉄道データを3Dで表示してみる

投稿日: 投稿者: mima3
Table of Content

国土数値情報の鉄道データで取得できる鉄道路線を3Dで表示してみます。

デモ
http://needtec.sakura.ne.jp/threemap/railroad3d_example.html

クライアント側 ソース
https://github.com/mima3/threemap

サーバ側 ソース
https://github.com/mima3/kokudo

路線2.png

サーバーサイドの処理

サーバーサイドの処理は格納済みの国土数値情報の鉄道データを要求された条件で返します。
この際、標高データを付与しています。

国土数値情報から鉄道データをSpatialiteにインポートする

国土数値情報の鉄道情報をSpatialiteにインポートします。

spatialiteの操作は以下のコードになります。
https://github.com/mima3/kokudo/blob/master/kokudo_db.py

これを国土数値情報の鉄道データのshapeファイルをデータベースにインポートするには次のコマンドを実行してください。

python import_railroad_section.py C:\tool\spatialite\mod_spatialite-4.2.0-win-x86\mod_spatialite.dll test.sqlite original_data\N02-13\N02-13_RailroadSection.shp

HTTP経由でGeoJsonを取得できるようにする。

Bottleを使用してHTTP経由で作成したデータベースから指定の路線のGeojsonを取得できるようにします。

呼び出し例
http://needtec.sakura.ne.jp/kokudo/json/get_railroad_section?operationCompany=%E6%9D%B1%E4%BA%AC%E6%80%A5%E8%A1%8C%E9%9B%BB%E9%89%84

この例では東京急行電鉄の路線を全て取得しています。

地理院地図の標高APIを使用して標高を求める

国土数値情報の鉄道データには経度、緯度はありますが、標高はありません。
これでは３Dで表示する意味がありません。

そこで、経度緯度から標高を取得できるようにします。
これには、地理院地図の標高APIを使用します。

http://portal.cyberjapan.jp/help/development/api.html

使用例：
http://cyberjapandata2.gsi.go.jp/general/dem/scripts/getelevation.php?lon=140.08531&lat=36.103543&outtype=JSON

結果 

{"elevation":25.3,"hsrc":"5m\uff08\u30ec\u30fc\u30b6\uff09"}

しかし、ながら、毎回このAPIを実行するのも無駄なのでデータベースに一度読み込んだ内容をキャッシュするようにします。

https://github.com/mima3/kokudo/blob/master/gsi_api.py 

# -*- coding: utf-8 -*-
import urllib
import urllib2
from peewee import *
from playhouse.sqlite_ext import SqliteExtDatabase
import json

database_proxy = Proxy()  # Create a proxy for our db.

def get_elevation_by_api(long, lat):
    """
    標高値の取得
    http://portal.cyberjapan.jp/help/development/api.html
    """
    url = ('http://cyberjapandata2.gsi.go.jp/general/dem/scripts/getelevation.php?lon=%f&lat=%f&outtype=JSON' % (long, lat))
    req = urllib2.Request(url)
    opener = urllib2.build_opener()
    conn = opener.open(req)
    cont = conn.read()
    ret = json.loads(cont)
    return ret

def str_isfloat(str):
    try:
        float(str)
        return True
    except ValueError:
        return False

class ElevationCache(Model):
    """
    標高のキャッシュテーブル
    """
    lat = FloatField(index=True)
    long = FloatField(index=True)
    hsrc = TextField(index=True)
    elevation = FloatField(null=True)

    class Meta:
        database = database_proxy

def connect(path):
    db = SqliteExtDatabase(path)
    database_proxy.initialize(db)

def setup(path):
    connect(path)
    database_proxy.create_tables([ElevationCache], True)

def get_elevation(long, lat):
    try:
        ret = ElevationCache.get((ElevationCache.long==long) & (ElevationCache.lat==lat))
        return {'elevation': ret.elevation, 'hsrc': ret.hsrc}

    except ElevationCache.DoesNotExist:
        ret = get_elevation_by_api(long, lat)
        elevation = ret['elevation']
        if not str_isfloat(elevation):
            elevation = None
        ElevationCache.create(
            long = long,
            lat = lat,
            elevation = elevation,
            hsrc = ret['hsrc']
        )
        return ret

class GsiConvertError(Exception):
    def __init__(self, value):
        self.value = value
    def __str__(self):
        return repr(self.value)

def convert_geojson(json):
    for feature in json['features']:
        if feature['geometry']['type'] == 'LineString':
            prop = {}
            start = feature['geometry']['coordinates'][0]
            end = feature['geometry']['coordinates'][len(feature['geometry']['coordinates'])-1]
            start_elevation = get_elevation(start[0], start[1])
            end_elevation = get_elevation(end[0], end[1])
            feature['properties']['start_elevation'] = start_elevation['elevation']
            feature['properties']['end_elevation'] = end_elevation['elevation']
        else:
            raise GsiConvertError('unexpected feature type')
    return json

if __name__ == '__main__':
    setup('elevation_cache.sqlite')
    #print get_elevation_by_api(133, 39)
    #print get_elevation_by_api(139.766084, 35.681382)

    print get_elevation(133, 39)
    print get_elevation(139.766084, 35.681382)
    with open('get_railroad_section.geojson' , 'rb') as f:
        cont = f.read()
        print convert_geojson(json.loads(cont))

get_elevation関数は経度、緯度から標高を求めます。この際、DBに値が格納済みならそれを返し、格納されていなければ標高APIを実行し、その結果をDBに格納して返します。

convert_geojson関数はLineStringのGeoJsonに標高を付与します。
線の開始点と終了点に対して標高を取得して、その結果をstart_elevation,end_elevationとしてプロパティに格納します。

標高データ付きの鉄道データのGeoJSONを取得する。

使用例
http://needtec.sakura.ne.jp/kokudo/json/get_railroad_section?operationCompany=%E6%9D%B1%E4%BA%AC%E6%80%A5%E8%A1%8C%E9%9B%BB%E9%89%84&embed_elevation=True

embed_elevationを付与することで標高が取得できます。

取得結果 

{
  "type": "FeatureCollection", 
  "features": [
    {
      "geometry": {
        "type": "LineString", 
        "coordinates": [[139.48677, 35.55760999999999], [139.4865599999999, 35.55839]]
      }, 
      "type": "Feature", 
      "properties": {
        "operationCompany": "\u6771\u4eac\u6025\u884c\u96fb\u9244",
        "end_elevation": 36.7, 
        "serviceProviderType": "4", 
        "railwayLineName": "\u3053\u3069\u3082\u306e\u56fd\u7dda",
        "railwayType": "12", 
        "start_elevation": 35.4
      }
    },// 略
  ]
}

クライアントサイドの処理

クライアント側の処理は大きく次の通りです。
１．取得したGeoJSONのz軸を付与する。この際、開始点と終了点はpropertiesで指定した標高、それ以外は、開始点と終了点を考慮した割合とする。
２．結合できる線は結合する
３．three.jsのTubeGeometryを用いて路線の描画

取得したGeoJSONのz軸を付与する

以下のようにpropertiesのstart_elevationとend_elevationを使用して全ての座標に対してz軸に標高を追加します。

    function expendElevation(features) {
      for (var i = 0; i < features.length; ++i) {
        var feature = features[i];
        var end_elevation = feature.properties.end_elevation;
        var start_elevation = feature.properties.start_elevation;
        var per_elevation = (end_elevation - start_elevation) / feature.geometry.coordinates.length;
        for (var j = 0; j < feature.geometry.coordinates.length; ++j) {
            feature.geometry.coordinates[j].push(start_elevation + (j * per_elevation));
        }
      }
      return features;
    }

結合できる線は結合する

たとえば、あるfeatureの開始点と別のfeatureの終了点が等しかった場合、これを結合して一つの線とみなします。

    function compressionLine(features) {
      var before = features.length;
      for (var i = 0; i < features.length; ++i) {
        for (var j = features.length -1; i < j; --j) {
          var f1Start = features[i].geometry.coordinates[0];
          var f1End = features[i].geometry.coordinates[features[i].geometry.coordinates.length-1];

          var f2Start = features[j].geometry.coordinates[0];
          var f2End = features[j].geometry.coordinates[features[j].geometry.coordinates.length-1];

          // f1の開始点がf2の終端と一致する場合、f1の前にf2がある
          if (f1Start[0] == f2End[0] && f1Start[1] == f2End[1]) {
            features[i].geometry.coordinates = features[j].geometry.coordinates.concat(features[i].geometry.coordinates);
            features.splice(j, 1);
            break;
          }
          // f1の終了点がf2の開始点と一致する場合、f1の後にf2がある
          if (f1End[0] == f2Start[0] && f1End[1] == f2Start[1]) {
            features[i].geometry.coordinates = features[i].geometry.coordinates.concat(features[j].geometry.coordinates);
            features.splice(j, 1);
            break;
          }
        }
      }
      if (features.length == before) {
        return features;
      }
      return compressionLine(features);
    }

TubeGeometryを用いて路線の描画

TubeGeometryを用いて路線を描画します。
開始点の座標をオフセットとして、mesh.positionの位置を決めます。
他の点については、TubeGeometryのパスとして追加する際に、開始点との相対座標を付与するようにします。

    function createRailroad(geodata) {
      console.log(geodata.length);
      geodata = expendElevation(geodata);
      geodata = compressionLine(geodata);
      console.log(geodata.length);
      var scaleElevation = 100;
      for (var i = 0 ; i < geodata.length ; i++) {
        var lineList = [];
        var geoFeature = geodata[i];
        var baseX;
        var baseY;
        for (var j = 0; j < geoFeature.geometry.coordinates.length; ++j) {
          var pt = reverseProjection([
            geoFeature.geometry.coordinates[j][0],
            geoFeature.geometry.coordinates[j][1]
          ]);
          if (j ==0) {
            baseX = pt[0];
            baseY = pt[1];
            lineList.push(new THREE.Vector3(0, 0, geoFeature.geometry.coordinates[j][2] / scaleElevation));
          } else {
            lineList.push(new THREE.Vector3(pt[0] - baseX, pt[1] - baseY, geoFeature.geometry.coordinates[j][2] /scaleElevation));
          }
        }
        var spline = new THREE.SplineCurve3(lineList);
        var tubeGeo = new THREE.TubeGeometry(spline, 32, 0.03, 8, false);
        var mesh = new THREE.Mesh(
          tubeGeo,
          new THREE.MeshLambertMaterial( { 
            color: 0xff0000,
            transparent: true,
            opacity: 0.9
          })
        );
        mesh.position.set(baseX, baseY, 0.1);
        scene.add(mesh);
      }
    }

この時、経度緯度から座標を取得するのに使用しているreverseProjectionは次のようにd3.geo.path()のprojectionを利用しています。

reverseProjectionの内容 

    //上下が反転しているのでムリくり合わせる
    var reverseProjection = function(x, y) {
      pt = projection(x, y);
      pt[1] *= -1;
      return pt;
    };

    //geoJSONのデータをパスに変換する関数を作成
    var path = d3.geo.path().projection(reverseProjection);

まとめ

国土数値情報の鉄道データをspatialiteに格納することで、動的にGeojsonが作成できるようになります。

地理院の標高APIを使えば、緯度、経度から標高が取得できます。この際、結果はキャッシュしておいたほうがいいです。

three.jsを使えば簡単に3D表現がおこなえ、この際に、d3のprojectionを使うと経度緯度の変換が楽に行えます。

これらを利用すれば、鉄道の路線を3Dで表現できます。でも地下鉄とか新幹線は、対応できないこのプログラムはできそこないだ、使えないよ！

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です