Go to Qiita Advent Calendar 2024 Top
LoginSignup
12
11

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@mima_ita

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

Posted at

国土数値情報の鉄道データで取得できる鉄道路線を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://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で表現できます。でも地下鉄とか新幹線は、対応できないこのプログラムはできそこないだ、使えないよ！

12
11
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
12
11

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?