政府統計窓口(eStat)の地域メッシュを3Dで表示してみる

以前、政府統計窓口の地域メッシュをWebブラウザで表示しました。

政府統計窓口(eStat)の地域メッシュをWebブラウザで表示する方法
http://qiita.com/mima_ita/items/38784095a146c87dcd23

せっかくなので、three.jsを用いて地域メッシュを3Dで表示してみます。

demo
http://needtec.sakura.ne.jp/threemap/mesh3d_example.html

ソースコード
https://github.com/mima3/threemap

3dmesh.png

指定範囲のデータをJSONで取得できるようにする。

政府統計窓口の地域メッシュと、国土数値情報の行政区域を範囲を指定して取得できるようにします。

「細かいことはいいんだよ!俺は3Dで表示したいんだ!」って人は、以下の呼び出し例をコピペすればJSONが取得できるので、それを利用してください。

指定の範囲の地域メッシュをJSONで取得できるようにする。

まず、政府統計窓口の地域メッシュをHTTP経由でJSONで取得できるようにします。

呼び出し例:
http://needtec.sakura.ne.jp/estat/json/get_mesh_stat_group_by_mesh?swlat=35.45&swlng=139.4&nelat=35.7&nelng=139.81&stat_id=T000608

取得内容:

{
  "type": "FeatureCollection", 
  "features": [
    {
      "geometry": {
        "type": "Polygon", 
        "coordinates": [
          [
            [139.4124999999999, 35.45833333333333], 
            [139.4249999999999, 35.45833333333333], 
            [139.4249999999999, 35.46666666666666], 
            [139.4124999999999, 35.46666666666666], 
            [139.4124999999999, 35.45833333333333]
          ]
        ]
      }, 
      "type": "Feature", 
      "properties": {
        "世帯総数": "5018", 
        "女": "6093", "\u4eba\u53e3\u7dcf\u6570": "12349", 
        "男": "6256", "area": "53391353"
      }
    },// 略 
  ]
} 

ソースコード
https://github.com/mima3/estat

この実装はPython+bottleで行っています。

指定の範囲の行政区域をGeoJSONで取得できるようにする。

同様に国土数値情報の行政区域も範囲を指定してGeoJSONで取得できるようにします。

呼び出し例:
http://needtec.sakura.ne.jp/kokudo/json/get_administrative_district_by_geometry?swlat=35.45&swlng=139.4&nelat=35.7&nelng=139.81

取得内容:

{
  "type": "FeatureCollection", 
  "features": [
    {
      "geometry": {
        "type": "Polygon", 
        "coordinates": [
          [
            [139.4478029999999, 35.42797099999999], 
            [139.445154, 35.428581], // 略
          ]
        ]
      }, 
      "type": "Feature", 
      "properties": {
        "prefectureName": "\u795e\u5948\u5ddd\u770c", 
        "administrativeAreaCode": "14218", 
        "cityName": "\u7dbe\u702c\u5e02", 
        "countryName": "                    ", 
        "subPrefectureName": "                    ", 
        "PK_UID": 94566
      }
    }, // 略
  ]
}

ソースコード
https://github.com/mima3/kokudo

この実装はPython+bottleで行っています。

three.jsを使用して3Dで表示する

D3.jsで読み込んだ行政区域のGeojsonを3D表示する。

GeojsonをD3.jsで読み込み、作成されたSVGをthree.jsで3D表示します。
この方法は、下記のページを参考にしました。

GUNMA GIS GEEK D3.jsで読み込んだgeoJSONデータをthree.jsを使って3D表示する。
http://shimz.me/blog/d3-js/3137

Render geographic information in 3D with Three.js and D3.js
http://www.smartjava.org/content/render-geographic-information-3d-threejs-and-d3js

まず、JSONをD3.jsを利用して取得します。
この際、async.jsを用いて並列でJSONをリクエストしています。

async.parallel([
  function (callback) {
    d3.json("http://needtec.sakura.ne.jp/kokudo/json/get_administrative_district_by_geometry?swlat=35.45&swlng=139.4&nelat=35.7&nelng=139.81", function(json) {
        callback(null, json);
    });
  },
  function (callback) {
    d3.json("http://needtec.sakura.ne.jp/estat/json/get_mesh_stat_group_by_mesh?swlat=35.45&swlng=139.4&nelat=35.7&nelng=139.81&stat_id=T000608", function(json) {
        callback(null, json);
    });
  }
], function(err, results) {
  // 略
});

この取得したgeojsonをsvg pathに変換して、その後、three.jsに表示します。

  function initScene() {
    container = $('#map');

    //シーンの追加
    scene = new THREE.Scene();

    //カメラの設定
    camera = new THREE.PerspectiveCamera(70, window.innerWidth / window.innerHeight, 1, 10000);
    camera.position.set(-6.005123939520757, -6.754277108005312, 3.1756149676051133);
    var target = new THREE.Vector3();
    target.x = -7.162140160089945;
    target.y = 0.42803104270646924;
    target.z = -2.063747120298727;
    //camera.lookAt(target); // TrackballControls使用中は効かない

    scene.add(camera);

    //ライティングの設定
    var light = new THREE.DirectionalLight(0xffffff, 2);
    light.position.set(-10, -20, 30).normalize();
    scene.add(light);

    light2 = new THREE.AmbientLight(0x333333);               
    scene.add(light2);  

    var projection = d3.geo.equirectangular()
       .scale(900)
       .translate([-2200, 560]);

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

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

    createAdministrativeDistrictMesh(results[0].features);
    function createAdministrativeDistrictMesh(geodata) {
      //geoJSON→svg path→three.js mesh変換    
      var countries = [];
      for (i = 0 ; i < geodata.length ; i++) {
        var geoFeature = geodata[i];
        var properties = geoFeature.properties;
        var feature = path(geoFeature);

        //svgパスをthree.jsのmeshに変換
        var gmesh = transformSVGPathExposed(feature);
        //console.log(properties, gmesh);
        countries.push({"data": properties, "mesh": gmesh});
      }

      //mesh追加
      for (i = 0 ; i < countries.length ; i++) {

          var material = new THREE.MeshLambertMaterial({color:0x888888});

          var shape3d = countries[i].mesh.extrude({
              amount: 0,
              bevelEnabled: false
          });

          var toAdd = new THREE.Mesh(shape3d, material);
          //toAdd.rotation.x = 60;
          scene.add(toAdd);
      }
    }

svg pathからthree.jsへの変換は下記のライブラリを使用します。

d3-threeD
https://github.com/asutherland/d3-threeD

このライブラリにはtransformSVGPathという内部関数があります。この関数がsvg pathからthree.jsのジオメトリに変換します。

この関数は内部関数のため、d3-threeD.jsを以下のように変更して外部から利用できるようにします。

d3-threeD.js

/* This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this file,
 * You can obtain one at http://mozilla.org/MPL/2.0/. */
var transformSVGPathExposed;

function d3threeD(exports) {
// 略
transformSVGPathExposed = transformSVGPath;
}

地域メッシュをBoxGeometryで表現する

地域メッシュのGeojsonからBoxGeometryを作成します。
geojsonの経度緯度をBoxGeometryのx,yに指定し、男、女の人口数をz軸に指定しています。

    function createPopulationMesh(geodata) {
      //geoJSON→svg path→three.js mesh変換    
      var countries = [];
      for (i = 0 ; i < geodata.length ; i++) {
        var geoFeature = geodata[i];
        console.log();
        var properties = geoFeature.properties;

        var extX = d3.extent(geoFeature.geometry.coordinates[0], function(d) { return d[0];});
        var extY = d3.extent(geoFeature.geometry.coordinates[0], function(d) { return d[1];});
        var ptMin = reverseProjection([extX[0], extY[0]]);
        var ptMax = reverseProjection([extX[1], extY[1]]);

        var vm = parseFloat(geoFeature.properties['男']) /10000;
        var vw = parseFloat(geoFeature.properties['女']) /10000;

        var geometry = new THREE.BoxGeometry(Math.abs(ptMax[0] - ptMin[0]), Math.abs(ptMax[1] - ptMin[1]), vm);
        var material = new THREE.MeshLambertMaterial({
          color:0x0000ff,
          transparent: true,
          opacity: 0.7
        });
        var mesh = new THREE.Mesh(geometry, material);
        mesh.position.set(ptMin[0], ptMin[1], vm/2);
        scene.add(mesh);

        var geometry = new THREE.BoxGeometry(Math.abs(ptMax[0] - ptMin[0]), Math.abs(ptMax[1] - ptMin[1]), vw);
        var material = new THREE.MeshLambertMaterial({
          color:0xff0000,
          transparent: true,
          opacity: 0.7
        });
        var mesh = new THREE.Mesh(geometry, material);
        mesh.position.set(ptMin[0], ptMin[1], vm + (vw/2));
        scene.add(mesh);
      }
    }

まとめ

three.jsを使用すれば簡単に3D表現ができます。
D3.jsを使用してgeojsonをSVGに簡単に変換できます。
D3.jsで作成したSVGをthree.jsで使用できるようにするには、d3-threeD.jsの内部関数、transformSVGPathを使用します。

地域メッシュは3Dで表現するより色で表現した方がわかりやすいってはっきりわかんだね。

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

国土数値情報の鉄道データで取得できる鉄道路線を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を使用して標高を求める

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

そこで、経度緯度から標高を取得できるようにします。
これには、地理院地図の標高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
      }
    },// 略
  ]
} 

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

クライアント側の処理は大きく次の通りです。
1.取得したGeoJSONのz軸を付与する。この際、開始点と終了点はpropertiesで指定した標高、それ以外は、開始点と終了点を考慮した割合とする。
2.結合できる線は結合する
3.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で表現できます。でも地下鉄とか新幹線は、対応できないこのプログラムはできそこないだ、使えないよ!

流体とか軟体を扱えるLiquidFunをJavaScriptで100%利用する

LiquidFunは流体とか軟体を扱うためのライブラリで、Box2dをベースに実装しています。
基本的にC++で実装されていますが、JavaまたはJavaScriptからも使用できます。
今回はJavaScriptでLiquidFunを利用する方法について説明します。

Liquidで表現したゆっくり
yukkuri.gif
http://needtec.sakura.ne.jp/box2d_yunyaa/yukkuri.html

LiquidFun
http://google.github.io/liquidfun/

導入方法とサンプルコード

以下からliquidfun-x.x.x.zipをダウンロードして任意のフォルダに解凍します。
https://github.com/google/liquidfun/releases

以下のフォルダがJavaScript関係のコードになります。
liquidfun-1.1.0\liquidfun\Box2D\lfjs

色々なコードが存在しますが、使うだけなら、「liquidfun.js」を以下のように取り込むだけで使えます。

<script type="text/javascript" src="liquidfun-1.1.0/liquidfun/Box2D/lfjs/liquidfun.js"></script>

index.htmlならびにtestbed以下はサンプルコードになっています。
LiquidFunはあくまで物理エンジンのライブラリであり、描画はおこなっていません。また、C++と異なり、JavaScript版にはデバッグ用の表示機能も有効になっていません。サンプルコードでは、three.jsを用いて描画を行っています。

もし、three.jsを使う学習コストが高いようでしたら、下記のようにD3.jsで描画するといいでしょう。

LiquidFunでJavaScript流体シミュレーション.
http://qiita.com/Quramy/items/578efec667267acf6871

サンプルコードの注意

もし、LiquidFunに付属するthree.jsでなく、最新のthree.jsを用いてサンプルを動作させる場合には注意があります。
renderer.js中のaddAttributeの呼び出し方法が古いので以下のように修正する必要があります。

liquidfun-1.1.0/liquidfun/Box2D/lfjs/testbed/renderer.js

