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

こんにちは。
[R-tree](https://en.wikipedia.org/wiki/R-tree) (rectangle tree) を利用し、矩形形状を持つデータの集合に対して指定矩形範囲内に含まれるものを検索してみました（地理座標系：下記図は [Braun投影図法](http://qiita.com/kkdd/items/864aa75041b77e2b4c28)）[[Demonstration](https://bl.ocks.org/kkdd/92a185b33e387aa1de19cbce2ba7464f)]。マウスでこの矩形検索範囲を変えられます（d3.js 利用[^1] [^2]）。

[^1]: d3.js の bounding rectangle について、"[Geographic Bounding Boxes](https://www.jasondavies.com/maps/bounds/)" に解説があります。
[^2]: 内部で、gojson の polygon 処理に、[rfc7946-to-d3.js](https://github.com/tyrasd/rfc7946-to-d3) を使っています。

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

[^3]: なお、Javascript 版 R-tree としては、[RBush](https://github.com/mourner/rbush) が有名のようです。
[^4]: GIS 処理系の中には正しい検索ができない（日付変更線への対処など） ものもあるようです。

```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](https://qiita-image-store.s3.amazonaws.com/0/54617/9a7e4098-ccd6-c8ff-f3ca-03f2689d6b2d.jpeg)

```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>
```

-----------------------