JavaScriptで空間演算

2015/03/04
先日FirefoxとQGISで行った位置情報の取得・加工を、JavaScriptで行ってみる。
スクレイピングにはnode-horsemanを使用した。Nightmareによく似た使いやすいモジュール。結果をGeoJSON形式で出力する。1ページ5秒のウエイトと、結果に5000件の制限を設けてある。
var fs = require('fs');

var Horseman = require('node-horseman');
var horseman = new Horseman();

var results = [];

function getResult() {
    return horseman
        .evaluate(function() {
            var features = [];
            $('table[id="searchResult"] tbody tr:gt(1)').each(function(item) {
                var lat = $.trim($(this).find('td:eq(4)').text()).split("°");
                var lat_deg = parseFloat(lat[0]);
                lat = lat[1].split("′");
                var lat_min = parseFloat(lat[0]) / 60;
                lat = lat[1].split("″");
                var lat_sec = parseFloat(lat[0]) / 3600;
                var lng = $.trim($(this).find('td:eq(5)').text()).split("°");
                var lng_deg = parseFloat(lng[0]);
                lng = lng[1].split("′");
                var lng_min = parseFloat(lng[0]) / 60;
                lng = lng[1].split("″");
                var lng_sec = parseFloat(lng[0]) / 3600;

                var feature = {
                    'type': 'Feature',
                    'geometry': {
                        'type': 'Point',
                        'coordinates': [
                            (lng_deg + lng_min + lng_sec),
                            (lat_deg + lat_min + lat_sec)
                        ]
                    },
                    'properties': {
                        'id' :            $.trim($(this).find('td:eq(0) a').text()),
                        'project' :       $.trim($(this).find('td:eq(1)').text()),
                        'investigation' : $.trim($(this).find('td:eq(2)').text()),
                        'organization' :  $.trim($(this).find('td:eq(3)').text()),
                        'length' :        $.trim($(this).find('td:eq(6)').text()),
                        'elevation' :     $.trim($(this).find('td:eq(7)').text()),
                        'fig' :           $.trim($(this).find('td:eq(8) a').attr("href")),
                        'result' :        $.trim($(this).find('td:eq(9) a').attr("href")),
                        'xml' :           $.trim($(this).find('td:eq(0) a').attr("href"))
                    }
                };
                features.push(feature);
            });
            return features;
        });
}

function hasNextPage() {
    return horseman.exists('td.nextLink a');
}

function scrape() {
    var result = getResult();
    results = results.concat(result);
    var page = horseman
        .evaluate(function() {
            return $('td.pageCount').text();
        });
    console.log(page);
    if (hasNextPage() && results.length < 5000){
        horseman
            .click('td.nextLink a')
            .waitForNextPage()
            .wait(5000);
        scrape();
    }
}

horseman
    .userAgent("Mozilla/5.0 (Windows NT 6.3; WOW64; rv:35.0) Gecko/20100101 Firefox/35.0")
    .viewport(1024, 768)
    .open('http://www.kunijiban.pwri.go.jp/jp/denshikokudo/DB_Search/boringsearch.php')
    .type('input[name="tl_lat_deg"]', '34')
    .type('input[name="tl_lat_min"]', '23')
    .type('input[name="tl_lat_sec"]', '52')
    .type('input[name="tl_lon_deg"]', '131')
    .type('input[name="tl_lon_min"]', '54')
    .type('input[name="tl_lon_sec"]', '03')

    .type('input[name="br_lat_deg"]', '32')
    .type('input[name="br_lat_min"]', '51')
    .type('input[name="br_lat_sec"]', '40')
    .type('input[name="br_lon_deg"]', '133')
    .type('input[name="br_lon_min"]', '52')
    .type('input[name="br_lon_sec"]', '03')
    .click('input[value="検索"]')
    .waitForSelector('td.nextLink a')
//    .screenshot('kunijiban.png');

scrape();

var GeoJSON = {
    "type": "FeatureCollection",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:EPSG::4612"
        }
    },
    "features": results
};

fs.writeFile('kunijiban.geojson', JSON.stringify(GeoJSON, null, "    "));

horseman.close();

次に行政区域の境界のバッファを作成して、空間演算により取得した検索結果のうち、愛媛県のみを抽出する。

行政区域の地形データはQGISで作成してもよいが、下記のGeoJSONを使用した。

dataofjapan/land


空間演算には下記のモジュールを使用した。

