诀窍是在构造边界单元时使用递归传递给$geoWithin 函数。以下是在高层次上发生的事情：
- 制作一个大的 (lon,lat) 正方形点并将其命名为
coords。你在 0 级。
- 执行管道。
- 如果有任何返回，然后
level++ 并将正方形分成 4 个较小的正方形，并递归调用 topleft、topright、bottomleft 和 bottomright 正方形的管道。
“就是这样。”当您将小细胞 和 浓度（即非均匀分布）与候选几何形状相交时，这种方法非常有效。很多时候，可以绕过几十个小小区分组，因为可以很快确定在接近 0 的水平上，大范围的地理是空的。
下面是代码。请注意，我将级别的结果保存在名为sqrbin 的集合中，因为用例需要放大/缩小分析并且需要非常快速的查询响应。它与 OPs 问题不完全匹配，但很接近。
/*
The Recursive Squarebin generator.
Algo:
Start with a square and level 0.
Capture square as NW point and a side size.
Given a NW point and level.
dive = false;
If res is too low (too big a square)
dive = true
else
run a geom for that square.
If there is material (numbers, etc.) then
post the material and level (e.g. 0)
dive = true;
if dive == true and level < maxDepth then
break up the square into 4 subsquares
recurse into each of the subsquares (top left, top right,
bottom left, bottom right) with level+1
That's it!
*/
db = db.getSiblingDB("myDBWithGeoCollections");
/*
units place in lon,lat is good for 111km
first decimal place in lon,lat is good for 11.1km
TX is 1200km
8 sqrbins ; so ... 150km units (eg. 90 to 91 is 150km)
18 sqrbin ; so ... 66km "half units?"
top left (canadian border & pacific)
-125.203941, 49.143621
top right
-66.928408, 49.143621
bottom left:
-125.203914, 25.85675
bottom right:
-66.928408, 25.856750
*/
//var maxN = 100; // safety switch for now
var maxN = -1; // -1 means let 'er rip!
var emptyTotal = 0;
function geoMe(nwpt, res, level) {
var coords = [];
coords.push( [ nwpt[0], nwpt[1] ] );
coords.push( [ nwpt[0]+res, nwpt[1] ] );
coords.push( [ nwpt[0]+res, nwpt[1]-res ] );
coords.push( [ nwpt[0], nwpt[1]-res ] );
coords.push( [ nwpt[0], nwpt[1] ] ); // close loop
c = db.site.aggregate([
{$match: { "loc": { $geoWithin: { $geometry:
{ type: "Polygon", coordinates: [ coords ] } }}
}}
,{$group: {_id:null, n: {$sum:1}}}
]);
var foundSomething = false;
// Even if the geo yields nothing, you STILL get 1 record back
// in the group where n = 0 so that's how we can tell...
c.forEach(function(d2) {
// ...but ONLY create things where tiv or totGrossLoss > 0
if(d2['n'] > 0) {
d2['loc'] = {
type: "Polygon",
coordinates: [ coords ]
};
d2['res'] = res;
d2['level'] = level;
db.sqrbin.insert(d2);
foundSomething = true;
maxN--;
} else {
emptyTotal++;
}
});
return foundSomething;
};
function doSquare(nwpt, res, level) {
var dive = false;
if(maxN == 0) {
return; // bail!
}
var nlevel = level;
if(res > resMin) {
print(res + " too high; dive");
dive = true;
// Do not incr level if res is too low.
// This permits level 0 to be the first one with material
} else {
// Try to do something at this level
dive = geoMe(nwpt, res, level);
nlevel = level+1; // ah!
}
if(dive == true && res > resMax) {
var nres = res/2.0; /// ah HA!
var nnwpt = [nwpt[0], nwpt[1]];
doSquare(nnwpt, nres, nlevel); // top left;
var nnwpt = [nwpt[0]+nres, nwpt[1]];
doSquare(nnwpt, nres, nlevel); // top right;
var nnwpt = [nwpt[0], nwpt[1]-nres];
doSquare(nnwpt, nres, nlevel); // bottom left
var nnwpt = [nwpt[0]+nres, nwpt[1]-nres];
doSquare(nnwpt, nres, nlevel); // bottom right
}
}
/* This is boundbox of USA minus HI:
var llon = -125.203941;
var blat = 25.040;
var rlon = -66.928408;
var tlat = 49.143621;
So -125 - -67 is 58. Call it 60 to be safe
Call it 60.
*/
var nwpt = [ -125.203941, 49.143621 ];
var res = 60;
// Don't create squares smaller than this on a side.
// 0.01 is good to 1km, or 3200 ft each side.
var resMax = 0.01;
// Don't create squares bigger than this on a side.
// 1.0 is good for 110 or about 62 miles on a side.
var resMin = 1.0;
var total_start = (new ISODate()).getTime();
db.sqrbin.drop();
doSquare(nwpt, res, 0);
var total_end = (new ISODate()).getTime();
print((total_end - total_start) + " overall millis");
print(emptyTotal + " squares searched and empty");
另一种方法是使用 Java (https://locationtech.github.io) 或类似语言中的 JTS 地理库，并仅使用 MongoDB 来获取外部大区域并且根本不执行任何分组，然后使用纯代码迭代 100 或 1000对结果的细胞。这全都在内存中，应该只需要几毫秒。当$geoWithin 快速产生
更新
如果您希望在正方形内对计数进行分组，还有另一种“简单方法”：直接在 lon 和 lat 上进行数学运算。这个想法是，给定左上角和度数（例如 .1 是 11 公里）和“矩阵分辨率”（例如 10 表示该正方形内的 10x10 网格，然后在主要 $match 之后提取点根据公式 (delta pt to corner)/resolution 计算行和列的平方。即使分辨率为 100（即 100x100 矩阵），这种方法也能执行，这将比对 $geoWithin 的 10,000 次离散调用更快。
// Inputs. leftLon,topLat is the top left point of our square:
leftLon = -76.487605;
topLat = 41.0215428;
degrees = .1; // Every .1 is 11km
rez = 10; // Make a 10x10 matrix inside the bounding square
// From here down all based on inputs.
// Create a bounding square for our query:
squareGeom = [
[leftLon, topLat],
[leftLon, topLat - degrees],
[leftLon + degrees, topLat - degrees],
[leftLon + degrees, topLat],
[leftLon, topLat] // close the loop
];
factor = degrees/rez;
c = coll.aggregate([
{'$match': {'loc': {'$geoWithin': {'$geometry': {'coordinates': [ squareGeom ], 'type': 'Polygon'}}}
}}
,{$addFields: {
row: {$floor:{$divide: [{$subtract: [ topLat, {$arrayElemAt: ["$loc.coordinates",1]}]}, factor]}},
col: {$floor:{$divide: [{$subtract: [ {$arrayElemAt: ["$loc.coordinates",0]}, leftLon]}, factor]}},
}}
,{$group: {_id: {row: "$row", col: "$col"}, n: {$sum:1}}}
,{$sort: {"_id.row":1, "_id.col":1}} // optional but handy
]);