Michael,

Thanks for the quick and awesome answer.  I was spending my time trying to
figure out how to store an association between the ID in the index and geom
and never thought about storing the actual geom in the index.  I had tried
to use Objects=True with RTree, but didn't get it to work.

I created a bunch of test scripts to determine which strategy was fastest.
I posted the code as gist's for anyone that wants to look at them.  There
is definitely room for cleaning up the code.  I was also going for
readability over compactness.  Because there really aren't a lot of
examples for fiona/shapely/rtree out there, I am interested in seeing how
other people might code it.  Functional style would be interesting too.

My test data set was 7142 points in 87 polygons with mildly overlapping
bboxes.

Here are my examples in order of processing speed (fastest to slowest):

RTree index and shapely optimizations - 5.01 seconds
https://gist.github.com/fawcett/5251327

Added shapely speedups and prepared geometries - 7.07 seconds
https://gist.github.com/fawcett/5251267

Basic point in poly, no optimization - 10.14 seconds
https://gist.github.com/fawcett/5251238

Rtree index, storing the geoms in the index as objects - 66.10
(only doing true intersects on points:polys where a point hits multiple
RTree leaves)
https://gist.github.com/fawcett/5251373

Rtree index, storing the geoms in the index as objects - 68.85
(doing intersects on all point RTree leaf hits)
https://gist.github.com/fawcett/5251407

Obviously something went very wrong on the last two.  They use the RTree,
but are an order of magnitude slower than the most basic RTree example.  I
am definitely curious if my code is bad or if it has to do with the way
that geometries and properties are store within the RTree.


And yes, you do get extra points for using fiona!  I am using it in my
scripts.  (We should really push this example to gis.se for more visibility
for fiona and shapely.)

David.

On Tue, Mar 26, 2013 at 4:52 PM, Michael Weisman <mweis...@gmail.com> wrote:

> Hi David,
>
> No idea if this is best practice or if hobu and Sean will be shaking their
> heads in disappointment in a second, but RTree allows one to store an
> object (anything pickleable I think) in the index. You can store a
> dictionary with attributes and the geometry from a fiona feature like this:
> geom = geometry.shape(feature['geometry'])
> idx.insert(int(feature['id']), geom.bounds, obj={'properties':
> feature['properties'], 'geom':geom})
>
> Now I can iterate over the results of the RTree intersection query with
> shapely:
> for riding in idx.intersection((point.y, point.x), objects="raw"):
>         if riding['geom'].contains(point):
>             return riding
>
> There may be a better way to do this, but in my experience the above works
> quite well.
>
> Do I get bonus points for also using Fiona in my solution?
>
> Thanks,
>
> Michael
>
> On 2013-03-26, at 2:17 PM, David Fawcett <david.fawc...@gmail.com> wrote:
>
> > I am thinking through the logic for efficiently overlaying ~10,000
> points on ~100 polys.  What I want in the end is a list of dicts holding
> the string ID for the point and string ID for the intersecting Poly.   I
> have an example with prepared geometries, but now I am working up an RTree
> example to benchmark.
> >
> > Here is the logic:
> >
> > - I use RTree to create an index on my polys
> > - at the same time, I create a dict to lookup the string ID (stationID)
> for the integer index value
> >
> > - I loop through the points in my collection and test them against the
> index.
> > - if only one intersection is returned, I grab the id and look up my
> countyID.
> >
> > - if my intersection operation on the index returns multiple poly IDs, I
> need to test against the actual geometries.
> >   This is where I am struggling.  I don't see a way to access a specific
> feature from a shapely collection based on a property.  My first idea is to
> store the polygon features in a dictionary after reading them from the
> shapefile.  This would give me access by ID.  At the same time, I am trying
> to work within the Shapely data structures.
> >
> > Does anyone have any suggestions on best practices for integrating RTree
> and Shapely?  Specifically, how you relate geometry collections to RTree
> indexes.
> >
> > I will post a code example after I modify my current example to not use
> external data sources.
> >
> > Thanks,
> >
> > David.
> >
> >
> >
> > _______________________________________________
> > Community mailing list
> > Community@lists.gispython.org
> > http://lists.gispython.org/mailman/listinfo/community
>
> _______________________________________________
> Community mailing list
> Community@lists.gispython.org
> http://lists.gispython.org/mailman/listinfo/community
>
_______________________________________________
Community mailing list
Community@lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community

Reply via email to