I've built a few fun demos using SpatiaLite and Datasette. Here's an API that tells you the timezone for a latitude longitude point:
https://timezones-api.datasette.io/timezones/by_point?longit... - implementation here: https://github.com/simonw/timezones-api
And here's the same thing for which US county a point is within: https://us-counties.datasette.io/counties/county_for_latitud... - implementation here: https://github.com/simonw/us-counties-datasette
I've also built an experimental Datasette plugin that lets you draw a shape on a Leaflet map to generate a GeoJSON polygon, then uses SpatiaLite to show you geometries that are contained by that drawn polygon. I wrote about that here: https://simonwillison.net/2021/Jan/24/drawing-shapes-spatial...
Last time I needed something along those lines, I implemented geohash  in the application level ontop of a time-series database. I assume there is a better, more sound way to achieve that, as geohash doesn't handle boundaries well.
For the uninitiated, geo-hash uses a space filling curve multiple times at different precision to map quadrants into strings. So imagine a map, top left is North America, bottom right Australia, we can name the quadrants 0,1,2,3, we then take each quadrant and again use a z filling curve to split it to more quadrants (well, quadrant of a quadrant) and keep the name we gave it, we repreat the process up to desired precision, finally, we take the all the names of quadrants with the same order and concatenate them, to create a string that repretents the position. Strings with the same prefix are in the same quadrant and depending on the length of the prefix, the quadrant they are in is either bigger or smaller, longer prefix, smaller quadrant.
We found that geohash, even at the application level, provided much faster queries than searching for coordinates within range on two fields.
A less generalized solution is to grid things up (quad tree and oct trees) as you implied above but the kd-tree is the generalized solution without edge cases (what if the majority of points of interest end up in a single grid square in your example?). A kd-tree is essentially what a sort list is for one dimensional data but instead supports multiple dimensions.
Geohash provides arbitrary precision, at 8 character length, we can measure difference in in the order of decameters iirc.
The edge cases are on the literal edges, where you could have two points, one at say -79.995, 35.0 and -80.001, 35.0, so the area is split at -80 and the first point has different prefix than the second even though they are very very close.
SQLite R*Tree Module - https://www.sqlite.org/rtree.html
The Geopoly Interface - https://www.sqlite.org/geopoly.html
Building QGIS extensions seems to be the only viable option I'm seeing.
It's great for basic things, but the geo index is considered a "special" index.
You can only use one "special" index at a time in a query, so (for example) you can't combine a full-text search with geo queries.
From the geopackage FAQs:
> What is the relationship between GeoPackage and SpatiaLite?
> SpatiaLite was a major influence on the vector portions of GeoPackage but after extensive
> discussions, the working group chose to diverge from the SpatiaLite format in a few subtle ways.
> SpatiaLite now supports GeoPackage as of version 4.2.0.
Spatialite is a spatial analysis database: it describes itself as "a complete and powerful Spatial DBMS (mostly OGC-SFS compliant".
They both use the same container, SQLite. And indeed there are other geospatial formats built on SQLite, such as mbtiles. But the intents are very different.