function Renderer() {
  // init large buffer geometry
  this.maxVertices = 31000;
  var geometry = new THREE.BufferGeometry();
  geometry.dynamic = true;
  //addAttributeの形式が古い
  //geometry.addAttribute('position', Float32Array, this.maxVertices, 3);
  //geometry.addAttribute('color', Float32Array, this.maxVertices, 3);
  geometry.addAttribute('position', new THREE.BufferAttribute(new Float32Array(this.maxVertices * 3 ), 3));
  geometry.addAttribute('color', new THREE.BufferAttribute(new Float32Array(this.maxVertices * 3), 3));

JavaScriptの拡張

使っているとわかるのですが、JavaScript版のLiquidFunはC++版の一部の機能しかサポートされていません。

たとえば、粒子と剛体の接触イベントはC++ではサポートされています。
http://google.github.io/liquidfun/API-Ref/html/classb2_contact_listener.html

しかしながら、これらの機能は、JavaScriptではサポートされていません。
JavaScriptでは剛体と剛体の接触しかイベントがコールされないのです。

この場合、以下の手順でC++のコードをJavaScriptから利用できるようにする必要があります。

作成までの流れ.png

まず、emscriptenを用いてC++のコードをJavaScriptのコードに変換します。
次に、複数のJavaScriptのファイルをClosureCompilerを用いて、1つのファイルに圧縮します。

このビルドについての公式資料は下記を参考にしてください。
https://google.github.io/liquidfun/Building/html/md__building_java_script.html

ここでは、実際にDebian7上でliquidfunのJavaScriptをビルドしてみます。

C++のliquidfunをビルドする

まず、C++版のliquidfunのコードをビルドします。
C++でビルドするためには下記のライブラリーがインストール済みである必要があります。

apt-get install cmake
apt-get install libglapi-mesa
apt-get install libglu1-mesa-dev
apt-get install libxi-dev

その後は次のようにコマンドを実行してください。

cd liquidfun/Box2D
cmake -G'Unix Makefiles'
make
make install

emscriptenの前準備

emscriptenを動作させるのにnode.jsとPythonが必要なのでインストールします。
PythonはLinux系ならデフォルトでインストールされているので、それを利用するといいでしょう。

node.jsは下記の方法でインストールできます。

Debian/UbuntuでNode.jsをインストールする(nvm)
http://qiita.com/tamurashingo@github/items/6348863668e1e3fd70c9

emscriptenの構築

emscriptenはCやC++のコードをJavaScriptに変換します。
たとえば、以下で紹介されているEM-DOSBOXはEmScriptenでMS-DOSをJavascript上で動作させています。

Internet Archiveが2000以上のMS-DOSゲームをWebに移植、ブラウザで遊べる
http://jp.techcrunch.com/2015/01/10/20150109internet-archive-brings-oregon-trail-prince-of-persia-lemmings-and-2200-other-ms-dos-games-to-your-browser/

emscriptenは一応クロスプラットフォーム対応なので、Windowsでも動作しますが、 emsdk-1.29.0のインストーラーでインストールしたら環境変数をぶち壊してシステム復元をするはめになった ので、おすすめしません。
Windowsしかもってない方も、VMPlayer+Debian7で仮想環境を作れるのでそっちをお勧めします。

そもそも、liquidfunの公式で「We use Emscripten on Linux, but you should be able to use the Emscripten SDK on Mac or Windows too, if you prefer. Note that Mac and Windows build environments have not been tested.」って言っているので、おとなしくLinux系でやったほうがいいでしょう。

まず、下記のページから「Portable Emscripten SDK for Linux and OS X (emsdk-portable.tar.gz)」をダウンロードします。
http://kripken.github.io/emscripten-site/docs/getting_started/downloads.html

圧縮ファイルを取得後は次のコマンドを実行すればEmscriptenのビルドが行われます。

gzip -dc emsdk-portable.tar.gz | tar xvf -
cd emsdk_portable
./emsdk update
./emsdk install latest

この際に、clang/fastcompをダウンロードしてビルドをしていますが、この処理にはかなり時間がかかります。また、このビルドを中断して、再実行した場合は注意が必要です。
emsdk installでは、フォルダの有無でインストール済みかを判断しているので、中途半端に処理を行っているとインストール済みとみなされます。
その場合は手動でfastcomp/build_masterをmakeしなおしましょう。

ビルドが完了したら、以下のコマンドを実行して使用するemsdkのバージョンと環境変数の登録を行います。

./emsdk activate latest
source ./emsdk_env.sh

ClosureCompilerの構築

ClosureComilerはJavascriptの圧縮を行うJavaで作成されたコマンドラインツールです。

まず、以下から compiler.jarをダウンロードしてください。
https://developers.google.com/closure/compiler/docs/gettingstarted_app

その後、環境変数でCLOSURE_JARにcomiler.jarのパスを指定します。

export CLOSURE_JAR=/share/closure/compiler.jar

Javaのバージョン確認

Javaのバージョンが7未満の場合、comiler.jarは以下のようなエラーを出力します。

Exception in thread "main" java.lang.UnsupportedClassVersionError: com/google/javascript/jscomp/CommandLineRunner : Unsupported major.minor version 51.0

この場合、Javaのバージョンを確認してみてください。

java --version 

DebianのJavaを7にあげるには以下の手順で行うと良いでしょう。

Installing Java 7 (Oracle) in Debian via apt-get [closed]
http://stackoverflow.com/questions/15543603/installing-java-7-oracle-in-debian-via-apt-get

echo "deb http://ppa.launchpad.net/webupd8team/java/ubuntu precise main" | tee -a /etc/apt/sources.list
echo "deb-src http://ppa.launchpad.net/webupd8team/java/ubuntu precise main" | tee -a /etc/apt/sources.list
apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys EEA14886
apt-get update
apt-get install oracle-java7-installer

liquidfunのJavaScriptをビルドする

liquidfunのJavaScriptをビルドします。
公式では下記のコマンドを実行するだけで作成できるといっていますし、実際作成できます。

cd Box2D/lfjs
make
./uglify.sh

しかし、ここで作成された、liquidfun.jsを使用するとエラーになります。
この場合は、Makefileを修正します。

$(EMSCRIPTEN)/emcc -I../ -o lf_core.js jsBindings/jsBindings.cpp $(OBJECTS) -s $(EXPORTS) -s TOTAL_MEMORY=33554432 -s ASSERTIONS=1 -O1 --js-library callbacks.js

ここでは、-O2から-O1に変更しました。また、-s ASSERTIONS=1とすることで、エラーの詳細を出力できます。-O2の場合、以下のエラーが発生していました。

"Assertion failed: you need to wait for the runtime to be ready (e.g. wait for main() to be called)"

emscriptenは開発中のようで、安定しておらず、同様の自称が下記で報告されています。
https://github.com/kripken/emscripten/issues/3082

とりあえず、ここまでで、liquidfun.jsのビルドが行えるようになりました。
既存のプログラムを新しいliquidfun.jsで動作させて動くことを確認しましょう。

拡張のコツ

ビルドまでも一苦労でしたが、実際、拡張するのも結構しんどいです。

lfjs/jsBindings以下のコードをよく読んで、既存のオブジェクトがどのようにJavaScriptで実装されているかつかみましょう。

offsets.jsを使用したクラス・構造体へのプロパティのアクセス方法

jsBindings/offsets.js には各オブジェクトの構造体やクラスのプロパティについてのオフセット値が格納されています。ここで指定したオフセットはプロパティのアクセスに利用します。

たとえば、b2ParticleContact構造体を実装する場合を例に挙げましょう。
b2ParticleContactの構造体はC++では次のようになります。

struct b2ParticleContact
{
// 略
private:
    b2ParticleIndex indexA, indexB;

    float32 weight;

    b2Vec2 normal;

    uint32 flags;
};
// 略

この構造体のweightを取得するJavaScriptコードは次のようになります。

lfjs/jsBindings/Particle/b2ParticleSystem.js

/**@constructor*/
function b2ParticleContact(ptr) {
  this.ptr = ptr;
  this.buffer = new DataView(Module.HEAPU8.buffer, ptr);
}

b2ParticleContact.prototype.GetWeight = function() {
  return this.buffer.getFloat32(Offsets.b2ParticleContact.weight, true);
};

このようにoffset.jsで指定されたoffset値とオブジェクト作成時に与えられたポインターを利用して構造体のメンバーにアクセスしています。

では、offset.jsのオフセット値はどのように指定するのでしょうか?
実は、このオフセット値を指定するためのヒントが、lfjs/jsBindings/jsBindings.cppにあります。

まず、lfjs/jsBindings/jsBindings.cppには次のようにコメントアウトされた行があるのでコレを有効にします。

lfjs/jsBindings/jsBindings.cpp

// TODO clean all of this up, and/or make it auto generated from the
// header files
void PrintOffsets(b2Body* b) {
  printf("\tb2Body: {\n");
  printf("\t\ttype: %u,\n", (unsigned int)&b->m_type - (unsigned int)b);
  printf("\t\tislandIndex: %u,\n", (unsigned int)&b->m_islandIndex - (unsigned int)b);
  printf("\t\txf: %u,\n", (unsigned int)&b->m_xf - (unsigned int)b);
  printf("\t\txf0: %u,\n", (unsigned int)&b->m_xf0 - (unsigned int)b);
  printf("\t\tsweep: %u,\n", (unsigned int)&b->m_sweep - (unsigned int)b);
  printf("\t\tlinearVelocity: %u,\n", (unsigned int)&b->m_linearVelocity - (unsigned int)b);
  printf("\t\tangularVelocity: %u,\n", (unsigned int)&b->m_angularVelocity - (unsigned int)b);
  printf("\t\tforce: %u,\n", (unsigned int)&b->m_force - (unsigned int)b);
  printf("\t\ttorque: %u,\n", (unsigned int)&b->m_torque - (unsigned int)b);
  printf("\t\tworld: %u,\n", (unsigned int)&b->m_world - (unsigned int)b);
  printf("\t\tprev: %u,\n", (unsigned int)&b->m_prev - (unsigned int)b);
  printf("\t\tnext: %u,\n", (unsigned int)&b->m_next - (unsigned int)b);
  printf("\t\tfixtureList: %u,\n", (unsigned int)&b->m_fixtureList - (unsigned int)b);
  printf("\t\tfixtureCount: %u,\n", (unsigned int)&b->m_fixtureCount - (unsigned int)b);
  printf("\t\tjointList: %u,\n", (unsigned int)&b->m_jointList - (unsigned int)b);
  printf("\t\tcontactList: %u,\n", (unsigned int)&b->m_contactList - (unsigned int)b);
  printf("\t\tmass: %u,\n", (unsigned int)&b->m_mass - (unsigned int)b);
  printf("\t\tinvMass: %u,\n", (unsigned int)&b->m_invMass - (unsigned int)b);
  printf("\t\tI: %u,\n", (unsigned int)&b->m_I - (unsigned int)b);
  printf("\t\tinvI: %u,\n", (unsigned int)&b->m_invI - (unsigned int)b);
  printf("\t\tlinearDamping: %u,\n", (unsigned int)&b->m_linearDamping - (unsigned int)b);
  printf("\t\tangularDamping: %u,\n", (unsigned int)&b->m_angularDamping - (unsigned int)b);
  printf("\t\tgravityScale: %u,\n", (unsigned int)&b->m_gravityScale - (unsigned int)b);
  printf("\t\tsleepTime: %u,\n", (unsigned int)&b->m_sleepTime - (unsigned int)b);
  printf("\t\tuserData: %u\n", (unsigned int)&b->m_userData - (unsigned int)b);
  printf("\t}\n");
}

void PrintOffsets(b2Contact* c) {
  printf("\tb2Contact: {\n");
  printf("\t\tflags: %u,\n", (unsigned int)&c->m_flags - (unsigned int)c);
  printf("\t\tprev: %u,\n", (unsigned int)&c->m_prev - (unsigned int)c);
  printf("\t\tnext: %u,\n", (unsigned int)&c->m_next - (unsigned int)c);
  printf("\t\tnodeA: %u,\n", (unsigned int)&c->m_nodeA - (unsigned int)c);
  printf("\t\tnodeB: %u,\n", (unsigned int)&c->m_nodeB - (unsigned int)c);
  printf("\t\tfixtureA: %u,\n", (unsigned int)&c->m_fixtureA - (unsigned int)c);
  printf("\t\tfixtureB: %u,\n", (unsigned int)&c->m_fixtureB - (unsigned int)c);
  printf("\t\tindexA: %u,\n", (unsigned int)&c->m_indexA - (unsigned int)c);
  printf("\t\tindexB: %u,\n", (unsigned int)&c->m_indexB - (unsigned int)c);
  printf("\t\tmanifold: %u,\n", (unsigned int)&c->m_manifold - (unsigned int)c);
  printf("\t\ttoiCount: %u,\n", (unsigned int)&c->m_toiCount - (unsigned int)c);
  printf("\t\ttoi: %u,\n", (unsigned int)&c->m_toi - (unsigned int)c);
  printf("\t\tfriction: %u,\n", (unsigned int)&c->m_friction - (unsigned int)c);
  printf("\t\trestitution: %u,\n", (unsigned int)&c->m_restitution - (unsigned int)c);
  printf("\t\ttangentSpeed: %u\n", (unsigned int)&c->m_tangentSpeed - (unsigned int)c);
  printf("\t},\n");
}

void PrintOffsets(b2Fixture* f) {
  printf("\tb2Fixture: {\n");
  printf("\t\tdensity: %u,\n", (unsigned int)&f->m_density - (unsigned int)f);
  printf("\t\tnext: %u,\n", (unsigned int)&f->m_next - (unsigned int)f);
  printf("\t\tbody: %u,\n", (unsigned int)&f->m_body - (unsigned int)f);
  printf("\t\tshape: %u,\n", (unsigned int)&f->m_shape - (unsigned int)f);
  printf("\t\tfriction: %u,\n", (unsigned int)&f->m_friction - (unsigned int)f);
  printf("\t\trestitution: %u,\n", (unsigned int)&f->m_restitution - (unsigned int)f);
  printf("\t\tproxies: %u,\n", (unsigned int)&f->m_proxies - (unsigned int)f);
  printf("\t\tproxyCount: %u,\n", (unsigned int)&f->m_proxyCount - (unsigned int)f);
  printf("\t\tfilter: %u,\n", (unsigned int)&f->m_filter - (unsigned int)f);
  printf("\t\tisSensor: %u,\n", (unsigned int)&f->m_isSensor - (unsigned int)f);
  printf("\t\tuserData: %u\n", (unsigned int)&f->m_userData - (unsigned int)f);
  printf("\t},\n");
}

void PrintOffsets(b2ParticleGroup* p) {
  printf("\tb2ParticleGroup: {\n");
  printf("\t\tsystem: %u,\n", (unsigned int)&p->m_system - (unsigned int)p);
  printf("\t\tfirstIndex: %u,\n", (unsigned int)&p->m_firstIndex - (unsigned int)p);
  printf("\t\tlastIndex: %u,\n", (unsigned int)&p->m_lastIndex - (unsigned int)p);
  printf("\t\tgroupFlags: %u,\n", (unsigned int)&p->m_groupFlags - (unsigned int)p);
  printf("\t\tstrength: %u,\n", (unsigned int)&p->m_strength - (unsigned int)p);
  printf("\t\tprev: %u,\n", (unsigned int)&p->m_prev - (unsigned int)p);
  printf("\t\tnext: %u,\n", (unsigned int)&p->m_next - (unsigned int)p);
  printf("\t\ttimestamp: %u,\n", (unsigned int)&p->m_timestamp - (unsigned int)p);
  printf("\t\tmass: %u,\n", (unsigned int)&p->m_mass - (unsigned int)p);
  printf("\t\tinertia: %u,\n", (unsigned int)&p->m_inertia - (unsigned int)p);
  printf("\t\tcenter: %u,\n", (unsigned int)&p->m_center - (unsigned int)p);
  printf("\t\tlinearVelocity: %u,\n", (unsigned int)&p->m_linearVelocity - (unsigned int)p);
  printf("\t\tangularVelocity: %u,\n", (unsigned int)&p->m_angularVelocity - (unsigned int)p);
  printf("\t\ttransform: %u,\n", (unsigned int)&p->m_transform - (unsigned int)p);
  printf("\t\tuserData: %u\n", (unsigned int)&p->m_userData - (unsigned int)p);
  printf("\t},\n");
}

void PrintOffsets(b2World* w) {
  printf("\tb2World: {\n");
  printf("\t\tbodyList: %u\n", (unsigned int)&w->m_bodyList - (unsigned int)w);
  printf("\t},\n");
}

void PrintOffsets(b2WorldManifold* wm) {
  printf("\tb2WorldManifold: {\n");
  printf("\t\tnormal: %u,\n", (unsigned int)&wm->normal - (unsigned int)wm);
  printf("\t\tpoints: %u,\n", (unsigned int)&wm->points - (unsigned int)wm);
  printf("\t\tseparations: %u\n", (unsigned int)&wm->separations - (unsigned int)wm);
  printf("\t},\n");
}

void PrintOffsets(b2ParticleContact* pc) {
  printf("\tb2ParticleContact: {\n");
  printf("\t\tindexA: %u,\n", (unsigned int)&pc->indexA - (unsigned int)pc);
  printf("\t\tindexB: %u,\n", (unsigned int)&pc->indexB - (unsigned int)pc);
  printf("\t\tweight: %u,\n", (unsigned int)&pc->weight - (unsigned int)pc);
  printf("\t\tnormal: %u,\n", (unsigned int)&pc->normal - (unsigned int)pc);
  printf("\t\tflags: %u\n", (unsigned int)&pc->flags - (unsigned int)pc);
  printf("\t},\n");
}

void PrintOffsets(b2ParticleBodyContact* pbc) {
  printf("\tb2ParticleBodyContact: {\n");
  printf("\t\tindex: %u,\n", (unsigned int)&pbc->index - (unsigned int)pbc);
  printf("\t\tbody: %u,\n", (unsigned int)&pbc->body - (unsigned int)pbc);
  printf("\t\tfixture: %u,\n", (unsigned int)&pbc->fixture - (unsigned int)pbc);
  printf("\t\tweight: %u,\n", (unsigned int)&pbc->weight - (unsigned int)pbc);
  printf("\t\tnormal: %u,\n", (unsigned int)&pbc->normal - (unsigned int)pbc);
  printf("\t\tmass: %u\n", (unsigned int)&pbc->mass - (unsigned int)pbc);
  printf("\t},\n");
}

# include <stdio.h>
extern "C" {
void GenerateOffsets() {
  printf("{\n");
  //create dummy body to generate offsets
  b2World world(b2Vec2(0, 0));
  b2BodyDef def;
  b2Body* body1 = world.CreateBody(&def);

  b2CircleShape s;
  b2FixtureDef d;
  d.shape = &s;
  b2Fixture* f =body1->CreateFixture(&d);

  b2ParticleSystemDef psd;
  b2ParticleSystem* ps = world.CreateParticleSystem(&psd);

  b2ParticleGroupDef pgd;
  b2ParticleGroup* pg = ps->CreateParticleGroup(pgd);

  b2WorldManifold worldManifold;
  PrintOffsets(body1);
  PrintOffsets(f);
  PrintOffsets(pg);
  PrintOffsets(&world);
  PrintOffsets(&worldManifold);

  b2ParticleContact pc;
  PrintOffsets(&pc);

  b2ParticleBodyContact pbc;
  PrintOffsets(&pbc);

  // need to instantiate contact differently
  //b2Contact contact;
  //PrintOffsets(&contact);

  printf("};\n");
}

これはGenerateOffsetsを経由して、オフセット値を出力するための関数です。
liquidfun.jsビルド後、Javascript側で下記のコードを実行すると、offset.jsに記述すべきオフセット値が出力されます。

console.log(_GenerateOffsets());

・・・ってことをしたかたんでしょうが、実はこのコードはビルドすらできません。
それは当然で、liquidfunのC++のコードのプライベートのメンバーにアクセスしているからです。この場合、ビルド時にエラーが発生した箇所をかたっぱしからpublicに指定して、offset.jsを作成後元に戻す必要があります。

粒子の衝突イベントのサポート

particleとparticleの衝突、particleとBodyの衝突イベントをサポートしたバージョンを下記に置きます。
http://needtec.sakura.ne.jp/box2d_yunyaa/lfjs.zip

Winmergeなりで、なにを替えたかを確認すれば、さらなる拡張のヒントになるかとおもいます。

このイベントを利用するには、particleのflagsにb2_particleContactListenerParticle とb2_fixtureContactListenerParticleを付与します。

    var psd = new b2ParticleSystemDef();
    console.log(psd);
    psd.radius = particleRadius;
    this.particleSystem = world.CreateParticleSystem(psd);
    circle = new b2PolygonShape();
    circle.SetAsBoxXYCenterAngle(r, r, new b2Vec2(x, y), 0);
    pgd = new b2ParticleGroupDef();
    pgd.flags = b2_elasticParticle  | b2_particleContactListenerParticle | b2_fixtureContactListenerParticle ;
    pgd.groupFlags = b2_solidParticleGroup;
    pgd.shape = circle;
    this.particleSystem.CreateParticleGroup(pgd);

その後、以下のようにコールバック関数を登録することで、該当フラグを付与したオブジェクトの衝突が検知できます。

    world = new b2World(gravity);
    var listener =  {
        BeginContactBody: function(contact) {
          console.log(contact.GetFixtureA());
        },
        EndContactBody: function(contact) {
            console.log(contact.GetFixtureA());
        },
        PostSolve: function(contact, impulse) {
            console.log('PostSolve');
        },
        PreSolve: function(contact, oldManifold) {

        },
        /**
         * ParticleオブジェクトとBodyの接触
         */
        BeginContactParticleBody : function(particleSystem, particleBodyContact) {
          console.log('----------------BeginContactParticleBody', 
            particleSystem,
            particleBodyContact.GetIndex(),
            particleBodyContact.GetNormal(),
            particleBodyContact.GetWeight(),
            particleBodyContact.GetMass(),
            particleBodyContact.GetBody(),
            particleBodyContact.GetFixture());
        },
        /**
         * Particleオブジェクト同士の接触
         */
        BeginContactParticle : function(particleSystem, particleContact) {
          console.log('----------------BeginContactParticleContact', particleSystem, particleContact.GetFlags(), particleContact.GetIndexA(), particleContact.GetIndexB(), particleContact.GetNormal(), particleContact.GetWeight());
        },
        EndContactParticleBody : function(fixture, particleSystem, index) {
          console.log('----------------EndContactParticleBody', fixture, particleSystem, index);
        },
        EndContactParticle : function(particleSystem, indexA, indexB) {
          console.log('----------------EndContactParticle', indexA, indexB);
        }
    }
    world.SetContactListener(listener);

以下のサンプルのConsoleLogを監視してみてください。
「ゆっくり」同士をぶつけた場合と離れた場合にメッセージが表示することが確認できます。

http://needtec.sakura.ne.jp/box2d_yunyaa/yukkuri.html

まとめ

LiquidFunを利用すると流体の物理演算をおこなってくれるので、流体やSoftBodyの表現がおこなえます。

基本的にC++で100%の実力を発揮できますが、頑張ればJavaScriptでもなんとかなります。
・・・とはいえ、くっそメンドクサイのです。

政府統計窓口(eStat)の地域メッシュをWebブラウザで表示する方法

概要

地域メッシュ統計とは,緯度・経度に基づき地域を隙間なく網の目(Mesh)の区域に分けて,統計データをそれぞれの区域に編成したものです。

メッシュ.png
※色が濃いほど人口が多い。

デモ:
http://needtec.sakura.ne.jp/estat/population
http://needtec.sakura.ne.jp/tokuraku/passenger.html

このデータの元になっている情報は政府統計窓口(estat)から取得しています。
以前はPythonでmatplotlib で描画していましたが、ブラウザで表示するようにしました。

政府統計の総合窓口(e-Stat)のAPIを使ってみよう
http://qiita.com/mima_ita/items/44f358dc1bc4000d365d

なお、地域メッシュの詳しい求め方は以下を参考にしてください。
http://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf

データ描画までの流れ

描画のたびに、政府統計窓口に問い合わせるっていうのは、遅すぎるのでデータ描画までの手順を以下のように分けます。

1.ある統計に紐づく地域メッシュ情報を政府統計窓口より、取得する
2.spatialiteに記録する
3.Webサーバーでリクエストに応じてgeojson形式でメッシュ情報を返すようにする。
4.画面側は必要な領域だけWebサーバーに要求して取得したgeojsonを元に色を塗る。

これらの手順をサーバーサイドはPython、クライアントサイドはJavaScriptで実施します。

実際のコードは下記から取得できます。
https://github.com/mima3/estat

spatialiteについての説明は下記を参考にしてください。
地図とかの空間情報をSQLiteに格納するSpatiaLiteを使用してみる
http://qiita.com/mima_ita/items/64f6c2b8bb47c4b5b391

地域メッシュを格納するデータベースの作成

ここで紹介する統計情報を格納するデータベースの実装は以下のコードになります。
https://github.com/mima3/estat/blob/master/estat_db.py

テーブル構成

estatの地域メッシュを格納するためのテーブルは以下のようになります。

database.png

Table名:Stat
統計に対する情報を格納するテーブル

列名 説明
stat_id TextField 主キー。統計情報のID
stat_name TextField 政府統計名
stat_name_code TextField 政府統計名コード
gov_org TextField 作成機関名
gov_org_code TextField 作成機関名コード
survey_date TextField 調査日
title TextField 統計表の表題

Table名:StatValue
各統計の値を格納する

列名 説明
id PrimaryKeyField 自動インクリメントのID
stat_id TextField 統計情報のID
value TextField 統計の値

Table名:StatValueAttr
各統計の値に紐づく属性値を格納する

列名 説明
id PrimaryKeyField 自動インクリメントのID
stat_id TextField 統計情報のID
stat_value_id INT StatValueのidの外部キー
attr_name TextField 属性名
attr_value TextField 属性値

Table名:MapArea
属性areaにおけるPolygonを格納する

列名 説明
id PrimaryKeyField 自動インクリメントのID
stat_id TextField 統計情報のID
stat_val_attr_id INT StatValueAttrのID
geometry Polygon ジオメトリ情報

Pythonでテーブルの作成方法

pythonでpeeweeを用いて作成する場合の実際のコードは以下のようになります。

estat_db.py

# -*- coding: utf-8 -*-
import os
from peewee import *
from playhouse.sqlite_ext import SqliteExtDatabase
from estat_go_jp import *
import jpgrid

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

SRID = 4326

class PolygonField(Field):
    db_field = 'polygon'

class PointField(Field):
    db_field = 'point'

class LineStringField(Field):
    db_field = 'linestring'

class MultiPolygonField(Field):
    db_field = 'multipolygon'

class MultiPointField(Field):
    db_field = 'multipoint'

class MultiLineStringField(Field):
    db_field = 'multilinestring'

class Stat(Model):
    """
    統計情報を格納するモデル
    """
    stat_id = TextField(primary_key=True)
    stat_name = TextField()
    stat_name_code = TextField()
    gov_org = TextField()
    gov_org_code = TextField()
    #statistics_name = TextField()
    #cycle = TextField()
    survey_date = TextField()
    #open_date = TextField()
    #small_area = TextField()
    title = TextField()

    class Meta:
        database = database_proxy

class StatValue(Model):
    """
    統計の値を格納するモデル
    """
    id = PrimaryKeyField()
    stat_id = TextField(index=True, unique=False)
    value = TextField()

    class Meta:
        database = database_proxy

class StatValueAttr(Model):
    """
    統計の属性値を格納するモデル
    """
    id = PrimaryKeyField()
    stat_id = TextField()
    stat_value = ForeignKeyField(
        db_column='stat_value_id',
        rel_model=StatValue,
        to_field='id'
    )
    attr_name = TextField()
    attr_value = TextField()

    class Meta:
        database = database_proxy
        indexes = (
            (('stat_id', 'stat_value'), False),
            (('stat_id', 'stat_value', 'attr_name'), False),
        )

class MapArea(Model):
    """
    統計の属性値 areaにおけるPolygonを格納するモデル
    """
    id = PrimaryKeyField()
    stat_id = TextField()
    stat_val_attr = ForeignKeyField(
        db_column='stat_val_attr_id',
        rel_model=StatValueAttr,
        to_field='id'
    )
    geometry = PolygonField()

    class Meta:
        database = database_proxy
        indexes = (
            (('stat_id', 'stat_val_attr'), True),
        )

class idx_MapArea_Geometry(Model):
    pkid = PrimaryKeyField()
    xmin = FloatField()
    xmax = FloatField()
    ymin = FloatField()
    ymax = FloatField()

    class Meta:
        database = database_proxy

def connect(path, spatialite_path, evn_sep=';'):
    """
    データベースへの接続
    @param path sqliteのパス
    @param spatialite_path mod_spatialiteへのパス
    @param env_sep 環境変数PATHの接続文字 WINDOWSは; LINUXは:
    """
    os.environ["PATH"] = os.environ["PATH"] + evn_sep + os.path.dirname(spatialite_path)
    db = SqliteExtDatabase(path)
    database_proxy.initialize(db)
    db.field_overrides = {
        'polygon': 'POLYGON',
        'point': 'POINT',
        'linestring': 'LINESTRING',
        'multipolygon': 'MULTIPOLYGON',
        'multipoint': 'MULTIPOINT',
        'multilinestring': 'MULTILINESTRING',
    }
    db.load_extension(os.path.basename(spatialite_path))

def setup(path, spatialite_path, evn_sep=';'):
    """
    データベースの作成
    @param path sqliteのパス
    @param spatialite_path mod_spatialiteへのパス
    @param env_sep 環境変数PATHの接続文字 WINDOWSは; LINUXは:
    """
    connect(path, spatialite_path, evn_sep)

    database_proxy.create_tables([Stat, StatValue, StatValueAttr], True)
    database_proxy.get_conn().execute('SELECT InitSpatialMetaData()')

    # Geometryのテーブルは直接実装する必要がある.
    database_proxy.get_conn().execute("""
        CREATE TABLE IF NOT EXISTS "MapArea" (
          "id" INTEGER PRIMARY KEY AUTOINCREMENT,
          "stat_id" TEXT,
          "stat_val_attr_id" INTEGER ,
          FOREIGN KEY(stat_val_attr_id) REFERENCES StatValueAttr(id));
    """)

    database_proxy.get_conn().execute("""
        CREATE INDEX IF NOT EXISTS "ix_MapArea_stat_id" ON MapArea(stat_id);
    """)

    database_proxy.get_conn().execute("""
        CREATE INDEX IF NOT EXISTS "ix_MapArea_stat_id_stat_val_attr_id" ON MapArea(stat_id, stat_val_attr_id);
    """)

    database_proxy.get_conn().execute("""
        Select AddGeometryColumn ("MapArea", "Geometry", ?, "POLYGON", 2);
    """, (SRID,))

    database_proxy.get_conn().execute("""
        SELECT CreateSpatialIndex("MapArea", "geometry")
    """)

setup関数でテーブルの作成を行っています。
ジオメトリを含まない場合、peeweeでは以下のコードでテーブルを生成します。

database_proxy.create_tables([Stat, StatValue, StatValueAttr], True)

しかし、ジオメトリーの列を含む場合、以下のように、自分で実装する必要があります

    # Geometryのテーブルは直接実装する必要がある.
    database_proxy.get_conn().execute("""
        CREATE TABLE IF NOT EXISTS "MapArea" (
          "id" INTEGER PRIMARY KEY AUTOINCREMENT,
          "stat_id" TEXT,
          "stat_val_attr_id" INTEGER ,
          FOREIGN KEY(stat_val_attr_id) REFERENCES StatValueAttr(id));
    """)

    database_proxy.get_conn().execute("""
        CREATE INDEX IF NOT EXISTS "ix_MapArea_stat_id" ON MapArea(stat_id);
    """)

    database_proxy.get_conn().execute("""
        CREATE INDEX IF NOT EXISTS "ix_MapArea_stat_id_stat_val_attr_id" ON MapArea(stat_id, stat_val_attr_id);
    """)

    database_proxy.get_conn().execute("""
        Select AddGeometryColumn ("MapArea", "Geometry", ?, "POLYGON", 2);
    """, (SRID,))

    database_proxy.get_conn().execute("""
        SELECT CreateSpatialIndex("MapArea", "geometry")
    """)

このコードではジオメトリー以外の列を含んだテーブルを作成したあとに、「AddGeometryColumn」でジオメトリーの列を追加し、それに対し「CreateSpatialIndex」でインデックスを作成します。

ジオメトリーに対して作られたインデックスは仮想テーブルですので、Peeweeから使用することができません。インデックスを利用する際は、自分でSQLを記述する必要があります。

Pythonで統計情報をデータベースにインポートする方法

eStatの地域メッシュでは、一つの統計について、第一次メッシュごとに分割されて統計IDをもっています。
たとえば、「平成22年国勢調査-世界測地系(1KMメッシュ)20101001」の統計情報では、以下のような統計IDが取得できます。

・T000608M3622
・T000608M3623
・T000608M3624
・T000608M3653

このIDは「T000608」と、第一次メッシュコードから作成されています。

つまり、地域メッシュの統計情報を取得するには次のような手順が必要になります。

1.eStatのgetStatsList APIで該当の統計をすべて取得する。
2.eStatのgetStatsData APIで統計ID毎に値を取得する。
3.DBに格納する際に、メッシュコードから経度緯度の範囲をもとめて、Polygonとして格納する。
4.全ての統計IDに対して2~3を繰り返す

getStatesList APIによる該当の統計を全て取得する

「平成22年国勢調査-世界測地系(1KMメッシュ)20101001」に関係する統計IDをすべて取得するためのgetStatesList APIを実行するための関数は以下のように実装します。

参照元:
https://github.com/mima3/estat/blob/master/estat_go_jp.py

estat_go_jp.py

def get_stats_list(api_key, search_kind, key_word):
    """
    統計情報の検索
    @param api_key APIキー
    @param search_kind 統計の種類
    @param key_word 検索キー
    """
    key_word = urllib.quote(key_word.encode('utf-8'))
    url = ('http://api.e-stat.go.jp/rest/1.0/app/getStatsList?appId=%s&lang=J&searchKind=%s&searchWord=%s' % (api_key, search_kind, key_word))
    req = urllib2.Request(url)
    opener = urllib2.build_opener()
    conn = opener.open(req)
    cont = conn.read()
    parser = etree.XMLParser(recover=True)
    root = etree.fromstring(cont, parser)
    ret = []
    data_list = root.find('DATALIST_INF')
    list_infs = data_list.xpath('.//LIST_INF')
    for list_inf in list_infs:
        item = {
            'id': trim_data(list_inf.get('id'))
        }
        stat_name = list_inf.find('STAT_NAME')
        if stat_name is not None:
            item['stat_name'] = trim_data(stat_name.text)
            item['stat_name_code'] = trim_data(stat_name.get('code'))

        gov_org = list_inf.find('GOV_ORG')
        if gov_org is not None:
            item['gov_org'] = trim_data(gov_org.text)
            item['gov_org_code'] = trim_data(gov_org.get('code'))

        statistics_name = list_inf.find('STATISTICS_NAME')
        if statistics_name is not None:
            item['statistics_name'] = trim_data(statistics_name.text)

        title = list_inf.find('TITLE')
        if title is not None:
            item['title'] = trim_data(title.text)

        cycle = list_inf.find('CYCLE')
        if cycle is not None:
            item['cycle'] = cycle.text

        survey_date = list_inf.find('SURVEY_DATE')
        if survey_date is not None:
            item['survey_date'] = trim_data(survey_date.text)

        open_date = list_inf.find('OPEN_DATE')
        if open_date is not None:
            item['open_date'] = trim_data(open_date.text)

        small_area = list_inf.find('SMALL_AREA')
        if small_area is not None:
            item['small_area'] = trim_data(small_area.text)

        ret.append(item)
    return ret

基本的にurllib2で取得した内容をlxmlで解析してディレクトリに格納しているだけです。

getStatsData APIで統計ID毎に値を取得する。

各統計ID毎の値を取得するには以下のget_stats_id_value()を実行します。

参照元:
https://github.com/mima3/estat/blob/master/estat_go_jp.py

estat_go_jp.py

def get_meta_data(api_key, stats_data_id):
    """
    メタ情報取得
    """
    url = ('http://api.e-stat.go.jp/rest/1.0/app/getMetaInfo?appId=%s&lang=J&statsDataId=%s' % (api_key, stats_data_id))
    req = urllib2.Request(url)
    opener = urllib2.build_opener()
    conn = opener.open(req)
    cont = conn.read()
    parser = etree.XMLParser(recover=True)
    root = etree.fromstring(cont, parser)
    class_object_tags = root.xpath('//METADATA_INF/CLASS_INF/CLASS_OBJ')
    class_object = {}

    for class_object_tag in class_object_tags:
        class_object_id = class_object_tag.get('id')
        class_object_name = class_object_tag.get('name')
        class_object_item = {
            'id': trim_data(class_object_id),
            'name': trim_data(class_object_name),
            'objects': {}
        }
        class_tags = class_object_tag.xpath('.//CLASS')
        for class_tag in class_tags:
            class_item = {
                'code': trim_data(class_tag.get('code')),
                'name': trim_data(class_tag.get('name')),
                'level': trim_data(class_tag.get('level')),
                'unit': trim_data(class_tag.get('unit'))
            }
            class_object_item['objects'][class_item['code']] = class_item
        class_object[class_object_id] = class_object_item
    return class_object

def _get_stats_id_value(api_key, stats_data_id, class_object, start_position, filter_str):
    """
    統計情報の取得
    """
    url = ('http://api.e-stat.go.jp/rest/1.0/app/getStatsData?limit=10000&appId=%s&lang=J&statsDataId=%s&metaGetFlg=N&cntGetFlg=N%s' % (api_key, stats_data_id, filter_str))
    if start_position > 0:
        url = url + ('&startPosition=%d' % start_position)
    req = urllib2.Request(url)
    opener = urllib2.build_opener()
    conn = opener.open(req)
    cont = conn.read()
    parser = etree.XMLParser(recover=True)
    root = etree.fromstring(cont, parser)
    ret = []
    row = {}
    datas = {}
    value_tags = root.xpath('//STATISTICAL_DATA/DATA_INF/VALUE')
    for value_tag in value_tags:
        row = {}
        for key in class_object:
            val = value_tag.get(key)
            if val in class_object[key]['objects']:
                text = class_object[key]['objects'][val]['name']
                row[key] = trim_data(text.encode('utf-8'))
            else:
                row[key] = val.encode('utf-8')
        row['value'] = trim_data(value_tag.text)
        ret.append(row)
    return ret

def get_stats_id_value(api_key, stats_data_id, filter):
    """
    統計取得
    @param api_key APIキー
    @param stats_data_id 統計表ID
    @param filter_str フィルター文字
    """
    filter_str = ''
    for key in filter:
        filter_str += ('&%s=%s' % (key, urllib.quote(filter[key].encode('utf-8'))))
    class_object = get_meta_data(api_key, stats_data_id)
    return _get_stats_id_value(api_key, stats_data_id, class_object, 1, filter_str), class_object

get_stats_id_valueはフィルターを掛けて検索できますが、今回はフィルターなしで全て取得するので、そこは考慮しないでください。

手順としては、get_meta_dataにより、統計情報のメタ情報を取得します。これにより、該当の統計がどのような属性を有しているか調べます。

その後、_get_stats_id_value()を用いて、統計値と、その属性の値を取得します。

DBに格納する際に、メッシュコードから経度緯度の範囲をもとめて、Polygonとして格納する

eStatから取得するデータには残念なことに経度緯度の情報は存在しません。
そこで、メッシュコードから経度、緯度を取得する必要があります。
Pythonの場合、以下のライブラリを使用すると楽でしょう。

python-geohash
https://code.google.com/p/python-geohash/

このライブラリは以下のように、easy_install等でインストール可能です。

easy_install python-geohash

使い方はメッシュコードを指定するだけ、該当のメッシュコードの範囲が、経度と緯度で表示されます。

import jpgrid
jpgrid.bbox('305463') #メッシュコードを指定する → {'s': 154.375, 'e': 20.583333333333332, 'w': 20.5, 'n': 154.5}

これらを利用して、データベースに地域メッシュの統計情報を格納するコードは次のようになります。

参照元:
https://github.com/mima3/estat/blob/master/estat_db.py

estat_db.py

def import_stat(api_key, stat_id):
    """
    統計情報のインポート
    @param api_key e-statのAPIKEY
    @param stat_id 統計ID
    """
    with database_proxy.transaction():
        MapArea.delete().filter(MapArea.stat_id == stat_id).execute()
        StatValueAttr.delete().filter(StatValueAttr.stat_id == stat_id).execute()
        StatValue.delete().filter(StatValue.stat_id == stat_id).execute()
        Stat.delete().filter(Stat.stat_id == stat_id).execute()

        tableinf = get_table_inf(api_key, stat_id)
        stat_row = Stat.create(
            stat_id=stat_id,
            stat_name=tableinf['stat_name'],
            stat_name_code=tableinf['stat_name_code'],
            title=tableinf['title'],
            gov_org=tableinf['gov_org'],
            gov_org_code=tableinf['gov_org_code'],
            survey_date=tableinf['survey_date']
        )
        values, class_object = get_stats_id_value(api_key, stat_id, {})
        for vdata in values:
            if not 'value' in vdata:
                continue
            value_row = StatValue.create(
                stat_id=stat_id,
                value=vdata['value']
            )
            for k, v in vdata.items():
                stat_val_attr = StatValueAttr.create(
                    stat_id=stat_id,
                    stat_value=value_row,
                    attr_name=k,
                    attr_value=v
                )
                if k == 'area':
                    # メッシュコード
                    meshbox = jpgrid.bbox(v)
                    database_proxy.get_conn().execute(
                        'INSERT INTO MapArea(stat_id, stat_val_attr_id, geometry) VALUES(?,?,BuildMBR(?,?,?,?,?))',
                        (stat_id, stat_val_attr.id, meshbox['s'], meshbox['e'], meshbox['n'], meshbox['w'], SRID)
                    )
        database_proxy.commit()

属性名がareaの場合、地域メッシュとして、MapAreaにジオメトリ情報を格納します。
R-Indexを用いている場合、PeeweeのORMで動作させるのが困難なので、ここは直接SQL文を記述しています。

インポート用のスクリプト

以上を利用したインポート用のスクリプトが次の通りです。

https://github.com/mima3/estat/blob/master/import_estat.py

使い方はAPI_KEY、統計のタイトル、mod_spatialite.dllへのパス、DBへのパスを指定して実行します。

python import_estat.py API_KEY 2 平成22年国勢調査-世界測地系(1KMメッシュ)20101001  C:\tool\spatialite\mod_spatialite-4.2.0-win-x86\mod_spatialite.dll estat.sqlite

このスクリプトはWindows用にcp932で文字コードをとりあつかっているので、お使いのターミナルに合わせて修正してください。

WebサーバーでGeoJSONとして地域メッシュの情報を返す

次にWebサーバーで指定の範囲の地域メッシュをGeoJSONとして返す方法について説明します。
ここでは、bottleを用いていますが、この辺りは好きなWebフレームワークで構わないと思いますし、Webフレームワークを使うまでもないかもしれません。

Bottle: Python Web Framework
http://bottlepy.org/docs/dev/index.html

まず、範囲を指定して、DB中から該当の地域メッシュを取得するようなコードを記述します。

estat_db.py

def get_mesh_stat(stat_id_start_str, attr_value, xmin, ymin, xmax, ymax):
    """
    地域メッシュの統計情報を取得する
    @param stat_id_start_str 統計IDの開始文字 この文字から始まるIDをすべて取得する.
    @param attr_value cat01において絞り込む値
    @param xmin 取得範囲
    @param ymin 取得範囲
    @param xmax 取得範囲
    @param ymax 取得範囲
    """
    rows = database_proxy.get_conn().execute("""
        SELECT
          statValue.value,
          AsGeoJson(MapArea.Geometry)
        FROM
          MapArea
          inner join idx_MapArea_Geometry ON pkid = MapArea.id AND xmin > ? AND ymin > ? AND xmax < ? AND ymax < ?
          inner join statValueAttr ON MapArea.stat_val_attr_id = statValueAttr.id
          inner join statValueAttr AS b ON b.stat_value_id = statValueAttr.stat_value_id AND b.attr_value = ?
          inner join statValue ON statValue.id = b.stat_value_id
        WHERE
          MapArea.stat_id like ?;
    """, (xmin, ymin, xmax, ymax, attr_value, stat_id_start_str + '%'))
    ret = []
    for r in rows:
        ret.append({
            'value': r[0],
            'geometory': r[1]
        })
    return ret

GeoJSONとして扱うために、AsGeoJson(MapArea.Geometry)でGeometryをGeoJSONに変換しています。
これを以下のように、結合して、一つのGeoJSONとしてクライアントに返します。

参照元:
https://github.com/mima3/estat/blob/master/application.py

application.py

@app.get('/json/get_population')
def getPopulation():
    stat_id = request.query.stat_id
    swlat = request.query.swlat
    swlng = request.query.swlng
    nelat = request.query.nelat
    nelng = request.query.nelng
    attrval = request.query.attr_value

    ret = estat_db.get_mesh_stat(stat_id, attrval, swlng, swlat, nelng, nelat)
    res = {'type': 'FeatureCollection', 'features': []}
    for r in ret:
        item = {
            'type': 'Feature',
            'geometry': json.loads(r['geometory']),
            'properties': {'value': r['value']}
        }
        res['features'].append(item)
    response.content_type = 'application/json;charset=utf-8'
    return json.dumps(res)

これにより、次のようなリクエストが処理できるようになります。

http://needtec.sakura.ne.jp/estat/json/get_population?stat_id=T000608&attr_value=%E4%BA%BA%E5%8F%A3%E7%B7%8F%E6%95%B0&swlat=35.503426100823496&swlng=139.53192492382811&nelat=35.83811583873688&nelng=140.08124133007811

レスポンス:

{"type": "FeatureCollection", "features": [{"geometry": {"type": "Polygon", "coordinates": [[[139.5374999999999, 35.6], [139.5499999999999, 35.6], [139.5499999999999, 35.60833333333333], [139.5374999999999, 35.60833333333333], [139.5374999999999, 35.6]]]},略 

GoogleMapを用いた地域メッシュの描画の例

ここではGoogleMapを用いた地域メッシュの描画例を記述します。

mesh.png

デモ:
http://needtec.sakura.ne.jp/estat/population

GoogleMapではaddGeoJSONを用いることで、任意のGeoJSONをGoogleMAP上に描画できます。

population.js

          features  = map.data.addGeoJson(result);
          var max = 0;
          for (var i = 0; i < features.length; i++) {
            if (max < features[i].getProperty('value')) {
              max = features[i].getProperty('value');
            }
          }
          map.data.setStyle(styleFeature(max));

この際、setStyleで、GeoJSONのスタイルを指定できます。この例ではプロパティvalueに応じて色の濃さを変更しています。

population.js

      var styleFeature = function(max) {
        var colorScale = d3.scale.linear().domain([0, max]).range(["#CCFFCC", "red"]);
        return function(feature) {
          return {
            strokeWeight : 1,
            fillColor: colorScale(+feature.getProperty('value')),
            fillOpacity: 0.5
          };
        };
      }

参考:
Google Map上にGeoJSONデータを表示する
http://shimz.me/blog/google-map-api/3445

あまり、大きな範囲を一度に描画すると重いので、拡大機能を無効にしています。

D3.jsを用いた地域メッシュの描画の例

ここではD3.jsを用いてSVGによる地域メッシュの描画例を記述します。

メッシュ.png

デモ:
http://needtec.sakura.ne.jp/tokuraku/passenger.html

passenger.js

$('#selMesh').change(function() {
  svgMeshGrp
     .attr('class', 'tracts')
     .selectAll('rect')
     .data([])
     .exit()
     .remove();

  var sel = $('#selMesh').val();
  if (!sel) {
    return;
  }

  function drawMesh(json) {
    console.log(json);
    var max = 0;
    for (var i = 0; i < json.features.length; ++i) {
      var v = parseInt(json.features[i].properties.value);
      if (max < v) {
        max = v;
      }
    }
    console.log(max);
    var colorScale = d3.scale.linear().domain([0, max]).range([0.0, 0.8]);
    svgMeshGrp
       .attr('class', 'tracts')
       .selectAll('rect')
       .data(json.features)
       .enter()
       .append('rect')
       .attr('x', function(d, i) {
         var extX = d3.extent(d.geometry.coordinates[0], function(d) { return d[0];});
         var extY = d3.extent(d.geometry.coordinates[0], function(d) { return d[1];});
         var pt = projection([extX[0], extY[0]]);
         return pt[0];
       })
       .attr('y', function(d) {
         var extX = d3.extent(d.geometry.coordinates[0], function(d) { return d[0];});
         var extY = d3.extent(d.geometry.coordinates[0], function(d) { return d[1];});
         var pt = projection([extX[0], extY[0]]);
         return pt[1];
       })
       .attr('width', function(d) {
         var extX = d3.extent(d.geometry.coordinates[0], function(d) { return d[0];});
         var extY = d3.extent(d.geometry.coordinates[0], function(d) { return d[1];});
         var ptMin = projection([extX[0], extY[0]]);
         var ptMax = projection([extX[1], extY[1]]);
         return Math.abs(ptMax[0] - ptMin[0]);
       })
       .attr('height', function(d) {
         var extX = d3.extent(d.geometry.coordinates[0], function(d) { return d[0];});
         var extY = d3.extent(d.geometry.coordinates[0], function(d) { return d[1];});
         var ptMin = projection([extX[0], extY[0]]);
         var ptMax = projection([extX[1], extY[1]]);
         return Math.abs(ptMax[1] - ptMin[1]);
       })
       .attr('fill-opacity', function(d) {
         console.log('color' , d.properties.value, colorScale(d.properties.value));
         return colorScale(d.properties.value);
       })
       .attr('fill' , '#00f');
  }

  if (stat_dict[sel].json) {
    drawMesh(stat_dict[sel].json);
  } else {
    $.blockUI({ message: '<img src="/railway_location/img/loading.gif" />' });
    var api = encodeURI(util.format('/estat/json/get_population?swlat=%d&swlng=%d&nelat=%d&nelng=%d&stat_id=%s&attr_value=%s',
      swlat, swlng, nelat, nelng, stat_dict[sel].id, stat_dict[sel].attrval
    ));
    d3.json(api, function(json) {
      drawMesh(json);
      stat_dict[sel].json = json;
      $.unblockUI();
    });
  }
}).keyup(function() {
  $(this).blur().focus();
});

ここでの注意点は、GeoJSONを使用して、pathとして描画せずに、rectで描画している点です。
svgにpathを1000個追加すると重くなりますが、rectであればそれよりも軽いです。
地域メッシュは必ず四角形なのでrectで描画した方がいいでしょう。

まとめ

政府総合の統計窓口eStatで地域メッシュを記述する方法を説明しました。
SpatialiteなどによるGeometry情報を格納できるDBに一旦格納することで、処理速度があがることが実証できたかと思います。

PHPで地図とかの空間情報をDBに格納しようとしてSpatialiteを使うと苦労する

目的

この記事ではPHPでSpatialiteを使う方法について解説します。

Spatialiteのビルドなどについては下記を参考にしてください。

地図とかの空間情報をSQLiteに格納するSpatiaLiteを使用してみる
http://qiita.com/mima_ita/items/64f6c2b8bb47c4b5b391

最終目的としては以下のコードが動作することを目標とします。

<?php

class MyDB extends SQLite3
{
    function __construct()
    {
        $this->open(':memory:');
    }
}

$ver = SQLite3::version();
var_dump($ver);

$db = new MyDB();
if (!$db->loadExtension('mod_spatialite.dll')) {
  exit();
}

# reporting some version info
$rs = $db->query('SELECT sqlite_version()');
while ($row = $rs->fetchArray())
{
  print "SQLite version: $row[0]" . PHP_EOL;
}

$rs = $db->query('SELECT spatialite_version()');
while ($row = $rs->fetchArray())
{
  print "SpatiaLite version: $row[0]" . PHP_EOL;
}

$db->exec('SELECT InitSpatialMetaData()');
$db->exec('CREATE TABLE foo (name TEXT)');
$db->exec('SELECT AddGeometryColumn("foo", "Geometry", 0, "POINT", 2)');

# SQLite 3.7.17以降じゃないとRTreeインデックスは作れない
if ($ver->versionNumber < 3007017) {
  print "RTree-Index is Not supported. (< 3.7.17)" . PHP_EOL;
} else {
  $db->exec('SELECT CreateSpatialIndex("foo", "Geometry")');
}
$db->exec('INSERT INTO foo VALUES("test1", GeometryFromText("POINT(140 30)"))');
$db->exec('INSERT INTO foo VALUES("test2", GeometryFromText("POINT(135 33)"))');

$rs = $db->query('SELECT name, AsGeoJson("Geometry") AS geo FROM foo');

while ($row = $rs->fetchArray())
{
  var_dump($row);
}

$db->close();

?>

前提

PHP5.5以降でないとloadExtensionは動作しません。
少なくとも、PHP5.3では基本無理(※)で、PHP5.6ではいけました。
なので、特に制約がなければPHPの最新を用意してください。

なお、基本無理というのは、機能を落として、なんやかんやすれば、なんとかなりはします。

あと、Windowsユーザの人は苦労します。LinuxとかMacのパスの区切りが"/"の環境だと比較的楽です。

php.iniの修正

loadExtensionを用いて拡張モジュールを使う場合、php.iniを修正して、拡張モジュールを読み込めるようにする必要があります。

Windowsの場合はデフォルトでSQLite3が無効になっている場合もあるのでphp.iniのコメントアウトを外してください。

php.ini

extension=php_sqlite3.dll

次に拡張モジュールを格納しているフォルダを指定します。

php.ini

[sqlite3]
sqlite3.extension_dir =C:\tool\spatialite\mod_spatialite-4.2.0-win-x86

PHPのSQLiteでは、ここで指定したフォルダ以外から拡張モジュールをロードすることはできません。

あとは環境変数のPATHに上記のフォルダを指定すれば、ラッキーガイな環境では動作します。お疲れ様でした。

以下のエラーが出る場合は、ここからが本当の地獄です。

PHP Warning:  SQLite3::loadExtension(): 指定されたプロシージャが見つかりません。

何でエラーになるか

おそらく、ここでエラーになるのはPHPが使用しているSQLiteのバージョンが古いか、Windowsユーザの人だと思います。

なぜエラーになるかの前にSQLiteの拡張モジュールのエントリーポイントが関係しています。

Run-Time Loadable Extensions
http://www.sqlite.org/loadext.html

本来のSQLiteのload_extensionはload_extension(X,Y)となっており、2つの引数を取ります。Xがモジュール名でYがエントリーポイントの関数名となります。
このYは省略可能で、PHPのloadExtension関数の場合、省略して実行されています。

この省略された場合の挙動は次のようになります。
・もし省略されている場合は sqlite3_extension_init を使用する
・これでロードできない場合は、最後の"/" 以降を使用する。この時、接頭語のlibは除去し、アルファベットのみを小文字にし、拡張子は考慮しない。
具体的には次のようになります。

"/usr/lib/libmathfunc-4.8.so" ⇒ "sqlite3_mathfunc_init".
"./SpellFixExt.dll" ⇒"sqlite3_spellfixext_init"

なお、mod_spatialite.dllのエントリーポイントは「sqlite3_modspatialite_init」となっています。

これについては、Dependency.Walkerで調べるかソースコード読んでください。

spatialite.png
http://www.dependencywalker.com/

つまり、PHPからloadExtensionを行った場合にエントリーポイントを「sqlite3_modspatialite_init」に変換できていないのが、このエラーの原因となります。

何故、古いPHPではsqlite3_modspatialite_initというエントリーポイントを取得できないか

PHP5.4以前の場合は話は簡単です。古いSQLiteの仕様にはこの機能がなかったものと推測できます。
実際、PHP5.4にバンドルされているSQLiteのコードに、該当の処理は存在しません。

https://github.com/php/php-src/blob/PHP-5.4/ext/sqlite3/libsqlite/sqlite3.c

static int sqlite3LoadExtension(
  sqlite3 *db,          /* Load the extension into this database connection */
  const char *zFile,    /* Name of the shared library containing extension */
  const char *zProc,    /* Entry point.  Use "sqlite3_extension_init" if 0 */
  char **pzErrMsg       /* Put error message here if not 0 */
){
// 略
  if( zProc==0 ){
    zProc = "sqlite3_extension_init";
  }

  handle = sqlite3OsDlOpen(pVfs, zFile);
  if( handle==0 ){
    if( pzErrMsg ){
      *pzErrMsg = zErrmsg = sqlite3_malloc(nMsg);
      if( zErrmsg ){
        sqlite3_snprintf(nMsg, zErrmsg, 
            "unable to open shared library [%s]", zFile);
        sqlite3OsDlError(pVfs, nMsg-1, zErrmsg);
      }
    }
    return SQLITE_ERROR;
  }
  xInit = (int(*)(sqlite3*,char**,const sqlite3_api_routines*))
                   sqlite3OsDlSym(pVfs, handle, zProc);
  if( xInit==0 ){
    if( pzErrMsg ){
      *pzErrMsg = zErrmsg = sqlite3_malloc(nMsg);
      if( zErrmsg ){
        sqlite3_snprintf(nMsg, zErrmsg,
            "no entry point [%s] in shared library [%s]", zProc,zFile);
        sqlite3OsDlError(pVfs, nMsg-1, zErrmsg);
      }
      sqlite3OsDlClose(pVfs, handle);
    }
    return SQLITE_ERROR;
  }else if( xInit(db, &zErrmsg, &sqlite3Apis) ){
    if( pzErrMsg ){
      *pzErrMsg = sqlite3_mprintf("error during initialization: %s", zErrmsg);
    }
    sqlite3_free(zErrmsg);
    sqlite3OsDlClose(pVfs, handle);
    return SQLITE_ERROR;
  }

見てわかる通り、sqlite3OsDlSymでxInitが取得できなかった場合にエラーとして返しています。

一方、PHP5.5では、xInitが取得できなかった場合に、代替のエントリーポイントを取得する処理が記述されています。

https://raw.githubusercontent.com/php/php-src/PHP-5.5/ext/sqlite3/libsqlite/sqlite3.c

static int sqlite3LoadExtension(
  sqlite3 *db,          /* Load the extension into this database connection */
  const char *zFile,    /* Name of the shared library containing extension */
  const char *zProc,    /* Entry point.  Use "sqlite3_extension_init" if 0 */
  char **pzErrMsg       /* Put error message here if not 0 */
){
// 略
  xInit = (int(*)(sqlite3*,char**,const sqlite3_api_routines*))
                   sqlite3OsDlSym(pVfs, handle, zEntry);

  /* If no entry point was specified and the default legacy
  ** entry point name "sqlite3_extension_init" was not found, then
  ** construct an entry point name "sqlite3_X_init" where the X is
  ** replaced by the lowercase value of every ASCII alphabetic 
  ** character in the filename after the last "/" upto the first ".",
  ** and eliding the first three characters if they are "lib".  
  ** Examples:
  **
  **    /usr/local/lib/libExample5.4.3.so ==>  sqlite3_example_init
  **    C:/lib/mathfuncs.dll              ==>  sqlite3_mathfuncs_init
  */
  if( xInit==0 && zProc==0 ){
    int iFile, iEntry, c;
    int ncFile = sqlite3Strlen30(zFile);
    zAltEntry = sqlite3_malloc(ncFile+30);
    if( zAltEntry==0 ){
      sqlite3OsDlClose(pVfs, handle);
      return SQLITE_NOMEM;
    }
    memcpy(zAltEntry, "sqlite3_", 8);
    for(iFile=ncFile-1; iFile>=0 && zFile[iFile]!='/'; iFile--){}
    iFile++;
    if( sqlite3_strnicmp(zFile+iFile, "lib", 3)==0 ) iFile += 3;
    for(iEntry=8; (c = zFile[iFile])!=0 && c!='.'; iFile++){
      if( sqlite3Isalpha(c) ){
        zAltEntry[iEntry++] = (char)sqlite3UpperToLower[(unsigned)c];
      }
    }
    memcpy(zAltEntry+iEntry, "_init", 6);
    zEntry = zAltEntry;
    xInit = (int(*)(sqlite3*,char**,const sqlite3_api_routines*))
                     sqlite3OsDlSym(pVfs, handle, zEntry);
  }

このように古いPHPだとバンドルされているSQLiteのエントリーポイントの取得方法が異なるために、エラーになるのです。

何故、Windowsではsqlite3_modspatialite_initというエントリーポイントを取得できないか

次にWindowsではどうして新しいPHPでもエントリーポイントを取得できないか説明します。
PHPはsqlite3.extension_dir で指定したフォルダ名と、loadExtensionで指定したファイル名を組み合わせて、sqlite3LoadExtensionを実行しています。

つまり、sqlite3LoadExtensionに渡されるzFileは次のようになります。

C:\tool\spatialite\mod_spatialite-4.2.0-win-x86\mod_spatialite.dll

UNIXの場合、パスの区切りが「/」なので適切にファイル名を抽出しますが、Windowsだと「\」なので抽出しません。この場合、エントリーポイントとして期待されるのが以下のようになってしまうのです。

ctoolspatialitemod_spatialitewinmod_spatialite

対応策

ようするにエントリーポイントを認識させればいいのです。
残念なことにPHPでは、以下のようなSQLを認めていません。

SELECT load_extension("mod_spatialite", "sqlite3_modspatialite_init")'

https://github.com/php/php-src/blob/PHP-5.5/ext/sqlite3/sqlite3.c

sqlite3.c

PHP_METHOD(sqlite3, loadExtension)
{
    php_sqlite3_db_object *db_obj;
    zval *object = getThis();
    char *extension, *lib_path, *extension_dir, *errtext = NULL;
    char fullpath[MAXPATHLEN];
    int extension_len, extension_dir_len;
    db_obj = (php_sqlite3_db_object *)zend_object_store_get_object(object TSRMLS_CC);

    SQLITE3_CHECK_INITIALIZED(db_obj, db_obj->initialised, SQLite3)

    if (FAILURE == zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "s", &extension, &extension_len)) {
        return;
    }

# ifdef ZTS
    if ((strncmp(sapi_module.name, "cgi", 3) != 0) &&
        (strcmp(sapi_module.name, "cli") != 0) &&
        (strncmp(sapi_module.name, "embed", 5) != 0)
    ) {     php_sqlite3_error(db_obj, "Not supported in multithreaded Web servers");
        RETURN_FALSE;
    }
# endif

    if (!SQLITE3G(extension_dir)) {
        php_sqlite3_error(db_obj, "SQLite Extension are disabled");
        RETURN_FALSE;
    }

    if (extension_len == 0) {
        php_sqlite3_error(db_obj, "Empty string as an extension");
        RETURN_FALSE;
    }

    extension_dir = SQLITE3G(extension_dir);
    extension_dir_len = strlen(SQLITE3G(extension_dir));

    if (IS_SLASH(extension_dir[extension_dir_len-1])) {
        spprintf(&lib_path, 0, "%s%s", extension_dir, extension);
    } else {
        spprintf(&lib_path, 0, "%s%c%s", extension_dir, DEFAULT_SLASH, extension);
    }

    if (!VCWD_REALPATH(lib_path, fullpath)) {
        php_sqlite3_error(db_obj, "Unable to load extension at '%s'", lib_path);
        efree(lib_path);
        RETURN_FALSE;
    }

    efree(lib_path);

    if (strncmp(fullpath, extension_dir, extension_dir_len) != 0) {
        php_sqlite3_error(db_obj, "Unable to open extensions outside the defined directory");
        RETURN_FALSE;
    }

    /* Extension loading should only be enabled for when we attempt to load */
    sqlite3_enable_load_extension(db_obj->db, 1);
    if (sqlite3_load_extension(db_obj->db, fullpath, 0, &errtext) != SQLITE_OK) {
        php_sqlite3_error(db_obj, "%s", errtext);
        sqlite3_free(errtext);
        sqlite3_enable_load_extension(db_obj->db, 0);
        RETURN_FALSE;
    }
    sqlite3_enable_load_extension(db_obj->db, 0);

    RETURN_TRUE;
}

loadExtensionをするまえにsqlite3_enable_load_extensionで拡張モジュールを許可して、処理が終了したら不許可にしていることがわかります。
つまり、PHPの思想としてはSELECTでload_extensionはさせないというものになります。

これ以外に回避する方法は3つあります。

1つ mod_spatialiteのソースコードでエントリーポイント名を変更する
2つ mod_spatialite.dllをバイナリエディタで開いて、エントリーポイント名を改ざんする
3つ sqliteのコードで「\」も考慮するようにする。

以上の3つです。

mod_spatialiteのソースコードでエントリーポイント名を変更する

これに関しては未検証です。
ただ、エントリーポイントを修正するのは楽でも、mod_spatialiteをWindowsでビルドする環境を作るのはしんどいと思います。

mod_spatialite.dllをバイナリエディタで開いて、エントリーポイント名を改ざんする

spatialite_bin.png

spatialite_bin2.png

sqlite3_modspatialite_initという文字をバイナリエディタで検索して、sqlite3_extension_initに置き換えます。この際、足りないところには00で埋めます。削除したりすると、アドレスが変わるので動作しなくなります。

この方法で対応した場合、RTreeインデックスを使用できなくなることに目をつぶればPHP5.3でもmod_spatialiteを利用できます。

sqliteのコードで「\」も考慮するようにする。

そもそも論として、sqlite3.cを以下のように修正すれば「\」でもファイル名のみを抽出します。

sqlite3.c

  if( xInit==0 && zProc==0 ){
    int iFile, iEntry, c;
    int ncFile = sqlite3Strlen30(zFile);
    zAltEntry = sqlite3_malloc(ncFile+30);
    if( zAltEntry==0 ){
      sqlite3OsDlClose(pVfs, handle);
      return SQLITE_NOMEM;
    }
    memcpy(zAltEntry, "sqlite3_", 8);
    // Windows 対応
    //for(iFile=ncFile-1; iFile>=0 && zFile[iFile]!='/'; iFile--){}
    for(iFile=ncFile-1; iFile>=0 && (zFile[iFile]!='/' && zFile[iFile]!='\\'); iFile--){}

あとはPHPのソースコードをWindowsでビルドする方法ですが、以下を参考にすると良いでしょう。

PHP拡張モジュールを Windows でビルド
http://ngyuki.hatenablog.com/entry/20120625/p1

PHP を Windows でビルド
http://ngyuki.hatenablog.com/entry/20120701/p1

php_sqlite3.dllを作成するには次のようなconfigureでいけました。

configure --disable-all --enable-cli --enable-cgi  --with-sqlite3=shared
nmake

一応x86でPHP5.6のSQLiteをビルドしたものを以下に置きます。
http://needtec.sakura.ne.jp/release/php_sqlite3.zip

まとめ

このように、SpatialiteをPHPで使うにはかなり苦労が必要です。

「公共施設等情報のオープンデータ実証」を使ってみよう

現在、福岡県の福岡市、糸島市において公共施設情報や統計情報を取得できるAPIが公開されています。

公共施設等情報のオープンデータ実証 開発者サイト
http://teapot.bodic.org/

この記事ではPHP、JavaScriptを用いてこのAPIを操作する事例を記述します。

これらを利用して以下のような物が作れます。

福岡市病院マップ
http://needtec.sakura.ne.jp/fukuoka_map/page/hospital_map?lang=ja

福岡市災害マップ
http://needtec.sakura.ne.jp/fukuoka_map/page/disaster_map?lang=ja

Github:
https://github.com/mima3/fukuoka_map/

提供中のAPI

現在提供されているAPIは大きく2種類存在します。
1つはSPARQLクエリーを用いた検索用のAPIで、もう一つは経度緯度より検索する地理的検索APIです。
これらのAPIは開発者サイトの以下のページから実験できます。

http://teapot.bodic.org/test_api.html

RDFについて

Resource Description Framework(RDF)はWeb上にあるリソースを記述するために統一された枠組みです。

RDFはトリプルと言われる主語、述語(プロパティともいう)、目的語の集合です。
主語、述語、目的語のトリプルの1組をRDFグラフとよび、ノードと有向グラフで表現できます。

sparql1.png

実際に「公共施設等情報のオープンデータ」のデータの一部のデータを見てみましょう。
この例では「ひかり保育園〒815_0082」を主語としたデータの一部を表示しています。

sparql2.png

主語と述語はURIで表現しています。
目的語は、「ttp://teapot.bodic.org/dataset/福岡市公共施設」というURIと「福岡市」、「ひかり保育園」、「2014-03-31」といったリテラルで表現しています。

次に複雑な構造をもっているデータの例をみてみましょう。
たとえば、所在地といった場合、県、市、区、町・・・といった様々な要素から構成されています。これらは空白ノードを用いることで、構造体として表現できます。

sparql3.png

この例では「ttp://teapot.bodic.org/predicate/所在地」は空白ノードを目的語としています。
そして、その空白ノードを主語として、県、市、区、町を表現します。

RDFについての詳細は以下の文献を参考にしてください。

RDF入門
http://www.asahi-net.or.jp/~ax2s-kmtn/internet/rdf/rdf-primer.html

RDF(Resource Description Framework):概念および抽象構文
http://www.asahi-net.or.jp/~ax2s-kmtn/internet/rdf/rdf-concepts.html

SPARQLを使用して「公共施設等情報のオープンデータ」を探索する

SPARQLはRDF用のクエリー言語で、RDFに対しての操作を行えます。
下記のURLにその詳細があります。

RDF用クエリ言語SPARQL
http://www.asahi-net.or.jp/~ax2s-kmtn/internet/rdf/rdf-sparql-query.html

この章では、実際に「公共施設等情報のオープンデータ」をSPARQLで探索してみましょう。

「公共施設等情報のオープンデータ」に対しては以下のページからSPARQLを実行できます。

APIを試す(公式)
http://teapot.bodic.org/test_api.html

teapotのSPARQL実行
http://needtec.sakura.ne.jp/fukuoka_map/page/teapot_sparql?lang=ja

今回は、「teapotのSPARQL実行」を用いて説明します。

もっとも単純な例

まずは、なんでもいいので、トリプルを取得するクエリーを実行します。

select * where { ?s ?p ?o . } LIMIT 1

ここでは、構文の詳細は説明しませんが、SQLと似た感じであることがわかります。
この実行結果は以下のようになります。

s p o
http://teapot.bodic.org/facility/福岡県財産9200375 (uri) http://teapot.bodic.org/predicate/用途 (uri) 体育館 (literal)

もしかしたら、データの内容は異なるかもしれませんが、列数と行数は同じになります。

なお、()の中身はデータの型を表します。以下のいずれかになります。

type 説明
uri URI
literal リテラル
bnode 空白ノード

では、where区を次のように変更して実行してみましょう。

select * where { ?a ?b ?c . } LIMIT 1

こんどは、列名が「a,b,c」となることがわかります。

a b c
http://teapot.bodic.org/facility/福岡県財産9200375 (uri) http://teapot.bodic.org/predicate/用途 (uri) 体育館 (literal)

where区には、主語、述語、目的語のトリプルを設定します。
この際、「?」ではじまる文字を指定した場合、それは変数となり、任意の主語、述語、目的語を取得します。

では、特定のuriまたはリテラルでデータを取得する場合、どうしたらいいのでしょうか?
これについては次で説明します。

特定の主語で検索する。

ここでは、主語が「ttp://teapot.bodic.org/facility/ひかり保育園〒815_0082」のデータを取得してみましょう。

select * where { <http://teapot.bodic.org/facility/ひかり保育園〒815_0082> ?p ?o . } 

where区の主語が変数ではなく、明示的にURIを指定します。この際、前後を<>で囲みます。
この結果は次のようになります。

p o
http://teapot.bodic.org/predicate/自治体 (uri) 福岡市 (literal)
http://teapot.bodic.org/predicate/施設名 (uri) ひかり保育園 (literal)
http://teapot.bodic.org/predicate/データセット (uri) http://teapot.bodic.org/dataset/福岡市公共施設 (uri)
http://teapot.bodic.org/predicate/地番(文) (uri) 福岡県福岡市南区大楠三丁目525 (literal)
http://teapot.bodic.org/predicate/種別 (uri) 保育園 (literal)
http://teapot.bodic.org/predicate/担当部局 (uri) 福岡市こども未来局保育課 (literal)
http://teapot.bodic.org/predicate/dataset (uri) http://teapot.bodic.org/dataset/福岡市保育園 (uri)
http://teapot.bodic.org/predicate/addressClean (uri) 福岡県福岡市南区大楠3丁目25番27号 (literal)
http://teapot.bodic.org/predicate/緯度 (uri) 33.5693597 (literal)
http://teapot.bodic.org/predicate/事業者・開設者等 (uri) 福岡市 (literal)
http://teapot.bodic.org/predicate/所在地 (uri) b0 (bnode)
http://www.w3.org/2003/01/geo/wgs84_pos#long (uri) 130.4126095 (literal)
http://teapot.bodic.org/predicate/address (uri) b1 (bnode)
http://www.w3.org/1999/02/22-rdf-syntax-ns#type (uri) http://teapot.bodic.org/type/保育園 (uri)
http://teapot.bodic.org/predicate/dataset (uri) http://teapot.bodic.org/dataset/福岡市公共施設 (uri)
http://teapot.bodic.org/predicate/fax番号 (uri) 521-7120 (literal)
http://www.w3.org/2003/01/geo/wgs84_pos#lat (uri) 33.5693597 (literal)
http://teapot.bodic.org/predicate/所在地(文) (uri) 福岡県福岡市南区大楠3丁目25-27 (literal)
http://teapot.bodic.org/predicate/定員 (uri) 170 (literal)
http://teapot.bodic.org/predicate/データセット (uri) http://teapot.bodic.org/dataset/福岡市保育園 (uri)
http://www.w3.org/2000/01/rdf-schema#label (uri) ひかり保育園 (literal)
http://teapot.bodic.org/predicate/fax (uri) 521-7120 (literal)
http://teapot.bodic.org/predicate/所在地(正) (uri) 福岡県福岡市南区大楠3丁目25番27号 (literal)
http://teapot.bodic.org/predicate/住居表示(文) (uri) 福岡県福岡市南区大楠3丁目25-27 (literal)
http://teapot.bodic.org/predicate/経度 (uri) 130.4126095 (literal)
http://teapot.bodic.org/predicate/郵便番号 (uri) 815-0082 (literal)
http://teapot.bodic.org/predicate/延長保育 (uri) 1時間 (literal)

このように、「ttp://teapot.bodic.org/facility/ひかり保育園〒815_0082」を主語とするデータが取得できました。

特定の述語で検索する。

ここでは、「ttp://teapot.bodic.org/predicate/延長保育」を述語とするデータを取得してみましょう。

select * where { ?s <http://teapot.bodic.org/predicate/延長保育> ?o . } 

今度はwhere区で述語の箇所にURIを明示します。

s o
http://teapot.bodic.org/facility/泰幸保育園〒811_1361 (uri) 1時間 (literal)
http://teapot.bodic.org/facility/こじか保育園〒819_0383 (uri) 1時間 (literal)
http://teapot.bodic.org/facility/ゆなの木保育園〒814_0123 (uri) 1時間 (literal)
http://teapot.bodic.org/facility/信和保育園〒814_0175 (uri) 1時間 (literal)
http://teapot.bodic.org/facility/西新保育園〒814_0004 (uri) 2時間 (literal)

このように、特定の述語(プロパティ)のトリプルの一覧を容易に取得できます。

特定の目的語で検索する

ここでは特定の目的語を検索してみます。
まずは、「ttp://teapot.bodic.org/dataset/福岡市保育園」が目的語のトリプルを取得してみましょう。

select * where { ?s ?p <http://teapot.bodic.org/dataset/福岡市保育園> . } 

この結果は次のようになります。

s p
http://teapot.bodic.org/facility/みとま保育園〒811_0201 (uri) http://teapot.bodic.org/predicate/データセット (uri)
http://teapot.bodic.org/facility/福岡舞鶴誠和保育園〒819_0366 (uri) http://teapot.bodic.org/predicate/データセット (uri)
http://teapot.bodic.org/facility/博多保育園〒813_0041 (uri) http://teapot.bodic.org/predicate/dataset (uri)
http://teapot.bodic.org/facility/さわら保育園〒811_1122 (uri) http://teapot.bodic.org/predicate/データセット (uri)

さて、目的語にはURIだけでなく、リテラルも格納されている可能性があります。
このリテラルを指定して検索する場合は次のようにダブルクォーテーションで値を囲んで実行します。

select * where { ?s ?p "保育園" . }
s p
http://teapot.bodic.org/facility/みどり保育園〒814_0022 (uri) http://teapot.bodic.org/predicate/種別 (uri)
http://teapot.bodic.org/facility/白百合保育園〒814_0104 (uri) http://teapot.bodic.org/predicate/種別 (uri)
http://teapot.bodic.org/facility/光薫寺保育園〒812_0015 (uri) http://teapot.bodic.org/predicate/種別 (uri)
http://teapot.bodic.org/facility/例:内野保育所〒811_1123 (uri) http://teapot.bodic.org/predicate/種別 (uri)

主語、述語、目的語のいずれか2つを指定した検索

ここまでで、主語、述語、目的語のいずれか1つを指定した検索を紹介しました。

主語、述語、目的語は同時に2つ指定可能です。
ここでは、「ttp://teapot.bodic.org/predicate/延長保育」が「3時間」の保育園の一覧を取得してみましょう。

select * where { ?s <http://teapot.bodic.org/predicate/延長保育> "3時間" . } 

この結果は次のようになります。

s
http://teapot.bodic.org/facility/まつぼっくり保育園〒812_0053 (uri)

このように、主語、述語、目的語を指定して検索することが確認できました。

結果を並び替える。

ORDER BY区を使用して結果を並び替えることができます。

select * where { 
  ?s <http://teapot.bodic.org/predicate/延長保育> ?o.
} order by DESC(?o)  

この結果は目的語が大きいもの順に並び替えて出力されます。

s o
http://teapot.bodic.org/facility/玉川保育園分園向野保育園〒815_0035 (uri) - (literal)
http://teapot.bodic.org/facility/荒江保育園分園〒814_0104 (uri) - (literal)
http://teapot.bodic.org/facility/西浦保育園〒819_0202 (uri) - (literal)
http://teapot.bodic.org/facility/観音寺保育園愛宕浜小学校分園〒819_0013 (uri) - (literal)
http://teapot.bodic.org/facility/あすなろ保育園〒811_1323 (uri) 4時間 (literal)
http://teapot.bodic.org/facility/どろんこ保育園〒812_0018 (uri) 4時間 (literal)
http://teapot.bodic.org/facility/舞鶴保育園〒810_0072 (uri) 4時間 (literal)
http://teapot.bodic.org/facility/まつぼっくり保育園〒812_0053 (uri) 3時間 (literal)
http://teapot.bodic.org/facility/あかつき保育園〒813_0036 (uri) 2時間 (literal)
http://teapot.bodic.org/facility/あゆみらい保育園〒819_0041 (uri) 2時間 (literal)

なお、以下のように、DESC()を省略した場合、昇順になります。

select * where { 
  ?s <http://teapot.bodic.org/predicate/延長保育> ?o.
} order by ?o  

複数の列を指定する場合は以下のように、byの後に列を列挙します。

select * where { 
  ?s <http://teapot.bodic.org/predicate/延長保育> ?o.
} order by DESC(?o) ?s

データを分割して取得する

大量の結果が返ってくるクエリーを実行するとエラーが発生します。
たとえば、次のようなクエリーを実行してみましょう。

select * where {
  ?s ?p <http://teapot.bodic.org/dataset/福岡市公共施設> . 
}

この時、以下のようなエラーが発生します。

{"@message":"The SPARQL result is too large. Consider using limit and offset (the server limit is at 2097152 bytes of string)."}

公共施設等情報のオープンデータのSPARQLのAPIでは2MBを超える結果の場合、エラーとして処理されます。
このような大量の結果を取得するには、ORDER BYで結果を並びかえて、OFFSETとLIMITで部分取得します。

select * where {
  ?s ?p <http://teapot.bodic.org/dataset/福岡市公共施設> . 
} ORDER BY ?s ?p OFFSET 0 LIMIT 5

この例では大量の結果のうち、5件のみ抽出しています。

s p
http://teapot.bodic.org/equipement/福岡市財産・NPO・ボランティア交流センター:建物153 (uri) http://teapot.bodic.org/predicate/dataset (uri)
http://teapot.bodic.org/equipement/福岡市財産・NPO・ボランティア交流センター:建物153 (uri) http://teapot.bodic.org/predicate/データセット (uri)
http://teapot.bodic.org/equipement/福岡市財産・あゆみ学園:建物394 (uri) http://teapot.bodic.org/predicate/dataset (uri)
http://teapot.bodic.org/equipement/福岡市財産・あゆみ学園:建物394 (uri) http://teapot.bodic.org/predicate/データセット (uri)
http://teapot.bodic.org/equipement/福岡市財産・かなたけの里公園:建物2953 (uri) http://teapot.bodic.org/predicate/dataset (uri)

次の5件を取得するにはOFFSETを増やします。

select * where {
  ?s ?p <http://teapot.bodic.org/dataset/福岡市公共施設> . 
} ORDER BY ?s ?p OFFSET 5 LIMIT 5
s p
http://teapot.bodic.org/equipement/福岡市財産・かなたけの里公園:建物2953 (uri) http://teapot.bodic.org/predicate/データセット (uri)
http://teapot.bodic.org/equipement/福岡市財産・こうろ館:建物6352 (uri) http://teapot.bodic.org/predicate/dataset (uri)
http://teapot.bodic.org/equipement/福岡市財産・こうろ館:建物6352 (uri) http://teapot.bodic.org/predicate/データセット (uri)
http://teapot.bodic.org/equipement/福岡市財産・こども総合相談センター:建物412 (uri) http://teapot.bodic.org/predicate/dataset (uri)
http://teapot.bodic.org/equipement/福岡市財産・こども総合相談センター:建物412 (uri) http://teapot.bodic.org/predicate/データセット (uri)

このように、OFFSETを増加させていき、データが取得できなくなるまで繰り返せば大量の結果であっても取得することが可能です。

列の絞込みと重複削除

大量のデータをLIMITとOFFSETで分割取得する方法を紹介しましたが、列の絞込みと、同じ行の重複を削除することで結果のサイズを削減することが可能です。

この場合、以下のようにDISTINCT区と、表示する列名を指定してください。

select distinct ?s where {
  ?s ?p <http://teapot.bodic.org/dataset/福岡市公共施設> . 
}

1行に複数の列を指定する

今までは最大で、主語、述語、目的語の3つの列数でした。
しかし、WHERE区を複数書くことで、1行に複数の列を指定することが可能になります。

以下の例では1行に施設名、担当部局、経度、緯度を表示するようにしています。

select distinct ?name ?tantou ?lat ?lng where { 
  ?s <http://teapot.bodic.org/predicate/施設名> ?name . 
  ?s <http://teapot.bodic.org/predicate/担当部局> ?tantou . 
  ?s <http://teapot.bodic.org/predicate/緯度> ?lat. 
  ?s <http://teapot.bodic.org/predicate/緯度> ?lng. 

} limit 100
name tantou lat lng
小呂中学校 (literal) 福岡市教育委員会 (literal) 33.869652 (literal) 33.869652 (literal)
防火女原 (literal) 福岡市消防局管理課 (literal) 33.5732291 (literal) 33.5732291 (literal)
道路百道浜2丁目 (literal) 福岡市道路下水道局道路管理課 (literal) 33.593386 (literal) 33.593386 (literal)
用悪水路立花寺2丁目 (literal) 福岡市道路下水道局下水道河川管理課 (literal) 33.5659894 (literal) 33.5659894 (literal)
道路荒津2丁目 (literal) 福岡市道路下水道局道路管理課 (literal) 33.5916838 (literal) 33.5916838 (literal)
千代1丁目 (literal) 福岡市道路下水道局下水道河川管理課 (literal) 33.6064121 (literal) 33.6064121 (literal)
住宅仮称吉塚西 (literal) 福岡市住宅都市局住宅計画課 (literal) 33.5994292 (literal) 33.5994292 (literal)
見上納骨堂 (literal) 福岡市市民局地域施策課 (literal) 33.5693551 (literal) 33.5693551 (literal)
下山門公衆便所 (literal) 福岡市環境局収集管理課 (literal) 33.580178 (literal) 33.580178 (literal)

データのフィルタリング

結果を減らす別の方法としてはフィルタリングが存在します。
以下の例では特定の「内科」または「アレルギー科」をもつ病院の一覧を取得します。

select distinct ?s where {
  ?s  <http://teapot.bodic.org/predicate/診療科目> ?o. 
   FILTER(?o="アレルギー科" || ?o="内科")
} 
s
http://teapot.bodic.org/facility/南川整形外科病院〒819_8533(医療法人南川整形外科病院) (uri)
http://teapot.bodic.org/facility/高宮外科内科医院〒814_0153(医療法人高宮外科内科医院) (uri)
http://teapot.bodic.org/facility/木村内科医院〒812_0041(医療法人木村内科医院) (uri)
http://teapot.bodic.org/facility/遠藤内科クリニック〒812_0861(医療法人遠藤内科クリニック) (uri)

FILTERには正規表現が使用できます。
以下の例では住所に中央区が含まれるデータを取得します。

select * where {
  ?s <http://teapot.bodic.org/predicate/addressClean> ?o.
  FILTER regex(?o, "中央区", "i").
} 
s o
http://teapot.bodic.org/facility/白金2丁目(福岡市財産5084) (uri) 福岡県福岡市中央区白金2丁目2番21号 (literal)
http://teapot.bodic.org/facility/メンタルクリニック桜坂〒810_0023(渡邊泰良) (uri) 福岡県福岡市中央区警固3丁目6番1号 (literal)
http://teapot.bodic.org/facility/道路平尾4丁目(福岡市財産3916) (uri) 福岡県福岡市中央区平尾4丁目14番31号 (literal)

FILTER区を二つ続けた場合はAND条件となります。
以下の例では診察科目に「内科」または「アレルギー科」をもつ病院で、かつ、住所が「中央区」のものを取得しています。

select distinct * where {
  ?s  <http://teapot.bodic.org/predicate/診療科目> ?kamoku. 
  ?s <http://teapot.bodic.org/predicate/addressClean> ?address.
  FILTER(?kamoku="アレルギー科" || ?kamoku="内科")
  FILTER regex(?address, "中央区", "i").
} 
s kamoku address
http://teapot.bodic.org/facility/みかどクリニック〒810_0041(三角大児) (uri) 内科 (literal) 福岡県福岡市中央区大名2丁目4番33号 (literal)
http://teapot.bodic.org/facility/せと山荘クリニック〒810_0014((医)せと山荘クリニック) (uri) 内科 (literal) 福岡県福岡市中央区平尾5丁目4番21号 (literal)
http://teapot.bodic.org/facility/浜の町病院〒810_8539(国家公務員共済組合連合会) (uri) アレルギー科 (literal) 福岡県福岡市中央区長浜3丁目3番1号 (literal)
http://teapot.bodic.org/facility/天神つじクリニック〒810_0001((医)鸛進会天神つじクリニック) (uri) 内科 (literal) 福岡県福岡市中央区天神3丁目10番11号 (literal)

空白ノードの取り扱い

以下のクエリを実行すると、b0、b1という空白ノードが表示されます。

select * where { <http://teapot.bodic.org/facility/ひかり保育園〒815_0082> ?p ?o . } 
p o
http://teapot.bodic.org/predicate/所在地 (uri) b0 (bnode)
http://teapot.bodic.org/predicate/address (uri) b1 (bnode)

空白ノードの名称は、結果作成時に一意にわりあてられるだけのものです。
たとえば、次のクエリーを実行してみてください。

select * where { <http://teapot.bodic.org/facility/ひかり保育園〒815_0082> ?p ?o . } order by ?o limit 1 

この結果、次の列が表示されます。

p o
http://teapot.bodic.org/predicate/所在地 (uri) b0 (bnode)

では、OFFSET区で続きを取得したら「b1」というノードが取得できるのでしょうか?

select * where { <http://teapot.bodic.org/facility/ひかり保育園〒815_0082> ?p ?o . } order by ?o offset 1 

実際には、以下のように「b0」となります。

p o
http://teapot.bodic.org/predicate/address (uri) b0 (bnode)

これにより、空白ノードの名称は、結果中で一意の名称となるように便宜的に付与されていることがわかります。

では、次に空白ノードに紐づくデータを取得するにはどうしたらいいでしょうか?
この例では、空白ノードの住所に紐ずく県、市などをどうやって取得するかを示します。

select ?p ?o where { 
  { 
    <http://teapot.bodic.org/facility/ひかり保育園〒815_0082> ?p ?o . 
    FILTER(!isBlank(?o))
  }
  union {
    <http://teapot.bodic.org/facility/ひかり保育園〒815_0082> ?x ?y . 
    ?y ?p ?o . 
    FILTER(isBlank(?y))
  } 
}

まず、空白ノード以外のデータを取得します。
次にUNIONにより、そのデータに空白ノードを主語とするデータを結合します。

この実行結果は以下のようになります。

p o
http://teapot.bodic.org/predicate/自治体 (uri) 福岡市 (literal)
http://teapot.bodic.org/predicate/施設名 (uri) ひかり保育園 (literal)
http://teapot.bodic.org/predicate/データセット (uri) http://teapot.bodic.org/dataset/福岡市公共施設 (uri)
http://teapot.bodic.org/predicate/地番(文) (uri) 福岡県福岡市南区大楠三丁目525 (literal)
http://teapot.bodic.org/predicate/種別 (uri) 保育園 (literal)
http://teapot.bodic.org/predicate/担当部局 (uri) 福岡市こども未来局保育課 (literal)
http://teapot.bodic.org/predicate/dataset (uri) http://teapot.bodic.org/dataset/福岡市保育園 (uri)
http://teapot.bodic.org/predicate/addressClean (uri) 福岡県福岡市南区大楠3丁目25番27号 (literal)
http://teapot.bodic.org/predicate/緯度 (uri) 33.5693597 (literal)
http://teapot.bodic.org/predicate/事業者・開設者等 (uri) 福岡市 (literal)
http://www.w3.org/2003/01/geo/wgs84_pos#long (uri) 130.4126095 (literal)
http://www.w3.org/1999/02/22-rdf-syntax-ns#type (uri) http://teapot.bodic.org/type/保育園 (uri)
http://teapot.bodic.org/predicate/dataset (uri) http://teapot.bodic.org/dataset/福岡市公共施設 (uri)
http://teapot.bodic.org/predicate/fax番号 (uri) 521-7120 (literal)
http://www.w3.org/2003/01/geo/wgs84_pos#lat (uri) 33.5693597 (literal)
http://teapot.bodic.org/predicate/所在地(文) (uri) 福岡県福岡市南区大楠3丁目25-27 (literal)
http://teapot.bodic.org/predicate/定員 (uri) 170 (literal)
http://teapot.bodic.org/predicate/データセット (uri) http://teapot.bodic.org/dataset/福岡市保育園 (uri)
http://www.w3.org/2000/01/rdf-schema#label (uri) ひかり保育園 (literal)
http://teapot.bodic.org/predicate/fax (uri) 521-7120 (literal)
http://teapot.bodic.org/predicate/所在地(正) (uri) 福岡県福岡市南区大楠3丁目25番27号 (literal)
http://teapot.bodic.org/predicate/住居表示(文) (uri) 福岡県福岡市南区大楠3丁目25-27 (literal)
http://teapot.bodic.org/predicate/経度 (uri) 130.4126095 (literal)
http://teapot.bodic.org/predicate/郵便番号 (uri) 815-0082 (literal)
http://teapot.bodic.org/predicate/延長保育 (uri) 1時間 (literal)
http://teapot.bodic.org/predicate/県 (uri) 福岡 (literal)
http://teapot.bodic.org/predicate/市 (uri) 福岡 (literal)
http://teapot.bodic.org/predicate/区 (uri) 南 (literal)
http://teapot.bodic.org/predicate/町 (uri) 大楠 (literal)
http://teapot.bodic.org/predicate/丁目 (uri) 3 (literal)
http://teapot.bodic.org/predicate/番 (uri) 25 (literal)
http://teapot.bodic.org/predicate/号 (uri) 27 (literal)
http://www.w3.org/1999/02/22-rdf-syntax-ns#type (uri) http://teapot.bodic.org/type/住居表示 (uri)
http://teapot.bodic.org/predicate/県 (uri) 福岡 (literal)
http://teapot.bodic.org/predicate/市 (uri) 福岡 (literal)
http://teapot.bodic.org/predicate/区 (uri) 南 (literal)
http://teapot.bodic.org/predicate/町 (uri) 大楠 (literal)
http://teapot.bodic.org/predicate/丁目 (uri) 3 (literal)
http://teapot.bodic.org/predicate/番 (uri) 25 (literal)
http://teapot.bodic.org/predicate/号 (uri) 27 (literal)
http://www.w3.org/1999/02/22-rdf-syntax-ns#type (uri) http://teapot.bodic.org/type/住居表示 (uri)

プログラミング言語からAPIを呼ぶ方法

以下のエンドポイントに対してクエリーを指定してgetまたはpostを行うだけです。

http://teapot-api.bodic.org/api/v1/sparql?query=[query]

このAPIについてのドキュメントは下記を参照してください。
http://teapot.bodic.org/api_doc.html

JavaScriptからの実行例

ここではJavaScriptから実行する例を紹介します。

    var sql = "SELECT * WHERE{?s ?p ?o} LIMIT 10"//
    $.post(
      "https://teapot-api.bodic.org/api/v1/sparql",
      {
        query: sql
      },
      function (res) {
        console.log(res);
      },
      "json"
    ).error(function(e){
       console.log("Error: " , e);
    });

実際コードを組む場合、文字列の連結がメンドクサイのでORMっぽく実行できるようにしました。

          var query = new teapot.Query();
          query.columns(['?s', '?p', '?o']);
          query.where('?s', '<http://www.w3.org/1999/02/22-rdf-syntax-ns#type>', '<' + item + '>')
               .where('?s', '?p', '?o')
               .filter('!isBlank(?o)')
               .filter('!isBlank(?s)')
               .union()
                 .where('?s', '<http://www.w3.org/1999/02/22-rdf-syntax-ns#type>', '<' + item + '>')
                 .where('?s', '?p', '?a')
                 .where('?a', '?p', '?o')
                 .filter('isBlank(?a)')
                 .filter('!isBlank(?s)')
               .union()
                 .where('?x', '<http://www.w3.org/1999/02/22-rdf-syntax-ns#type>', '<' + item + '>')
                 .where('?x', '?p', '?o')
                 .where('?s', '?p', '?x')
                 .filter('!isBlank(?o)')
                 .filter('isBlank(?x)')
               .union()
                 .where('?x', '<http://www.w3.org/1999/02/22-rdf-syntax-ns#type>', '<' + item + '>')
                 .where('?x', '?p', '?a')
                 .where('?a', '?p', '?o')
                 .where('?x', '?p', '?a')
                 .where('?s', '?p', '?x')
                 .filter('isBlank(?a)')
                 .filter('!isBlank(?x)')
               .orderby('?s');
          query.executeSpilit(5000, function(err, res) {
            callback(err, res);
          });

このライブラリは以下から取得できます。
https://github.com/mima3/fukuoka_map/blob/master/js/teapot.js

デモ:
http://needtec.sakura.ne.jp/fukuoka_map/page/teapot_sparql?lang=ja
http://needtec.sakura.ne.jp/fukuoka_map/page/teapot_explore?lang=ja

PHPからの実行例

JavaScriptと同様にPHPでもORMっぽく実行できるようにしてみました。

        $query = $query->where('<http://teapot.bodic.org/facility/さくら病院〒814_0142(医療法人社団江頭会さくら病院)>', '?p', '?o')
                       ->filter('!isBlank(?o)')
                       ->orderby('?p ?o');

以下のソースコードをダウンロードしてプロジェクトに取り込むことで使用できます。

https://github.com/mima3/fukuoka_map/tree/master/src/MyLib/TeapotWhere.php
https://github.com/mima3/fukuoka_map/tree/master/src/MyLib/TeapotQuery.php
https://github.com/mima3/fukuoka_map/tree/master/src/MyLib/TeapotCtrl.php
https://github.com/mima3/fukuoka_map/tree/master/src/MyLib/ApiCtrlBase.php

まとめ

このページでは公共施設等情報のオープンデータ実証が提供しているデータの取得方法を説明しました。
また、SPARQLの簡単な使い方と、PHP,JavaScriptからどのように実行するかも説明しました。

これにより、公共施設等情報のオープンデータ実証が提供するデータを自由に探索できるかと思います。

路面管理のビッグデータを利用して路面状況を取得してみる

現在、「路面管理の高度化における実証実験」として、バス,タクシーなどの公共交通車両に取り付けられた安価なセンサーによって取得されたデータを収集し,それらを集約したビッグデータが公開されています。

路面管理ビッグデータのオープンデータコンテスト
http://micrms.force.com/apis

今回はこのデータを利用して地図上で路面の状況を確認できるようにしました。

路面状況Viewer
http://needtec.sakura.ne.jp/road_mgt/road_mgt.html

GitHub
https://github.com/mima3/road_mgt

路面3.png

路面管理の高度化における実証実験 APIについて

APIの仕様書は下記からダウンロードできます。
http://micrms.force.com/developer

データとしては大きく以下の3つに分けることができます。
・路面計測データ i-DRIMSという機械で計測したデータ。加速度、角速度、位置のデータ
・路面状況データ 路面計測データを分析した路面の状況のデータ
・路面緒元データ 自治体毎に短い区間で設定した道路管理情報

通常は路面状況データと路面緒元データを使用することになると思います。

路面状況データの取得方法

路面状況を取得するには以下のURLに対してGETを行います。

http://roadmgt.herokuapp.com/api/v1/datapoints?data_type=road-surface

この結果は次のようなXMLとなります。

<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:geo="http://www.w3.org/2003/01/geo/wgs84_pos#" xmlns:rm="http://roadmgt.herokuapp.com/vocab/rm#" xmlns:rdfs="http://www.w3.org/2000/01/rdf-schema#">
  <rdf:Description rdf:about="http://roadmgt.herokuapp.com/api/v1/datapoints/road-surface_282">
    <rm:iri>2.65</rm:iri>
    <rm:step>0.0</rm:step>
    <rm:patching_num>0</rm:patching_num>
    <rm:municipality_id>1</rm:municipality_id>
    <rm:too_slow_fast_flag>0</rm:too_slow_fast_flag>
    <rm:subsidence_and_puddle>0.0</rm:subsidence_and_puddle>
    <rm:section_id>-1</rm:section_id>
    <rm:speed>37.55711977</rm:speed>
    <rdf:type>http://roadmgt.herokuapp.com/vocab/rm#road-surface</rdf:type>
    <geo:alt>0.0</geo:alt>
    <rm:rms>21.14314346</rm:rms>
    <geo:long>140.1434769</geo:long>
    <rm:rutting_amount>3.75</rm:rutting_amount>
    <rm:pothole_num>0</rm:pothole_num>
    <rm:speed_fluctuation_flag>0</rm:speed_fluctuation_flag>
    <rdfs:label>road-surface_282</rdfs:label>
    <rm:cracking_rate>5.3</rm:cracking_rate>
    <geo:lat>35.61753776</geo:lat>
    <rm:analysis_timestamp>2014-12-30 17:00:02.0</rm:analysis_timestamp>
    <rm:distance>220</rm:distance>
  </rdf:Description>
  <rdf:Description rdf:about="http://roadmgt.herokuapp.com/api/v1/datapoints/road-surface_294">
    <rm:speed_fluctuation_flag>1</rm:speed_fluctuation_flag>
    <rm:patching_num>0</rm:patching_num>
    <rm:iri>2.5</rm:iri>
    <rm:step>0.0</rm:step>
    <rm:too_slow_fast_flag>0</rm:too_slow_fast_flag>
    <rm:subsidence_and_puddle>0.0</rm:subsidence_and_puddle>
    <rm:section_id>-1</rm:section_id>
    <rm:municipality_id>1</rm:municipality_id>
    <rm:distance>340</rm:distance>
    <geo:alt>0.0</geo:alt>
    <rm:cracking_rate>7.5</rm:cracking_rate>
    <geo:long>140.1464737</geo:long>
    <rdf:type>http://roadmgt.herokuapp.com/vocab/rm#road-surface</rdf:type>
    <rm:rms>26.23188158</rm:rms>
    <geo:lat>35.61759593</geo:lat>
    <rm:pothole_num>0</rm:pothole_num>
    <rm:speed>42.15859739</rm:speed>
    <rm:analysis_timestamp>2014-12-30 17:00:02.0</rm:analysis_timestamp>
    <rm:rutting_amount>5.3</rm:rutting_amount>
    <rdfs:label>road-surface_294</rdfs:label>
  </rdf:Description>
  略
</rdf:RDF>

路面状況を全データ取得する方法

APIを実行した時の最大取得件数はデフォルトで300件です。
以降のデータを取得するには、以下のようにoffsetを利用する必要があります。

(1)まず先頭を取得
http://roadmgt.herokuapp.com/api/v1/datapoints?data_type=road-surface&limit=300&offset=0

(2)正常にデータが取得できたら、offsetを変更して実行
http://roadmgt.herokuapp.com/api/v1/datapoints?data_type=road-surface&limit=300&offset=300

(3)これを以下のレスポンスになるまで繰り返す。

404 notFound

路面状況のパラメータ

路面状況はパラメータを指定してデータをフィルタリングすることが可能です。
指定可能なパラメータは次の通りです。

名称 説明
municipality xsd:int 自治体ID
自治体のIDらしいが当然の権利のように説明が一切ないので意味がわかりません。
1と11が入っていてます。おそらく、1が千葉市で11が豊中市だとは思います。
section xsd:int 区間ID。おそらく、道路の区間を内部で一意のIDをふっているのでしょうがよくわかりません。-1とか入っています。
date_start xsd:date 分析日。終了の指定とセットでない場合はエラー(YYYY-MM-DD)
date_end xsd:date 分析日。開始の指定とセットでない場合はエラー(YYYY-MM-DD)
lat_start xsd:double 緯度開始。終了の指定とセットでない場合はエラー
lat_end xsd:double 緯度終了。開始の指定とセットでない場合はエラー
lon_start xsd:double 経度開始。終了の指定とセットでない場合はエラー
lon_end xsd:double 経度終了。開始の指定とセットでない場合はエラー

経度緯度を指定した実行例:

http://roadmgt.herokuapp.com/api/v1/datapoints?data_type=road-surface&lat_start=34.66802946&lon_start=135.362425&lat_end=35.67436512&lon_end=140.2823427&offset=300&limit=300

注意 
ここで紹介したパラメータ並びに、offset,limit以外のパラメータを指定するとエラーとなります。
このことは、以下のようなコードだとエラーになることを意味します。

$.ajax({
  type : 'GET',
  url : 'http://roadmgt.herokuapp.com/api/v1/datapoints',
  cache: false // キャッシュを回避するとエラーになる。
  data : {
    data_type : 'road-surface',
    lat_start  : lat_start,
    lon_start : lon_start,
    lat_end  : lat_end,
    lon_end  : lon_end,
    offset : offset,
    limit : limit
  },
  success : function (data) {
  },
  'error' :  function(xhr, textStatus, error) {
  }
});

このことはIE系でデータが更新されていても、内容が変わらないというリスクを含むことになります。

cahce=falseの場合、パラメータにユニークな一時的な値を指定することで、別のリクエストとして認識させ、キャッシュの使用を回避しているのですが、これが利用できなくなります。

http://stackoverflow.com/questions/4303829/how-to-prevent-a-jquery-ajax-request-from-caching-in-internet-explorer

ターゲットの直接指定

rdf:aboutに記述してある、URIを指定することで直接ターゲットを取得できます。

実行例:
http://roadmgt.herokuapp.com/api/v1/datapoints/road-surface_282

取得結果:

<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:geo="http://www.w3.org/2003/01/geo/wgs84_pos#" xmlns:rm="http://roadmgt.herokuapp.com/vocab/rm#" xmlns:rdfs="http://www.w3.org/2000/01/rdf-schema#">
<rdf:Description rdf:about="http://roadmgt.herokuapp.com/api/v1/datapoints/road-surface_282">
<rm:iri>2.65</rm:iri>
<rm:step>0.0</rm:step>
<rm:patching_num>0</rm:patching_num>
<rm:municipality_id>1</rm:municipality_id>
<rm:too_slow_fast_flag>0</rm:too_slow_fast_flag>
<rm:subsidence_and_puddle>0.0</rm:subsidence_and_puddle>
<rm:section_id>-1</rm:section_id>
<rm:speed>37.55711977</rm:speed>
<rdf:type>http://roadmgt.herokuapp.com/vocab/rm#road-surface</rdf:type>
<geo:alt>0.0</geo:alt>
<rm:rms>21.14314346</rm:rms>
<geo:long>140.1434769</geo:long>
<rm:rutting_amount>3.75</rm:rutting_amount>
<rm:pothole_num>0</rm:pothole_num>
<rm:speed_fluctuation_flag>0</rm:speed_fluctuation_flag>
<rdfs:label>road-surface_282</rdfs:label>
<rm:cracking_rate>5.3</rm:cracking_rate>
<geo:lat>35.61753776</geo:lat>
<rm:analysis_timestamp>2014-12-30 17:00:02.0</rm:analysis_timestamp>
<rm:distance>220</rm:distance>
</rdf:Description>
</rdf:RDF>

ターゲット名の後にプロパティ名を指定することで、プロパティーのみの取得ができます。
この際、「:」を「_」に変更する必要があります。たとえば、「rm:step」を取得したい場合は「rm_step」と指定します。

実行例:
http://roadmgt.herokuapp.com/api/v1/datapoints/road-surface_282/rm_step

取得結果

<rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:geo="http://www.w3.org/2003/01/geo/wgs84_pos#" xmlns:rm="http://roadmgt.herokuapp.com/vocab/rm#">
<rdf:Description rdf:about="http://roadmgt.herokuapp.com/api/v1/datapoints/road-surface_282">
<rm:step>0.0</rm:step>
</rdf:Description>
</rdf:RDF>

レスポンスの説明

主なプロパティを以下に説明します。

プロパティ 名称 説明
rm:analysis_timestamp 分析日時 xsd:dateTime 「2015-01-14 09:01:57」という形式で入っています。
rm:iri IRI xsd:double International Roughness Index 舗装の平坦性(乗り心地)を客観的. に評価する尺度
http://www.kandoken.jp/huroku/130913_2_kouenshiryo.pdf
rm:pothole_num ポットホール数 xsd:int 道路の舗装表面が陥没してできた穴の数
rm:patching_num パッチング数 xsd:int 舗装の損傷に対する応急処置を実施した数
http://www.mlit.go.jp/road/ir/ir-council/pdf/roadstock05.pdf
rm:cracking_rate ひび割れ率 xsd:double http://www.town.oki.fukuoka.jp/file/gyousei/nyuusatsu/FILE_346_45.pdf
rm:rutting_amount わだち掘れ深さ xsd:double http://www.fuzita.org/cod/rut_.html
rm:step 段差 xsd:double
rm:subsidence_and_puddle 沈下・水たまり xsd:double
geo:lat 緯度 xsd:double
geo:long 経度 xsd:double
geo:alt 高度 xsd:double

路面緒元データの取得方法

基本的に路面状況と同じです。
以下のURLに対してGETを行います。

http://roadmgt.herokuapp.com/api/v1/datapoints?data_type=road-master&lat_start=34.66802946&lon_start=135.362425&lat_end=35.67436512&lon_end=140.2823427&offset=300&limit=300

現在は以下の赤線のデータしか存在していません。
路面.png

拡大するとわかりますが、一応、車線毎にデータが存在するようです。

路面2.png

プログラム言語で使用する

Python

ここではPythonで路面状況を取得する方法を記述します。

road_mgt.py

# -*- coding: utf-8 -*-
import urllib
import urllib2
from lxml import etree
import csv
from collections import defaultdict

def get_road_surface(offset, limit):
    """
    路面状況データ
    """
    url = ('http://roadmgt.herokuapp.com/api/v1/datapoints?data_type=road-surface&offset=%d&limit=%d' % (offset, limit))

    req = urllib2.Request(url)
    opener = urllib2.build_opener()
    conn = opener.open(req)
    cont = conn.read()
    parser = etree.XMLParser(recover=True)
    root = etree.fromstring(cont, parser)
    namespaces = {
        'geo' : 'http://www.w3.org/2003/01/geo/wgs84_pos#',
        'rdf' : 'http://www.w3.org/1999/02/22-rdf-syntax-ns#', 
        'rdfs' : 'http://www.w3.org/2000/01/rdf-schema#', 
        'rm' : 'http://roadmgt.herokuapp.com/vocab/rm#'
    }

    row = []
    datas = {}
    value_tags = root.xpath('//rdf:Description', namespaces=namespaces)
    value_tags = root.xpath('//rdf:Description', namespaces=namespaces)
    for value_tag in value_tags:
        label = value_tag.find('rdfs:label', namespaces).text
        step = value_tag.find('rm:step', namespaces).text
        alt = value_tag.find('geo:alt', namespaces).text
        analysis_timestamp = value_tag.find('rm:analysis_timestamp', namespaces).text
        rutting_amount = value_tag.find('rm:rutting_amount', namespaces).text
        municipality_id = value_tag.find('rm:municipality_id', namespaces).text
        speed_fluctuation_flag = value_tag.find('rm:speed_fluctuation_flag', namespaces).text
        section_id = value_tag.find('rm:section_id', namespaces).text
        distance = value_tag.find('rm:distance', namespaces).text
        long = value_tag.find('geo:long', namespaces).text
        iri = value_tag.find('rm:iri', namespaces).text
        cracking_rate = value_tag.find('rm:cracking_rate', namespaces).text
        pothole_num = value_tag.find('rm:pothole_num', namespaces).text
        subsidence_and_puddle = value_tag.find('rm:subsidence_and_puddle', namespaces).text
        speed = value_tag.find('rm:speed', namespaces).text
        rms = value_tag.find('rm:rms', namespaces).text
        lat = value_tag.find('geo:lat', namespaces).text
        too_slow_fast_flag = value_tag.find('rm:too_slow_fast_flag', namespaces).text
        patching_num = value_tag.find('rm:patching_num', namespaces).text

        row.append({
            'label' : label,
            'step' : step,
            'alt' : alt,
            'analysis_timestamp' : analysis_timestamp,
            'rutting_amount' : rutting_amount,
            'municipality_id' : municipality_id,
            'speed_fluctuation_flag' : speed_fluctuation_flag,
            'section_id' : section_id,
            'distance' : distance,
            'long' : long,
            'iri' : iri,
            'cracking_rate' : cracking_rate,
            'pothole_num' : pothole_num,
            'subsidence_and_puddle' : subsidence_and_puddle,
            'speed' : speed,
            'rms' : rms,
            'lat' : lat,
            'too_slow_fast_flag' : too_slow_fast_flag,
            'patching_num' :patching_num
        })
    return row

def get_meas_locale(offset, limit):
    """
    路面計測時の位置データ
    """
    url = ('http://roadmgt.herokuapp.com/api/v1/datapoints?data_type=meas-locale&offset=%d&limit=%d' % (offset, limit))

    req = urllib2.Request(url)
    opener = urllib2.build_opener()
    conn = opener.open(req)
    cont = conn.read()
    parser = etree.XMLParser(recover=True)
    root = etree.fromstring(cont, parser)

    namespaces = {
        'geo' : 'http://www.w3.org/2003/01/geo/wgs84_pos#',
        'rdf' : 'http://www.w3.org/1999/02/22-rdf-syntax-ns#', 
        'rdfs' : 'http://www.w3.org/2000/01/rdf-schema#', 
        'rm' : 'http://roadmgt.herokuapp.com/vocab/rm#'
    }

    row = []
    datas = {}
    value_tags = root.xpath('//rdf:Description', namespaces=namespaces)
    for value_tag in value_tags:
        label = value_tag.find('rdfs:label', namespaces).text
        gpstimestamp = value_tag.find('rm:gpstimestamp', namespaces).text
        course = value_tag.find('rm:course', namespaces).text
        measurement_start_timestamp = value_tag.find('rm:measurement_start_timestamp', namespaces).text
        mounting_direction = value_tag.find('rm:mounting_direction', namespaces).text
        car_model = value_tag.find('rm:car_model', namespaces).text
        long = value_tag.find('geo:long', namespaces).text
        lat = value_tag.find('geo:lat', namespaces).text
        alt = value_tag.find('geo:alt', namespaces).text
        model_number = value_tag.find('rm:model_number', namespaces).text
        car_number = value_tag.find('rm:car_number', namespaces).text
        speed = value_tag.find('rm:speed', namespaces).text
        vertical_accuracy = value_tag.find('rm:vertical_accuracy', namespaces).text
        measurement_timestamp = value_tag.find('rm:measurement_timestamp', namespaces).text
        horizontal_accuracy = value_tag.find('rm:horizontal_accuracy', namespaces).text
        row.append({
            'label' : label,
            'gpstimestamp' : gpstimestamp,
            'course' : course,
            'measurement_start_timestamp' : measurement_start_timestamp,
            'mounting_direction' : mounting_direction,
            'car_model' : car_model,
            'long' : long,
            'lat' : lat,
            'alt' : alt,
            'model_number' : model_number,
            'car_number' : car_number,
            'speed' : speed,
            'vertical_accuracy' : vertical_accuracy,
            'measurement_timestamp' : measurement_timestamp,
            'horizontal_accuracy' : horizontal_accuracy
        })
    return row

def get_road_surface_all():
    ret = []
    limit = 300
    offset = 0
    try:
        while True:
            print ('get_road_surface_all %d' % offset)
            r = get_road_surface(offset, limit)
            ret.extend(r)
            offset += limit
    except urllib2.HTTPError, ex:
        if ex.code == 404:
            return ret
        else:
            raise

get_road_surface_all()を使用すれば、全路面状況が取得できます。
プログラムとしてはurllibでAPIを実行して、レスポンスをlxmlで解析しているだけです。

JavaScript

ここではJavaScriptで路面状況を取得する方法を記述します。

  /**
   * 路面情報のAPIを実行
   */
  function getDataByRange(data_type, lat_start, lon_start, lat_end, lon_end, callback) {
    var offset = 0;
    var limit = 300;
    var obj = {};

    function loadRoadOffset(lat_start, lon_start, lat_end, lon_end, offset, limit, obj, cb) {

      console.log(offset, limit);
      $.ajax({
        type : 'GET',
        url : 'http://roadmgt.herokuapp.com/api/v1/datapoints',
        cache: true , //, cacheをfalseにするとエラーになる。おそらく変なパラメータを無視するのでなく、はじいている
        data : {
          data_type : data_type,
          lat_start  : lat_start,
          lon_start : lon_start,
          lat_end  : lat_end,
          lon_end  : lon_end,
          offset : offset,
          limit : limit
        },
        success : function (data) {
          console.log(data);
          var root = data.documentElement;
          var attrs = root.attributes;
          var records = root.childNodes;
          for(var i=0; i<records.length; i++){
            if(records[i].nodeName.match(/rdf:Description/i)){
              var s = records[i].attributes["rdf:about"].value;
              var props = records[i].childNodes;
              for(var j=0; j<props.length; j++){
                if(props[j].nodeType == 1){
                  var p = props[j].nodeName;
                  var o = props[j].textContent;
                  if (!obj[s]) {
                    obj[s] = {};
                  }
                  if (obj[s][p]) {
                    if (!Array.isArray(obj[s][p])) {
                      var tmp = arys[s][p];
                      obj[s][p] = [];
                      obj[s][p].push(tmp);
                    }
                    obj[s][p].push(o);
                  } else {
                    obj[s][p] = o;
                  }
                }
              }
            }
          }
          loadRoadOffset(lat_start, lon_start, lat_end, lon_end, offset + limit, limit, obj, cb);
        },
        'error' :  function(xhr, textStatus, error) {
          console.log(xhr);
          var err = xhr.responseText;
          if (err == '404 notFound') {
            cb(null , obj);
          } else {
            cb(err , null);
          }
        }
      });
    }
    loadRoadOffset(lat_start, lon_start, lat_end, lon_end, offset, limit, obj, function(err, res) {
      callback(err, res);
    });
  }

この関数は次のように実行します。

getDataByRange('road-master', surfaceRange.lat_min, surfaceRange.long_min, surfaceRange.lat_max, surfaceRange.long_max, function(err, data) {
      console.log('loadRoad', err, data);
});

この処理は指定の範囲の道路緒元データを全て取得します。

まとめ

「路面管理の高度化における実証実験」のAPIを実行することで路面の状況を取得することができます。

これにより道路の劣化状況などの可視化が行えます。

PythonのORMのPeeweeを使ってデータベースを操作してみる(2019年8月版)

概要

peeweeはPython用のORMです。
以下のデータベースの操作が行えます。
・SQLite
・MySQL
・Postgres
・AWSP
・BerkeleyDatabase(未検証)

詳細は下記のドキュメントを参照してください。
https://peewee.readthedocs.org/en/latest/index.html

更新履歴

2019.08.03 サンプルコードを変更 Python2.7 → Python3.7 & Peewee 3.9.6

環境

Windows10
Python 3.7.4(32bit)
Peewee 3.9.6

インストール方法

pip install peewee

簡単な例

SQLiteによる簡単なサンプルを以下に示します。

from peewee import *
from datetime import date

db = SqliteDatabase(':memory:')

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db # This model uses the "people.db" database.

class Pet(Model):
    owner = ForeignKeyField(Person, related_name='pets')
    name = CharField()
    animal_type = CharField()

    class Meta:
        database = db # this model uses the people database

try:
    db.create_tables([Person, Pet])
    with db.transaction():
        # オブジェクトを作ってSaveすることでINSERTする
        uncle_bob = Person(name='Bob', birthday=date(1960, 1, 15), is_relative=True)
        uncle_bob.save() # bob is now stored in the database

        # createでINSERTする
        grandma = Person.create(name='Grandma', birthday=date(1935, 3, 1), is_relative=True)
        herb = Person.create(name='Herb', birthday=date(1950, 5, 5), is_relative=False)

        bob_kitty = Pet.create(owner=uncle_bob, name='Kitty', animal_type='cat')
        herb_fido = Pet.create(owner=herb, name='Fido', animal_type='dog')
        herb_mittens = Pet.create(owner=herb, name='Mittens', animal_type='cat')
        herb_mittens_jr = Pet.create(owner=herb, name='Mittens Jr', animal_type='cat')

        print ("全部取得-----------------")
        for person in Person.select():
            print(person.name, person.is_relative)

        print("catのみ取得-----------------")
        query = Pet.select().where(Pet.animal_type == 'cat')
        for pet in query:
            print(pet.name, pet.owner.name)

        print("Joinの例-----------------")
        query = (Pet
                 .select(Pet, Person)
                 .join(Person)
                 .where(Pet.animal_type == 'cat'))
        for pet in query:
            print(pet.name, pet.owner.name)

        print("更新の例-----------------")
        update_pet = Pet.get(Pet.name=='Kitty')
        update_pet.name = 'Kitty(updated)'
        update_pet.save() 

        query = (Pet
                 .select(Pet, Person)
                 .join(Person)
                 .where(Pet.animal_type == 'cat'))
        for pet in query:
            print(pet.name, pet.owner.name)

        print("削除の例-----------------")
        del_pet = Pet.get(Pet.name=='Mittens Jr')
        del_pet.delete_instance() 

        query = (Pet
                 .select(Pet, Person)
                 .join(Person)
                 .where(Pet.animal_type == 'cat'))
        for pet in query:
            print(pet.name, pet.owner.name)

        db.commit()

except IntegrityError as ex:
    print (ex)
    db.rollback()

この例のようにpeeweeを用いればSQL文を記述することなく、データベースの操作が行えます。
その他、グループ化や、max,min関数を使用したクエリーについては下記を参考にしてください。

https://peewee.readthedocs.org/en/latest/peewee/querying.html

SUBSTRなどの関数を使う場合

SUBSTRなどの関数を使う場合は、SQL()を使用します。

        print ('SUBSTRの例')
        for r in Pet.select(SQL('SUBSTR(name,1,1)').alias('a1')):
            print (r.a1)

この例ではname列の1文字目をSUBSTR関数で取得して、それにa1という列名を与えています。

MAXを使って最大値を取得する場合

fn.MAXで関数を指定して,クエリーからscalar()を使用して値を取得する。

max = LiveChatMessage.select(
    fn.MAX(LiveChatMessage.offset_time_msec)
).where(LiveChatMessage.video_id  == video_id).scalar()
print(max)

大量データの作成

一度に大量のデータを作成する場合は、insert_manyを使用します。

from peewee import *
from datetime import date

db = SqliteDatabase(':memory:')

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db

data_source = [
    {'name' : 'test1' , 'birthday' : date(1960, 1, 15), 'is_relative' : True},
    {'name' : 'test2' , 'birthday' : date(1960, 2, 15), 'is_relative' : True},
    {'name' : 'test3' , 'birthday' : date(1960, 3, 15), 'is_relative' : True},
    {'name' : 'test4' , 'birthday' : date(1960, 4, 15), 'is_relative' : True},
    {'name' : 'test5' , 'birthday' : date(1960, 5, 15), 'is_relative' : True},
]

db.create_tables([Person])
try:
    with db.transaction():
        Person.insert_many(data_source).execute()

    print ("全部取得-----------------")
    for person in Person.select():
        print (person.name, person.birthday, person.is_relative)

except IntegrityError as ex:
    print (ex)
    db.rollback()

DISTINCTの使用方法

DISTINCTは下記のように利用します。

    q = Person.select(Person.name).distinct()
    print (q.sql())
    for person in q:
        print (person.name, person.birthday, person.is_relative)

作成されるSQLは下記のようになります。

'('SELECT DISTINCT "t1"."name" FROM "person" AS "t1"', [])

条件を動的に組み立てる

下記の例では、operation_companyとrailway_lineを条件にクエリを取得しています。

def get_dynamic_sql(name = None, is_relative = None):
    ret = []
    query = Person.select()
    cond = None
    if not name is None:
        cond = (Person.name == name)
    if not is_relative is None:
        if cond:
            cond = cond & (Person.is_relative == is_relative)
        else:
            cond = (Person.is_relative == is_relative)
    rows = query.where(cond)
    for r in rows: # ここでSQLを発行する
        ret.append(r.name)
    return ret

この例ではパラメータの指定方法により4種類のSQLが作成されます。

name is_relative 作られるSQL
None None SELECT "t1"."id", "t1"."name", "t1"."birthday", "t1"."is_relative" FROM "person" AS "t1
None以外 None SELECT "t1"."id", "t1"."name", "t1"."birthday", "t1"."is_relative" FROM "person" AS "t1" WHERE ("t1"."name" = ?)
None None以外 SELECT "t1"."id", "t1"."name", "t1"."birthday", "t1"."is_relative" FROM "person" AS "t1" WHERE ("t1"."is_relative" = ?)
None以外 None以外 SELECT "t1"."id", "t1"."name", "t1"."birthday", "t1"."is_relative" FROM "person" AS "t1" WHERE (("t1"."name" = ?) AND ("t1"."is_relative" = ?))

また、作成されたSQLが実際に発行されるタイミングは使用する時です。
この挙動を利用することにより、動的に複雑な条件のクエリを組み立てることができます。

JOINについて

使用できるJOIN

peewee 2.4.5ではRIGHTやFULL JOINはできません。
INNER JOINかLEFT OUTER JOINのみ使用できます。

query = Curve.select(Curve, RailRoadSection).join(RailRoadSection, JOIN_FULL)

JOIN_FULLを指定した場合、例外が発生します。

peewee.OperationalError: RIGHT and FULL OUTER JOINs are not currently supported

LEFT OUTER JOINの例

取得したレコードに結合中のテーブル名があります。dir関数でレコードの内容をするといいでしょう。

# モデル
from peewee import *
from datetime import date

db = SqliteDatabase(':memory:')

class Group(Model):
    name = CharField()

    class Meta:
        database = db # This model uses the "people.db" database.

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()
    group = ForeignKeyField(Group, related_name='group')

    class Meta:
        database = db # This model uses the "people.db" database.

class Pet(Model):
    owner = ForeignKeyField(Person, related_name='pets', null=True)
    name = CharField()
    animal_type = CharField()

    class Meta:
        database = db # this model uses the people database

try:
    db.create_tables([Person, Pet, Group])
    with db.transaction():
        # 
        grp1 = Group(name='花組')
        grp1.save()
        grp2 = Group(name='奇面組')
        grp2.save()

        # オブジェクトを作ってSaveすることでINSERTする
        uncle_bob = Person(name='Bob', birthday=date(1960, 1, 15), is_relative=True, group=grp1)
        uncle_bob.save() # bob is now stored in the database

        # createでINSERTする
        grandma = Person.create(name='Grandma', birthday=date(1935, 3, 1), is_relative=True, group=grp2)
        herb = Person.create(name='Herb', birthday=date(1950, 5, 5), is_relative=False, group=grp1)

        bob_kitty = Pet.create(owner=herb, name='Kitty', animal_type='cat')
        herb_fido = Pet.create(owner=herb, name='Fido', animal_type='dog')
        herb_mittens = Pet.create(owner=herb, name='Mittens', animal_type='cat')
        herb_mittens_jr = Pet.create(owner=herb, name='Mittens Jr', animal_type='cat')
        ginga = Pet.create(owner=None, name='Mittens Jr', animal_type='cat')

        print ("全部取得-----------------")
        for group in Group.select():
            print(group.name)
        for person in Person.select():
            print(person.name, person.is_relative, person.group)
        for pet in Pet.select():
            print(pet.owner, pet.name, pet.animal_type)

        print("inner Joinの例-----------------")
        query = (Pet
                 .select(Pet, Person)
                 .join(Person))
        for pet in query:
            print(pet.name, pet.owner.name)

        print("left outer Joinの例-----------------")
        query = (Pet
                 .select(Pet, Person)
                 .join(Person, JOIN.LEFT_OUTER))
        for pet in query:
            print(pet.name, pet.owner)

        #print("right outer Joinの例-----------------")
        #query = (Pet
        #         .select(Pet, Person)
        #         .join(Person, JOIN.FULL))
        #for pet in query:
        #    print(pet.name, pet.owner)

except IntegrityError as ex:
    print (ex)
    db.rollback()

複数のテーブルをJOINする場合

複数のテーブルをJOINする場合は、switchを使用してどのテーブルに結合するかを明示しなければなりません。
http://stackoverflow.com/questions/22016778/python-peewee-joins-multiple-tables

# http://stackoverflow.com/questions/22016778/python-peewee-joins-multiple-tables
query = (TimeTableItem
    .select(TimeTableItem, TimeTable, BusStop)
    .join(TimeTable, on = (TimeTableItem.timeTable << list(timetableids.keys())))
    .switch(TimeTableItem)
    .join(BusStop, on=(TimeTableItem.busStop == BusStop.id))
)
for r in query:
    print (r.busStop.stopName)

自己結合

自己結合を行う場合、aliasで別名のオブジェクトを作成しておき、それを利用する

    fromBusStop = BusStopOrder.alias()
    toBusStop = BusStopOrder.alias()
    query = (fromBusStop
        .select(fromBusStop, toBusStop, BusStop)
        .join(
            toBusStop,
            on=((toBusStop.route == fromBusStop.route) & (toBusStop.stopOrder > fromBusStop.stopOrder))
            .alias('toBusStopOrder')
        )
        .switch(toBusStop)
        .join(BusStop, on=(toBusStop.busStop==BusStop.id))
        .where((fromBusStop.busStop==from_bus_stop))
    )
    for r in query:
        print (r.toBusStopOrder.busStop.id)

モデルの作成

ここではモデルの作成について説明します。
詳細は下記を参照してください。
https://peewee.readthedocs.org/en/latest/peewee/models.html

列の型

先の例のようにDateField、CharFieldをModelクラスに指定することでフィールドを設定できます。
ここで使用できるフィールどは下記の通りになります。

Field Type Sqlite Postgresql MySQL
CharField varchar varchar varchar
TextField text text longtext
DateTimeField datetime timestamp datetime
IntegerField integer integer integer
BooleanField smallint boolean bool
FloatField real real real
DoubleField real double precision double precision
BigIntegerField integer bigint bigint
DecimalField decimal numeric numeric
PrimaryKeyField integer serial integer
ForeignKeyField integer integer integer
DateField date date date
TimeField time time time
BlobField blob bytea blob
UUIDField not supported uuid not supported

field作成時のパラメータにより、デフォルト値の指定や重複の有無を指定できます。
詳細は下記を参照してください。

https://peewee.readthedocs.org/en/latest/peewee/api.html#fields

主キーの指定

peeweeで主キーを指定する方法について以下で説明します。

主キーを設定しない場合

PrimaryKeyを明示していない場合は、自動インクリメントのidという主キーを作成します。

class Test1(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db

作成されるSQLiteのデータベース

CREATE TABLE IF NOT EXISTS "test1" ("id" INTEGER NOT NULL PRIMARY KEY, "name" VARCHAR(255) NOT NULL, "birthday" DATE NOT NULL, "is_relative" INTEGER NOT NULL

特定のフィールドを主キーに指定する場合

フィールドを作成する際にprimary_key=Trueと指定することで、該当のフィールドを主キーとします。

class Test2(Model):
    name = CharField(primary_key=True)
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db

作成されるSQLiteのデータベース

CREATE TABLE IF NOT EXISTS "test2" ("name" VARCHAR(255) NOT NULL PRIMARY KEY, "birthday" DATE NOT NULL, "is_relative" INTEGER NOT NULL)

複数のフィールドを主キーにする

CompositeKeyを用いることで複数のフィールドを主キーにできます。

class Test3(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db
        primary_key = CompositeKey('name', 'birthday')

作成されるSQLiteのデータベース

CREATE TABLE IF NOT EXISTS "test3" ("name" VARCHAR(255) NOT NULL, "birthday" DATE NOT NULL, "is_relative" INTEGER NOT NULL, PRIMARY KEY ("name", "birthday"))

インデックスの指定

peeweeでインデックスを指定する方法について以下で説明します。

指定の単一フィールドをインデックスにする

フィールドを作成する際、index=Trueとすることでインデックスにできます

class Test4(Model):
    name = CharField(index=True)
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db

作成されるインデックス

CREATE TABLE IF NOT EXISTS "test4" ("id" INTEGER NOT NULL PRIMARY KEY, "name" VARCHAR(255) NOT NULL, "birthday" DATE NOT NULL, "is_relative" INTEGER NOT NULL)
CREATE INDEX IF NOT EXISTS "test4_name" ON "test4" ("name")

指定のフィールドを組み合わせてインデックスにする

Meta クラスにてindexesを指定することで、複数のキーを組み合わせてインデックスを作成できます

class Test5(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db
        indexes = (
            # 末尾に,がないとエラーになる
            # 複数指定も可能
            (('name', 'birthday'), False),
        )

作成されるインデックス

CREATE TABLE IF NOT EXISTS "test5" ("id" INTEGER NOT NULL PRIMARY KEY, "name" VARCHAR(255) NOT NULL, "birthday" DATE NOT NULL, "is_relative" INTEGER NOT NULL)

CREATE INDEX IF NOT EXISTS "test5_name_birthday" ON "test5" ("name", "birthday")

重複を禁止するインデックスにする

フィールドを作成する際にunique=Trueを指定するか、Metaクラスのindexesの第二引数にTrueを指定することで重複を禁止するインデックスを作成できます。

class Test6(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db
        indexes = (
            # 末尾に,がないとエラーになる
            # 複数指定も可能
            (('name', 'birthday'), True),
        )

作成されるインデックス

CREATE TABLE IF NOT EXISTS "test6" ("id" INTEGER NOT NULL PRIMARY KEY, "name" VARCHAR(255) NOT NULL, "birthday" DATE NOT NULL, "is_relative" INTEGER NOT NULL)
CREATE UNIQUE INDEX IF NOT EXISTS "test6_name_birthday" ON "test6" ("name", "birthday")

外部キーについて

ForeignKeyField()を使用すると外部キーが指定できます。

to_field で主キー以外を指定できますが、このキーは主キーのいずれかであるか、一意制約を持つ必要があります。

http://peewee.readthedocs.org/en/latest/peewee/api.html#ForeignKeyField

様々なデータベースへの接続方法

peeweeでは様々なデータベースを使用できます。
詳細は下記を参照してください。
https://peewee.readthedocs.org/en/latest/peewee/database.html

SQLiteの接続例

memoryまたはファイルを指定して接続ができます。

from peewee import *
from datetime import date

db = SqliteDatabase(':memory:')

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db

class Pet(Model):
    owner = ForeignKeyField(Person, related_name='pets')
    name = CharField()
    animal_type = CharField()

    class Meta:
        database = db

db.create_tables([Person, Pet], True)

APSWの接続例

APSWについての詳細とインストール方法は下記を参照してください。

SQLiteが本気を出せるPythonのライブラリのAPSWを使用してみる
https://needtec.sakura.ne.jp/wod07672/?p=9243

下記でインストールできます。

pip install --user https://github.com/rogerbinns/apsw/releases/download/3.28.0-r1/apsw-3.28.0-r1.zip --global-option=fetch --global-option=--version --global-option=3.28.0 --global-option=--all --global-option=build --global-option=--enable-all-extensions

memoryまたはファイルを指定して接続ができます。

また、withを抜けるときにCommitが実行されるので最後のコミットはコメントアウトしてます(コメントアウトしないとエラーになる)

from peewee import *
from datetime import date
from playhouse.apsw_ext import APSWDatabase
from playhouse.apsw_ext import DateField

db = APSWDatabase('apswdatabase.sqlite')

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db 

class Pet(Model):
    owner = ForeignKeyField(Person, related_name='pets')
    name = CharField()
    animal_type = CharField()

    class Meta:
        database = db 

db.create_tables([Person, Pet])

# db.set_autocommit(False)
with db.transaction():
    # birthday=date(1960, 1, 15) ...
    uncle_bob = Person(name='Bob', birthday=date(1960, 1, 15), is_relative=True)
    uncle_bob.save() # bob is now stored in the database

    grandma = Person.create(name='Grandma', birthday=date(1960, 1, 5), is_relative=True)
    herb = Person.create(name='Herb', birthday='1950-05-05', is_relative=False)

    bob_kitty = Pet.create(owner=uncle_bob, name='Kitty', animal_type='cat')
    herb_fido = Pet.create(owner=herb, name='Fido', animal_type='dog')
    herb_mittens = Pet.create(owner=herb, name='Mittens', animal_type='cat')
    herb_mittens_jr = Pet.create(owner=herb, name='Mittens Jr', animal_type='cat')

    herb_mittens.delete_instance() 
    print ("-----------------")
    for person in Person.select():
        print (person.name, person.is_relative)

    print ("-----------------")
    query = Pet.select().where(Pet.animal_type == 'cat')
    for pet in query:
        print (pet.name, pet.owner.name)

    print ("-----------------")
    query = (Pet
             .select(Pet, Person)
             .join(Person)
             .where(Pet.animal_type == 'cat'))
    for pet in query:
        print (pet.name, pet.owner.name)
    #db.commit() # Not required for APSWDatabase 

MySQLの接続例

データベース、ユーザー、パスワード、ホスト、ポートを指定して接続できます。

from peewee import *
from datetime import date

db = MySQLDatabase(
    database='testdb',
    user='test',
    password="test",
    host="192.168.80.131",
    port=3306)

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db # This model uses the "people.db" database.

class Pet(Model):
    owner = ForeignKeyField(Person, related_name='pets')
    name = CharField()
    animal_type = CharField()

    class Meta:
        database = db # this model uses the people database

db.create_tables([Person, Pet])

接続エラーが出る場合

下記のエラーが発生する場合はmysqlのドライバーをインストールしてください。

peewee.ImproperlyConfigured: MySQL driver not installed!

インストールの例:

pip3 install mysqlclient

https://github.com/coleifer/peewee/issues/1569

もしwindowsを使用していてエラーが出る場合は下記を参照。
https://stackoverflow.com/questions/51294268/pip-install-mysqlclient-returns-fatal-error-c1083-cannot-open-file-mysql-h

環境に応じたwhlファイルをダウンロードしてpip installでそのファイルを指定する。

Postgresの接続例

データベース、ユーザー、パスワード、ホスト、ポートを指定して接続できます。

from peewee import *
from datetime import date
from playhouse.postgres_ext import PostgresqlExtDatabase

# peeweeではデフォルトでhstoreという機能を使うようになっている。
db = PostgresqlExtDatabase(
    database='peewee_test',
    user='postgres',
    password="",
    host="192.168.80.131",
    port=5432,
    register_hstore=False)

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db # This model uses the "people.db" database.

class Pet(Model):
    owner = ForeignKeyField(Person, related_name='pets')
    name = CharField()
    animal_type = CharField()

    class Meta:
        database = db # this model uses the people database

db.create_tables([Person, Pet])

接続エラーの場合

 ERROR: Could not find a version that satisfies the requirement psqlclient (from versions: none)
ERROR: No matching distribution found for psqlclient

psycopg2をインストールする

pip install psycopg2

SpatiaLiteSQLへの接続方法

空間情報を扱うSQLiteの拡張であるSpatiaLiteへの接続について説明します。
http://qiita.com/mima_ita/items/64f6c2b8bb47c4b5b391

これにはplayhouse.sqlite_extのSqliteExtDatabaseを使用して次のように行います。

import os
from peewee import *
from playhouse.sqlite_ext import SqliteExtDatabase

# mod_spatialiteのあるフォルダをPATHに加える
os.environ["PATH"] = os.environ["PATH"] + ';C:\\tool\\spatialite\\mod_spatialite-NG-win-x86'
db = SqliteExtDatabase('testspatialite.sqlite')

class PolygonField(Field):
    db_field = 'polygon'
db.field_overrides = {'polygon': 'POLYGON'}

# mod_spatialiteの読み込み
db.load_extension('mod_spatialite.dll')

class GeometryTable(Model):
  pk_uid  = PrimaryKeyField()
  n03_001 = CharField()
  n03_002 = CharField()
  n03_003 = CharField()
  n03_004 = CharField()
  n03_007 = CharField()
  Geometry = PolygonField()

  class Meta:
      database = db

for r in GeometryTable.select(GeometryTable.n03_001 ,  SQL('AsText(Geometry)').alias('Geo')).limit(10):
    print (r.n03_001, r.Geo)

ポイントは以下の通りです。
・load_extensionを使用してmod_spatialite.dll/soを呼び出す。
・POINTやPOLYGONといったspatialiteの列はFieldクラスを継承して定義しておき、db.field_overrides でコードで指定したdb_fieldとDBの型名を対応付ける
・AsTextなどのspatialite固有の関数はR()を用いて利用する

なお、RTreeIndexのテーブルにアクセスすると下記のようなエラーが出ます。

TypeError: 'idx_Station_geometry' object does not support indexing

この場合は、直接SQL実行してください。(APSWDatabase使用してもダメだった)

    rows = database_proxy.connection().execute("""
        SELECT 
          statValue.value,
          AsGeoJson(MapArea.Geometry)
        FROM 
          MapArea 
          inner join idx_MapArea_Geometry ON pkid = MapArea.id AND xmin > ? AND ymin > ? AND xmax < ? AND ymax < ?
          inner join statValueAttr ON MapArea.stat_val_attr_id = statValueAttr.id 
          inner join statValueAttr AS b ON b.stat_value_id = statValueAttr.stat_value_id AND b.attr_value = ?
          inner join statValue ON statValue.id = b.stat_value_id
        WHERE 
          MapArea.stat_id like ?;
    """,(xmin, ymin, xmax, ymax, attr_value, stat_id_start_str + '%'))

実行時に接続先を選択する方法

いままでの方法では、実行時に接続先を指定することができません。
そこで、Proxy() を使用することで、後から、設定ファイルの情報に合わせた接続を作成することもできます。

from peewee import *
from datetime import date

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

class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = database_proxy

class Pet(Model):
    owner = ForeignKeyField(Person, related_name='pets')
    name = CharField()
    animal_type = CharField()

    class Meta:
        database = database_proxy

# あとから実際のデータベースを指定できる
db = SqliteDatabase(':memory:', autocommit=False)
database_proxy.initialize(db)

db.create_tables([Person, Pet])

Read Slaves やコネクションプール

Databaseによっては、Read Slavesを指定して、読み取り専用のDBにアクセスするようにしたり、コネクションプールを用いて複数の接続を取り扱うことも可能です。

https://peewee.readthedocs.org/en/latest/peewee/database.html#read-slaves
https://peewee.readthedocs.org/en/latest/peewee/database.html#connection-pooling

スキーマーのマイグレーション

Peeweeではスキーマの移行がサポートされています。
他のスキーマー移行ツールとことなり、バージョン管理は行っていませんが、移行を行うためのヘルパー関数が提供されています。

from peewee import *
from datetime import date
from playhouse.migrate import *

db = SqliteDatabase('mig.sqlite')

## 元のDBの作成
class Person(Model):
    name = CharField()
    birthday = DateField()
    is_relative = BooleanField()

    class Meta:
        database = db

db.create_tables([Person])
data_source = [
    {'name' : 'test1' , 'birthday' : date(1960, 1, 15), 'is_relative' : True},
    {'name' : 'test2' , 'birthday' : date(1960, 2, 15), 'is_relative' : True},
    {'name' : 'test3' , 'birthday' : date(1960, 3, 15), 'is_relative' : True},
    {'name' : 'test4' , 'birthday' : date(1960, 4, 15), 'is_relative' : True},
    {'name' : 'test5' , 'birthday' : date(1960, 5, 15), 'is_relative' : True},
    {'name' : 'test1' , 'birthday' : date(1960, 1, 15), 'is_relative' : True},
]
Person.insert_many(data_source).execute()

## スキーマーの移行
migrator = SqliteMigrator(db)

title_field = CharField(default='')
status_field = IntegerField(null=True)

with db.transaction():
    migrate(
        migrator.add_column('Person', 'title', title_field),
        migrator.add_column('Person', 'status', status_field),
        migrator.drop_column('Person', 'is_relative'),
    )

詳細は下記を参考にしてください。
https://peewee.readthedocs.org/en/latest/peewee/playhouse.html#migrate

既存のデータベースからオブジェクトを作成する方法

pwizを用いて既存のデータベースからオブジェクトを作成できます。

python -m pwiz --engine=sqlite mig.sqlite

このコマンドを実行することでpeople.dbからオブジェクトを標準出力します。

出力例

from peewee import *

database = SqliteDatabase('mig.sqlite')

class UnknownField(object):
    def __init__(self, *_, **__): pass

class BaseModel(Model):
    class Meta:
        database = database

class Person(BaseModel):
    birthday = DateField()
    name = CharField()
    status = IntegerField(null=True)
    title = CharField()

    class Meta:
        table_name = 'person'

SQLのロギング

peeweeが発行したSQLをロギングするには以下の通りにします。

import logging
logger = logging.getLogger('peewee')
logger.setLevel(logging.DEBUG)
logger.addHandler(logging.StreamHandler())

これによりpeeweeでSQLを発行するたびに次のようなログがstderrに出力されます。

('SELECT name FROM sqlite_master WHERE type = ? ORDER BY name;', ('table',))
('CREATE TABLE "test1" ("id" INTEGER NOT NULL PRIMARY KEY, "name" VARCHAR(255) NOT NULL, "birthday" DATE NOT NULL, "is_relative" SMALLINT NOT NULL)', [])

まとめ

以上のように、Peeweeを用いることでPythonから様々なDBの操作をSQLを記述せずに行えることが確認できました。

国土数値情報の鉄道データを活用してみる

目的

国土数値情報 ダウンロードサービスでは、国土交通省が管理するデータが取得できます。
今回は、鉄道データを用いて、その座標をGoogleMapにプロットしてみます。

デモ:
http://needtec.sakura.ne.jp/railway_location/railway

GIT:
https://github.com/mima3/railway_location

データについて

鉄道データは下記のページからダウンロードできます。

国土数値情報 鉄道データ
http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N02-v2_2.html

ダウンロードしたファイル中のXMLの使用については下記を参考にしてください。
http://nlftp.mlit.go.jp/ksj/gml/product_spec/KS-PS-N02-v2_1.pdf

簡単に説明すると、鉄道データには、路線の形を示す情報と、駅情報が格納されています。
ここで、重要な要素は以下の通りです。
・gml:Curve 曲線の情報
・ksj:railroadSection 鉄道区間の情報
・ksj:Station 駅情報

Curveには座標情報が格納されています。
railroadSectionとStationのlocation要素にて、Curveへのリンクが格納されています。

データベースへの格納

大量のデータをXMLで取り扱うのは困難なので、一旦リレーショナルデータベースに格納しています。

この時、大きなサイズのXMLを解析しますが、XMLファイル全体を一旦、メモリに格納してパースするというやりかたをすると、メモリ使用量が激増し、処理しきれなくなります。
そのため、lxml.etree.iterparseを用いて、順次処理するようにします。

しかしながら、N02-XX.xmlをlxml.etree.itreparseにて解析するとエラーが発生します。
これは、XML中に以下のような行があるからです。

 xmlns:schemaLocation="http://nlftp.mlit.go.jp/ksj/schemas/ksj-app KsjAppSchema-N02-v2_0.xsd">

lxmlは、ここで指定しているURIを不正なURIとみなし、エラーを出力します。
これを回避するにはlxmlにて、XMLを解析する際に、recover=Trueを指定する必要があります。
http://stackoverflow.com/questions/18692965/how-do-i-skip-validating-the-uri-in-lxml

回避例:

        context = etree.iterparse(
            xml,
            events=('end',),
            tag='{http://www.opengis.net/gml/3.2}Curve',
            recover=True
        )

iterparseにおける、この引数は、lxml==3.4.1以降に導入されているため、バージョンを指定してlxmlをインストールする必要があります。

easy_install lxml==3.4.1

以上のことを踏まえて、鉄道データのXMLをデータベースにインポートする処理は以下のようになります。

railway_db.py

# -*- coding: utf-8 -*-
import sqlite3
import sys
import os
# easy_install lxml==3.4.1
from lxml import etree
from peewee import *

database_proxy = Proxy()
database = None

class BaseModel(Model):
    """
    モデルクラスのベース
    """
    class Meta:
        database = database_proxy

class Curve(BaseModel):
    """
    曲線情報モデル
    """
    curve_id = CharField(index=True, unique=False)
    lat = DoubleField()
    lng = DoubleField()

class RailRoadSection(BaseModel):
    """
    鉄道区間情報モデル
    """
    gml_id = CharField(primary_key=True)
    # 外部キーは主キーまたは一意制約をもっていないといけないので、
    # 複数あるデータに関しては外部キーとしては指定できない。
    location = CharField(index=True)
    railway_type = IntegerField()
    service_provider_type = IntegerField()
    railway_line_name = CharField(index=True)
    operation_company = CharField(index=True)

class Station(BaseModel):
    """
    駅情報モデル
    """
    gml_id = CharField(primary_key=True)
    # 外部キーは主キーまたは一意制約をもっていないといけないので、
    # 複数あるデータに関しては外部キーとしては指定できない。
    location = CharField(index=True)
    railway_type = IntegerField()
    service_provider_type = IntegerField()
    railway_line_name = CharField(index=True)
    operation_company = CharField(index=True)
    station_name = CharField(index=True)
    railroad_section = ForeignKeyField(
        db_column='railroad_section_id',
        rel_model=RailRoadSection,
        to_field='gml_id',
        index=True
    )

def setup(path):
    """
    データベースのセットアップ
    @param path データベースのパス
    """
    global database
    database = SqliteDatabase(path)
    database_proxy.initialize(database)
    database.create_tables([Curve, RailRoadSection, Station], True)

def import_railway(xml):
    """
    国土数値院のN02-XX.xmlから路線と駅情報をインポートする
    TODO:
      外部キーがらみのインポートの効率がわるい
    @param xml XMLのパス
    """
    commit_cnt = 2000  # ここで指定した数毎INSERTする
    f = None
    contents = None
    namespaces = {
        'ksj': 'http://nlftp.mlit.go.jp/ksj/schemas/ksj-app',
        'gml': 'http://www.opengis.net/gml/3.2',
        'xlink': 'http://www.w3.org/1999/xlink',
        'xsi': 'http://www.w3.org/2001/XMLSchema-instance'
    }

    with database.transaction():
        insert_buff = []
        context = etree.iterparse(
            xml,
            events=('end',),
            tag='{http://www.opengis.net/gml/3.2}Curve',
            recover=True
        )
        for event, curve in context:
            curveId = curve.get('{http://www.opengis.net/gml/3.2}id')
            print (curveId)
            posLists = curve.xpath('.//gml:posList', namespaces=namespaces)
            for posList in posLists:
                points = posList.text.split("\n")
                for point in points:
                    pt = point.strip().split(' ')
                    if len(pt) != 2:
                        continue
                    insert_buff.append({
                        'curve_id': curveId,
                        'lat': float(pt[0]),
                        'lng': float(pt[1])
                    })
                    if len(insert_buff) >= commit_cnt:
                        Curve.insert_many(insert_buff).execute()
                        insert_buff = []
        if len(insert_buff):
            Curve.insert_many(insert_buff).execute()
        insert_buff = []
        context = etree.iterparse(
            xml,
            events=('end',),
            tag='{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}RailroadSection',
            recover=True
        )
        for event, railroad in context:
            railroadSectionId = railroad.get(
                '{http://www.opengis.net/gml/3.2}id'
            )
            locationId = railroad.find(
                'ksj:location',
                namespaces=namespaces
            ).get('{http://www.w3.org/1999/xlink}href')[1:]
            railwayType = railroad.find(
                'ksj:railwayType', namespaces=namespaces
            ).text
            serviceProviderType = railroad.find(
                'ksj:serviceProviderType',
                namespaces=namespaces
            ).text
            railwayLineName = railroad.find(
                'ksj:railwayLineName',
                namespaces=namespaces
            ).text
            operationCompany = railroad.find(
                'ksj:operationCompany',
                namespaces=namespaces
            ).text
            insert_buff.append({
                'gml_id': railroadSectionId,
                'location': locationId,
                'railway_type': railwayType,
                'service_provider_type': serviceProviderType,
                'railway_line_name': railwayLineName,
                'operation_company': operationCompany
            })
            print (railroadSectionId)
            if len(insert_buff) >= commit_cnt:
                RailRoadSection.insert_many(insert_buff).execute()
                insert_buff = []
        if len(insert_buff):
            RailRoadSection.insert_many(insert_buff).execute()

        insert_buff = []
        context = etree.iterparse(
            xml,
            events=('end',),
            tag='{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}Station',
            recover=True
        )
        for event, railroad in context:
            stationId = railroad.get('{http://www.opengis.net/gml/3.2}id')
            locationId = railroad.find(
                'ksj:location', namespaces=namespaces
            ).get('{http://www.w3.org/1999/xlink}href')[1:]
            railwayType = railroad.find(
                'ksj:railwayType',
                namespaces=namespaces
            ).text
            serviceProviderType = railroad.find(
                'ksj:serviceProviderType',
                namespaces=namespaces
            ).text
            railwayLineName = railroad.find(
                'ksj:railwayLineName',
                namespaces=namespaces
            ).text
            operationCompany = railroad.find(
                'ksj:operationCompany',
                namespaces=namespaces
            ).text
            stationName = railroad.find(
                'ksj:stationName',
                namespaces=namespaces
            ).text
            railroadSection = railroad.find(
                'ksj:railroadSection',
                namespaces=namespaces
            ).get('{http://www.w3.org/1999/xlink}href')[1:]
            print (stationId)
            insert_buff.append({
                'gml_id': stationId,
                'location': locationId,
                'railway_type': railwayType,
                'service_provider_type': serviceProviderType,
                'railway_line_name': railwayLineName,
                'operation_company': operationCompany,
                'station_name': stationName,
                'railroad_section': RailRoadSection.get(
                    RailRoadSection.gml_id == railroadSection
                )
            })
            if len(insert_buff) >= commit_cnt:
                Station.insert_many(insert_buff).execute()
                insert_buff = []
        if len(insert_buff):
            Station.insert_many(insert_buff).execute()

データベースに格納してしまえば、あとは簡単に活用できます。

使用上の注意

国土数値情報(鉄道データ)を取り扱う場合に気付いた点を以下に記述します。

・路線名だけでは絞りこめない。
 たとえば、「1号線」といった場合に、「横浜市」が保持する場合と、「千葉モノレール」が保持する場合がある。
 そのため、「運行会社」と「路線名」で絞りこむ必要がある。

・いつも利用している名称とことなる名称の場合がある。
 JR東日本は東日本旅客鉄道になるし、東京メトロは東京地下鉄になっている。

・いつも利用している路線と違う路線になっている場合がある。
 たとえば、通常の路線図では「東京」は「中央線」に含まれている。
 しかし、国土数値情報としては「東京」は「中央線」に含まれていない。
 「東京」~「神田」は「東北線」とみなされる。
 これは、東京駅 - 神田駅間は東北本線上に敷設された専用線路を走行しているためだと思われる。

VBAのFixやIntの計算誤差は浮動小数点レジスタの精度がかかわっている

問題

Double型のデータをIntやFixで切り捨てた場合に誤差がでます。
浮動小数点なので当然誤差は出ますが、CDblのキャストのタイミングでも誤差が発生します。

ExcelVBAのイミデイトウィンドウで下記を実行してください。

?Fix(0.6# * 10 )
5
?Fix(CDbl(0.6# * 10))
6

環境によってはFix関数に渡す前にCDblでキャストした場合と、しない場合で計算結果が変わります。

再現環境
 Office2010 32bit
 Windows7 64biy
 CPU Intel(R) Core(TM) i5 CPU M450 2.4GHz

原因

この問題は、下記のページでも議論されています。(このページではVB6が話題ですが、基本的に同じです)

http://hanatyan.sakura.ne.jp/logbbs1/wforum.cgi?mode=allread&no=3276&page=0#3276

この原因には、計算途中の浮動小数点の精度とDouble型の精度に違いがあるために発生しています。

浮動小数点の演算の途中経過が64ビットの仮数を持つ80ビットのレジスタで計算しているが、Dobule型は53ビットの仮数を持つ64bitの浮動小数点なので、結果として誤差が表示してしまいます。

「VC++とかで同じようなコードを書いても再現しないじゃないか!やっぱりVBAは壊れているじゃないか!(憤怒)」 と怒る兄貴、姉貴もおられると思いますが、浮動小数点レジスタの仮数を64ビットの仮数で計算するか、53ビットで計算するかはプログラム中で変更することができます。

このことを確認するには、_controlfp_sの関数を実行して調べることができます。

_controlfp_s
http://msdn.microsoft.com/ja-jp/library/c9676k6h%28v=vs.80%29.aspx

この関数から取得できる値に_MCW_PCをマスクした結果で浮動小数点の精度が確認できます。

_MCW_PC (精度制御)

Mask 定数 value
0x00030000 _PC_24 (24 ビット)
_PC_53 (53 ビット)
_PC_64(64 ビット)
0x00020000
0x00010000
0x00000000

つまり、VC++のデフォルトでは、_MCW_PCの_PC_53がたっている状態であり、53ビットの仮数を持つ64bitの浮動小数点になっているため計算誤差がでないのです。

VBAでの浮動小数点レジスタの確認方法

残念ながら、VBAで直接、_controlfp_s関数を呼ぶことはできません。
しかし、VBAではDLLを実行できるので、DLL経由で_controlfp_sを実行することで、確認することができます。

まず、VC++でDLL側のコードを記述します。

DLL側のコード

# include <windows.h>
# include <float.h>
int __stdcall GetControlFP() 
{ 
  int err,sts;
  err=_controlfp_s(&sts,0,0 );
  return sts;
} 
int WINAPI DllMain(HINSTANCE hInst, DWORD fdwReason, PVOID pvReserved) 
{ 
    return TRUE; 
} 

VisualStdio2008SP1 で32bitでコンパイルした結果
http://needtec.sakura.ne.jp/vba/checkfloat/VBADll.zip

※もし、64bitのVBAを使用している場合、64bitでビルドして64bitのDLLを作成する必要があります。

作成したDLLをVBAから呼び出します。

Option Explicit
' パスを修正すること
Declare Function GetControlFP Lib "C:\Users\xxxxx\Documents\Visual Studio 2008\Projects\VBADll\Release\VBADll.dll" () As Long

Public Sub TestCon()
    '
    '  _MCW_PC (精度制御) 0x00030000
    '
    ' _PC_24 (24 ビット) 0x00020000
    ' _PC_53 (53 ビット) 0x00010000
    ' _PC_64 (64 ビット) 0x00000000
    Msgbox hex$(GetControlFP())
End Sub

先にあげた検証環境では「C001F」という結果がかえってきました。
これにより、精度制御が_PC_64ビットであり、浮動小数点の演算の途中経過が64ビットの仮数を持つ80ビットのレジスタで計算していることが確認できました。

じゃあ、VBAたんは無罪なのね!

・・・じつは、参照しているCOMやDLLなどにより変更される可能性があり、この結果は全ての状況で正しいとは限りません。

たとえば、Jetを使うと途中の精度がかわったりするバグがあったりもします。
http://support.microsoft.com/kb/308702/ja

・・・やっぱり壊れているじゃないか!(憤怒)