LoginSignup
4
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@kkdd

R-tree 矩形検索（円筒図法地理座標）

Last updated at Posted at 2017-06-06

こんにちは。
R-tree (rectangle tree) を利用し、矩形形状を持つデータの集合に対して指定矩形範囲内に含まれるものを検索してみました（地理座標系：下記図は Braun投影図法）[Demonstration]。マウスでこの矩形検索範囲を変えられます（d3.js 利用1 2）。

ここで使っている r-tree.js は、"JavaScript (TypeScript) implementation of a Hilbert Packed R-Tree"3 を基にして機能拡張して作ったものです（日付変更線への対処など4）。利用方法を要約すると、

r-tree_usage.js
var tree = new RTree(maxNodes);
tree.batchInsert(rectsData);
var options = {xPeriod: 360, includedOnly: true, searchIndex: true}; // coping with antimeridian-crossing case
var resultIndexes = tree.search(boudingRect, options); // returns the indexes of data array

braun2.jpg

r-tree.html
<!DOCTYPE html>
<meta charset="utf-8">
<title>Braun projection</title>
<!--demonstration: https://bl.ocks.org/kkdd/92a185b33e387aa1de19cbce2ba7464f -->
<style>
.frame {
  fill: none;
  stroke: #000;
  stroke-width: 2px;
}

.fill {
  fill: #fff;
}

.graticule {
  fill: none;
  stroke: #aaa;
  stroke-width: 0.5px;
  stroke-opacity: 0.5;
}

.land {
  fill: #ddd;
}

.boundary {
  fill: none;
  stroke: #fff;
  stroke-width: 0.5px;
}

.rect {
  fill: #8f0;
  fill-opacity: .25;
  stroke: #0f0;
  stroke-width: 0.4px;
  pointer-events: none;
}

.rects {
  fill: #44f;
  fill-opacity: .6;
  stroke: #00f;
  stroke-width: 0.2px;
  pointer-events: none;
}
</style>
<div id="canvas" width="600" height="400"></div>
<script src="https://d3js.org/d3.v4.min.js"></script>
<script src="https://unpkg.com/topojson-client@3"></script>
<script src="./r-tree.js"></script>
<!-- https://github.com/kkdd/javascript-hilbert-r-tree/tree/master/r-tree -->
<script>
var rect, points = [[90, -10], [140, 30]]; // bounding rectangle
var rects, rectsData = [
  {x: 105, y: -45, width: 40, height: 20},
  {x: 155, y: -35, width: 30, height: 40},
  {x: 145, y:  25, width: 40, height: 30},
  {x: 115, y: -15, width: 20, height: 30},
  {x: 105, y: 55, width: 20, height: 10},
  {x: 15,  y: -15, width: 30, height: 30},
  {x: -165, y: -25, width: 60, height: 20},
  {x: -125, y: 35, width: 80, height: 10},
  {x: -155, y: 55, width: 20, height: 10},
  {x: -155, y: 15, width: 20, height: 20},
  {x: -95, y: -15, width: 50, height: 20},
];

for (var i=0;i<1000; i++) {
  var pt = [Math.random()*2-1, Math.random()*2-1];
  if (pt[0]*pt[0]+pt[1]*pt[1] < 1) {
	  rectsData.push({x: Math.floor(pt[0]*70+30), y: Math.floor(pt[1]*70), width: Math.floor(Math.random()*6+2), height: Math.floor(Math.random()*6+2)});
  }
}

var maxNodes = 4;   // Maximum number of nodes within a bounding rectangle
var tree     = new RTree(maxNodes);
tree.batchInsert(rectsData);

var svg = d3.select("#canvas")
  .append("svg")
  .attr("width", 600)
  .attr("height", 400);

function braunRaw(lambda, phi) {
  return [lambda, Math.tan(phi/2)*2];
}
braunRaw.invert = function(x, y) {
  return [x, Math.atan(y/2)*2];
};

var width = +svg.attr("width"), height = +svg.attr("height");
var projection = d3.geoProjection(braunRaw)
    .scale((width - 3) / (2 * Math.PI))
    .translate([width/2, height/2])
    .precision(0.1);

var path = d3.geoPath()
    .projection(projection);

var graticule = d3.geoGraticule();

svg.append("defs").append("path")
    .datum(graticule.outline())
    .attr("id", "sphere")
    .attr("d", path);

svg.append("use")
    .attr("class", "frame")
    .attr("xlink:href", "#sphere");

