Recently I was asked about options to do ‘find nearest’ geospatial calculations for a set with lots of polygon features, and whether something like GeoSPARQL, a (spatial) database or something else like Python and GeoPandas would be a good fit. The exact question was finding the nearest top 3 Natura2000 areas for all Natura2000 areas. This came up in some of the project work Sensing Clues is doing in cooperation with their field partners.
My advise after some discussion was not to do geospatial calculations, at least not more often than necessary, meaning don’t do them runtime but run the operation upfront and store the results as properties with the data (e.g. additional triples in a graph database, some optional fields in a document oriented data platform, or just a good-old lookup table in a relational database).
I started with the Natura2000 data in Esri Shapefile format, wrote some Python code to load and process it with GeoPandas and ran the job. No optimisations, not split to run stuff in parallel or whatever. After about 2-3 days only 4000 areas were done (single process execution on my 2015 MacBook Pro), and I did not want to wait anymore so changed the approach and switched to PostgreSQL/PostGIS.
Decided to use the Natura2000 Geopackage instead of the Shapefile because running find nearest after Shapefile import in PostGIS seemed to be slower than using the Geopackage. Disclaimer: no idea yet why imported Shapefile is slower, no deep investigation done, maybe number of decimals after import, or difference in complexity of the (multi) polygon features, or the way either data from Shapefile or Geopackage are spatially indexed by PostGIS. Of course ogr2ogr helped with loading into PostGIS.
ogr2ogr -f PostgreSQL PG:"dbname='sensingclues'" /Users/emilzegers/GitHub/taatuut/NATURA9000/data/Natura2000_end2021_rev1.gpkg
Also started looping over the Natura2000 features instead of one time join all, wrote a procedure for this (prefer procedure over function because in procedure you can commit in between and track progress easier). Some features were still slow to process because they are very big and complex, and/or so are their neighbors. Some Natura2000 areas can span complete rivers and waterway systems covering tens to hundreds of square kilometers.
With this approach more than 10.000 features were done in less than a day, still with the slower ones consuming a fair amount of time. In less than two days the whole set of 27020 areas finished. A quick examination of the results revealed that the results of PostgreSQL/PostGIS and Python/Geopandas find nearest operation are the same most of the tim e. Again disclaimer that I did not check where and why they differ. But I’m sure PostGIS is ok 🙂
One optimisation that might make calculations way simpler and faster, is doing rubberband or envelope on these areas with ST_ConvexHull or ST_OrientedEnvelope, then do nearest on these ‘bounding polygons’.
In this case doing geospatial calculations only once and then add the results as properties or in a lookup table works because the set of Natura2000 areas does not change often. So redoing the whole set when needed is way cheaper then running nearest operations every time where doing analytics for specific areas. And in real life there are many more of these cases where preparing the data by creating smart lookups works fine so it can be a time and resource saver to consider.
The full description and code used can be found at https://github.com/taatuut/NATURA9000 It also includes some example code to visualise the results. Of course quick&dirty code but it should work for you too in case you want to test it.
If you find this interesting and have ideas how to improve this, or want to help Sensing Clues anyway with the work they do then let me know, and check out https://sensingclues.org/ anyway. Also looking for volunteers with GeoServer skills!