1

I have a huge shapefile (over 400MB) of an entire city that contains many polygons of the residential area of this city. I can check it on QGIS as the following:

enter image description here

This shapefile contains a huge number of polygons, and I'd like cut a small shapefile from it in a circular or squared area around a point. So, let's say I give a point to a system, I want to return all the polygons around this area inside an offset (just so I can send only this segmentation of the data to the frontend of my application).

I know I could convert this shapefile to GeoJSON and find what are the polygons that exist inside the area that I have created (using turf.js for example), but that would be a brute force solution and it wouldn't be efficient.

On the other hand, I know I could install PostgreSQL and PostGIS and implement solutions like this one.

My doubt is that if I can achieve the goal of finding all polygons inside a polygon efficiently using only the command line, Python or JavaScript without implementing a spatial index for that, working directly with the shapefile. Is it possible?

BERA
  • 72,339
  • 13
  • 72
  • 161
raylight
  • 1,241
  • 4
  • 19

1 Answers1

6

You can achieve this using ogr2ogr with

# SHP from Natural Earth Data
ogr2ogr out.shp ne_110m_admin_0_countries.shp -dialect SQLite -sql "SELECT * FROM ne_110m_admin_0_countries WHERE ST_intersects(ST_Buffer(ST_GeomFromText('POINT (1 43)', 4326), 10), ne_110m_admin_0_countries.geometry);"

To sort out execution with index, you may want to add after ogr2ogr the option --debug ON

PS: the buffer units are degrees here because using EPSG 4326 units. The points coordinates are longitude 1 and latitude 43. You may need to reproject if you want to use meters or directly use a layer that use a projection based on meters.

ThomasG77
  • 30,725
  • 1
  • 53
  • 93
  • Using the time command... I've realized that the time it takes for doing this operation is the same with the index file .qix generated or not (6 seconds in my case). So, creating the spatial index is necessary in this case? Is it used automatically? – raylight Apr 29 '21 at 16:07
  • It seems here the index is not used as when using SQLite dialect, can't use the shp index according to https://gis.stackexchange.com/questions/153472/using-a-shapefile-spatial-index-with-ogr2ogr-and-ogr-datasource-executesql – ThomasG77 Apr 29 '21 at 16:48
  • Made an edit. Not sure it's faster (small dataset on my test so not sure if major speed difference) – ThomasG77 Apr 29 '21 at 17:39
  • If you want to try directly with my dataset it's this one. I've been testing these two lines here. The second one is actually slower, it takes around 20 seconds while the first takes 6 seconds on my computer. However, when I use the debug mode I see that the index is being created as you said. I think the problem is that it's creating the index on the run instead of using the .qix that we've created. – raylight Apr 29 '21 at 22:31
  • So, in theory, if there's a way of creating the index before executing the second command and using it instead of creating it on the run... it could make it faster. – raylight Apr 29 '21 at 22:33
  • What is your projection? No prj in your link – ThomasG77 Apr 29 '21 at 22:45
  • Tested with GPKG (without using SHP as an intermediate) and a spatial index. Take 4 s. Not sure if I miss something. Done comparison with PostGIS 120ms with a spatial index Tested on the same machine... – ThomasG77 Apr 29 '21 at 23:30
  • Well, in case performance becomes an issue I'll use the PostGIS solution since the query of your second sample with an index is so much faster... All questions here that I find and the documentation don't seem to give any extra procedure besides creating the .qix index file with the command ogrinfo... – raylight Apr 30 '21 at 00:28
  • Quite annoying. Does not seem to be able to use the spatial index correctly using only GPKG... – ThomasG77 Apr 30 '21 at 00:32
  • Unhappy because at the end, partial answer with gaps in reasoning. Removed SHP index creation part of the recipe as qix can't be used with SQLite dialect – ThomasG77 Apr 30 '21 at 00:34