bjornharrtell/jsts
JSTS Example

下記のhtmlをローカルに置いてブラウザで読みこめば、同じディレクトリにあるjapan.geojson, kunijiban.geojsonを読み込んで地図上に表示する。ただしFirefoxとChromeのみ(Chromeはローカルファイルが読めるように、"--allow-file-access-from-files"オプションを付けて起動する必要がある)。空間演算が割と重くて時間がかかるので、ハングしたようになるかも。
<!DOCTYPE html>
<html>

<head>
    <title>ボーリング柱状図</title>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">

    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css" />

    <script src="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.js"></script>
    <script src="http://code.jquery.com/jquery-2.1.0.min.js"></script>
    <script src="https://rawgithub.com/bjornharrtell/jsts/master/lib/javascript.util.js"></script>
    <script src="https://rawgithub.com/bjornharrtell/jsts/master/lib/jsts.js"></script>

    <link rel="kunijiban" type="application/json" href="kunijiban.geojson">
    <link rel="japan" type="application/json" href="japan.geojson">

    <style>
        html,
        body,
        #map {
            height: 100%;
            width: 100%;
            padding: 0px;
            margin: 0px;
        }
    </style>
</head>

<body>
    <div id="map"></div>

    <script>
        var map = L.map('map').setView([33.523342, 132.863776], 10);
        L.control.scale({'position':'bottomleft','metric':true,'imperial':false}).addTo(map);

        var distance = 5 / 111.12; // 5km
        var geoReader = new jsts.io.GeoJSONReader(),
            geoWriter = new jsts.io.GeoJSONWriter();

        var cyber = L.tileLayer('http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png', {
            maxZoom: 18,
            attribution: '<a href="http://www.gsi.go.jp/kikakuchousei/kikakuchousei40182.html" target="_blank">国土地理院</a>'
        }).addTo(map);
        map.addLayer(cyber);

        var gbank = L.tileLayer.wms("https://gbank.gsj.jp/ows/seamlessgeology200k_b", {
            layers: ['area,line,label'],
            format: 'image/png',
            transparent: true,
            opacity: 0.5,
            attribution: '<a href="https://gbank.gsj.jp/seamless/" target="_blank">シームレス地質図</a>'
        }).addTo(map);
        map.addLayer(gbank);

        L.control.layers({'地理院地図': cyber}, {'シームレス地質図': gbank}).addTo(map);

        function onEachFeature(feature, layer) {
            var popupContent = '<p><a href="' + feature.properties.fig + '" target="_blank">';
            popupContent += feature.properties.id + '</a><br />';
            popupContent += feature.properties.investigation + '</p>';

            layer.bindPopup(popupContent);
        }

        function createRegionBuffer() {
            return new Promise(function (resolve, reject) {
                $.getJSON($('link[rel="japan"]').attr("href"), function(data) {
                    var district;
                    data.features.forEach(function(item, index) {
                        if (item.properties.nam_ja == '愛媛県') district = item;
                    });
                    resolve(geoReader.read(district.geometry).buffer(distance));
                });
            });
        }

        function readBoringData(region) {
            return new Promise(function (resolve, reject) {
                $.getJSON($('link[rel="kunijiban"]').attr("href"), function(data) {
                    var inRegion = data.features.filter(function(item, index) {
                        return geoReader.read(item.geometry).within(region);
                    });
                    resolve(inRegion);
                });
            });
        }

        function renderRegionBuffer(region) {
            return new Promise(function (resolve, reject) {
                var buffer = geoWriter.write(region);
                L.geoJson(buffer, {
                    style: {
                        weight: 2,
                        color: "#999",
                        opacity: 1,
                        fillColor: "#B0DE5C",
                        fillOpacity: 0.8
                    }
                }).addTo(map);
                resolve(region);
            });
        }

        function renderBoringData(points) {
            L.geoJson(points, {
                onEachFeature: onEachFeature,

                pointToLayer: function(feature, latlng) {
                    return L.circleMarker(latlng, {
                        radius: 8,
                        fillColor: "#ff7800",
                        color: "#000",
                        weight: 1,
                        opacity: 1,
                        fillOpacity: 0.8
                    });
                }
            }).addTo(map);
        }

        createRegionBuffer()
//        .then(renderRegionBuffer)
        .then(readBoringData)
        .then(renderBoringData);
    </script>
</body>

</html>

マップは国土地理院のタイルと、日本シームレス地質図を使用している。

Read more ...