svg.append("use")
    .attr("class", "fill")
    .attr("xlink:href", "#sphere");

svg.append("path")
    .datum(graticule)
    .attr("class", "graticule")
    .attr("d", path);

var drag = d3.drag()
    .on("drag", function (d, i) {
      var pt = [d3.event.x, d3.event.y];
      points[i] = projection.invert(pt);
      drawGeom();
      d3.select(this)
        .attr("cx", d.x = pt[0])
        .attr("cy", d.y = pt[1]);
    });

svg.selectAll("circle")
	.data(points).enter()
	.append("circle")
	.attr("cx", function (d) {return projection(d)[0]})
	.attr("cy", function (d) {return projection(d)[1]})
	.attr("r", "2.5px")
	.attr("stroke-width", "1px")
	.attr("stroke", function(d, i) { return ["#00f", "#f00"][i]; })
	.attr("fill", "#fff")
    .call(drag);

drawGeom();

function drawGeom() {
    if (rect) {rect.remove()}
    if (rects) {rects.remove()}

    rect = svg.append("path")
		.datum(boundRect(points))
		.attr("class", "rect")
		.attr("d", path);

    var options = {xPeriod: 360, includedOnly: true, searchIndex: true}; // coping with antimeridian-crossing case (overlapping/intersecting rectangles)
    var resultIndexes = tree.search(bound2rect(points), options);	
    var resultRects = resultIndexes.map(function (i) {
        return boundRect(rect2bound(rectsData[i]));
    });

    rects = svg.selectAll("rects")
		.data(resultRects).enter()
		.append("path")
		.attr("class", "rects")
		.attr("d", path);
}

function rect2bound(rect) {
	return [[rect.x, rect.y], [rect.x+rect.width, rect.y+rect.height]];
}

function bound2rect(bb) {
	var height = bb[1][1] - bb[0][1],
	    width = moduloPascal(bb[1][0] - bb[0][0], 360);
	return {x: bb[0][0], y: bb[0][1], width: width, height: height};
}

function moduloPascal(a, b) {
	return (a%b+b)%b;
}

d3.json("https://unpkg.com/world-atlas@1/world/50m.json", function(error, world) {
  svg.insert("path", ".graticule")
      .datum(topojson.feature(world, world.objects.land))
      .attr("class", "land")
      .attr("d", path);

  svg.insert("path", ".graticule")
      .datum(topojson.mesh(world, world.objects.countries, function(a, b) { return a !== b; }))
      .attr("class", "boundary")
      .attr("d", path);
});

function boundRect(pnts) {
  var coords = [[pnts[0][0],pnts[0][1]]]
      .concat(parallel(pnts[0][0], pnts[1][0], pnts[0][1]))
      .concat(parallel(pnts[0][0], pnts[1][0], pnts[1][1]).reverse());
  return rfc7946tod3({type: "Polygon", coordinates: [coords]})
}

function parallel(λ0, λ1, φ)  {
  if (λ0 > λ1) λ1 += 360;
  var dλ = λ1 - λ0,
      step = dλ / Math.ceil(dλ);
  return d3.range(λ0, λ1 + 0.5 * step, step).map(function(λ) { return [normalise(λ), φ]; });
}

function normalise(x) {
  return (x + 180) % 360 - 180;
}

// https://github.com/tyrasd/rfc7946-to-d3
function rfc7946tod3(geojson) {
    switch ((geojson && geojson.type) || null) {
        case 'FeatureCollection':
            geojson.features.forEach(rfc7946tod3)
            break
        case 'GeometryCollection':
            geojson.geometries.forEach(rfc7946tod3)
            break
        case 'Feature':
            rfc7946tod3(geojson.geometry)
            break
        case 'MultiPolygon':
            geojson.coordinates.forEach(function(polygon) {
                polygon.forEach(function(ring) {
                    ring.reverse()
                })
            })
            break
        case 'Polygon':
            geojson.coordinates.forEach(function(ring) {
                ring.reverse()
            })
            break
    }
    return geojson
}
</script>

  1. d3.js の bounding rectangle について、"Geographic Bounding Boxes" に解説があります。 ↩

  2. 内部で、gojson の polygon 処理に、rfc7946-to-d3.js を使っています。 ↩

  3. なお、Javascript 版 R-tree としては、RBush が有名のようです。 ↩

  4. GIS 処理系の中には正しい検索ができない（日付変更線への対処など） ものもあるようです。 ↩

4
4
0

Register as a new user and use Qiita more conveniently

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

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?