> On 23 Dec 2022, at 11:33 pm, Thierry Florac <tflo...@gmail.com> wrote: > > I'm not sure to have time to test before next year! :/ > Can you just specify the "format" you use to store geometries? > Maybe it's documented in your package but I didn't have any time to check for > it… >
Hi Thierry, The index requires a shapely (https://shapely.readthedocs.io/) geometry from which it calculates a bounding box The - having the full geometry allows for it to do an exact intersection. Any of the https://shapely.readthedocs.io/en/stable/geometry.html#geometry-types should work. I’ve got a very simple demo (below) which takes in data from GADM (https://gadm.org/download_country.html) and there is also some more code in tests for some simplified Belgium geometries. I’m going to write up some more docs to go along with this in the repo but hopefully this gives you enough to get going after Christmas and the new year. If this is getting a bit into the weeds for others on a more general pylons list, say the word and Thierry you can ping me directly. Peter W. from BTrees.IOBTree import BTree from hypatia.spatial import SpatialIndex from persistent import Persistent from shapely.geometry.base import BaseGeometry from hypatia.catalog import Catalog @dataclass(order=True, eq=True) class Region(Persistent): fid: int name: str geometry: BaseGeometry = field(compare=False) @dataclass class Regions(Persistent): catalog: Catalog = field(default_factory=Catalog) container: BTree = field(default_factory=BTree) def __post_init__(self): self.catalog["geometry"] = SpatialIndex("geometry") def add(self, fid, name, geometry): region = Region(fid, name, geometry) self.container[fid] = region self.catalog.index_doc(fid, region) And load data with something like (this pulls in 350k + geometries from the whole world download): from shapely import shape import fiona regions = Regions() for v in filter(None, fiona.open('gadm_410.gpkg')): regions.add(int(v['id']), v['properties']['NAME_1'], shape(v['geometry'])) >From there you can search against any Shapely geometry something like: paris = shapely.geometry.Point(2.349014, 48.864716) near_paris = paris.buffer(0.1) list(regions.catalog['geometry'].intersects(near_paris).execute(resolver=regions.container.get).all()) Which returns 81 areas like this and takes about 1.5ms on my laptop against the 356508 geometries imported above: [Region(fid=76354, name='Île-de-France', geometry=<MULTIPOLYGON (((2.275 48.741, 2.284 48.748, 2.276 48.757, 2.288 48.761, 2.2...>), Region(fid=76355, name='Île-de-France', geometry=<MULTIPOLYGON (((2.319 48.788, 2.315 48.789, 2.31 48.785, 2.303 48.785, 2.29...>), Region(fid=76356, name='Île-de-France', geometry=<MULTIPOLYGON (((2.32 48.771, 2.312 48.772, 2.308 48.779, 2.31 48.785, 2.315...>), Region(fid=76358, name='Île-de-France', geometry=<MULTIPOLYGON (((2.303 48.811, 2.304 48.805, 2.292 48.796, 2.272 48.794, 2.2...>), Region(fid=76359, name='Île-de-France', geometry=<MULTIPOLYGON (((2.229 48.774, 2.227 48.776, 2.227 48.782, 2.235 48.786, 2.2...>), ... I will include a method to quickly do things like within 5km etc, maybe something like: regions.catalog['geometry'].near(paris, 5) -- You received this message because you are subscribed to the Google Groups "pylons-discuss" group. To unsubscribe from this group and stop receiving emails from it, send an email to pylons-discuss+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/pylons-discuss/69582D9C-BF72-4A71-82B7-44AC7B25B0F5%40thirdfloor.com.au.
