Q: MongoDB group by geospatial location I am looking to group geospatial data in MongoDB using an aggregate query. I have a collection containing around a million documents, with geojson with a 2D Sphere index. For a given radius or polygon, I would like to group and count objects by n-meter squares, returning a "mesh" of groups that would cover my given radius.
I am able to loop through predetermined polygons using a $geoNear query (each square of the mesh one at a time) one at a time to get the result I'm looking for, but it is inefficient. I would like to somehow query once and group the results by lat/long with some tolerance.
Something like this would show in my ui (except with counts)
https://www.researchgate.net/profile/Luca-Bedogni/publication/317988792/figure/fig1/AS:652923713376258@1532680552765/MGRS-encoding-of-Northern-Italy-Figure-2a-shows-the-subdivision-of-GZIs-in-100-km.png
Not like this (uneven division of points)
https://docs.mongodb.com/charts/saas/images/charts/geospatial-scatter-coordinates-example.png
It seems like this type of information is stored within the MongoDB geospatial index itself, but I don't see a way to access it.  Any help with aggregate query strategy would be greatly appreciated...is it even possible?
Simply, can I query a collection of geojson Points using a 10000m square, then group/count the results based on 100x100 squares?

A: The trick is to use recursion when constructing the bounding cell to pass to the $geoWithin function.   Here is what happens at a high level:

*

*Make a big (lon,lat) square of points and call it coords.  You at level 0.

*Perform the pipeline.

*If anything comes back, then level++ and divide the square into 4 smaller ones and recursively call the pipeline for topleft, topright, bottomleft, and bottomright squares.

"That's it."  This approach work very nicely when you have a combination of small cells and concentrations (i.e. not uniformly distributed) intersections with candidate geometries.  Very often dozens of small cell groupings can be bypassed because the very quickly it can be determined at a level closer to 0 that vast ranges of geo are empty.
Below is the code.  Note that I save the results of the level in a collection named sqrbin because the use case demanded zoom in/out analysis and needed VERY rapid query response.   It does not exactly match the OPs question but it is close.
/*

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");



An alternate approach is to use the JTS geo library in Java (https://locationtech.github.io) or similar and use MongoDB only to fetch the big outer region and perform no grouping at all, then use pure code to iterate thru the 100s or 1000s of cells against the results.  This is all in-memory and should take only milliseconds.  It works well when the one $geoWithin quickly yields << than the total collection count; in this case, say 50000 or less (better than an order of magnitude and heading toward 2).
UPDATED
There is one more "easy approach" if you wish to group count within a square: do math directly on lon and lat.  The idea is that given a top left corner and an extent in degrees (e.g. .1 is 11km) and the "matrix resolution" (e.g. 10 means a 10x10 grid inside that square, then after the main $match to extract points in the square you compute the row and col based on the formula (delta pt to  corner)/resolution.  This approach performs even with a resolution of 100 (i.e. a 100x100 matrix) which will be implicity faster than 10,000 discrete calls to $geoWithin.
// 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                                              
]);

