8

I'm working on adding some functionality to a QGIS plugin. It currently uses a QgsSpatialIndex to efficiently compute whether the features of a vector layer intersect (potentially millions) of QgsRectangles representing tile bounds. This is reasonably straightforward:

spatialIndex = QgsSpatialIndex()
canvas = self.iface.mapCanvas()
geomTransform = QgsCoordinateTransform(canvas.mapRenderer().destinationCrs(),\
                QgsCoordinateReferenceSystem('EPSG:4326'))
# layer is a QgsVectorLayer, from a combobox selection
vals = {feature.id(): feature for (Feature) in layer.getFeatures()}.values()
for v in vals:
    v.geometry().transform(geomTransform)
map(spatialIndex.insertFeature, vals)

Then later I can check whether a given QgsRectangle intersects features participating in this index:

# includes pseudo-code
if spatialIndex.intersects(rectangle):
    # Feature intersects the rectangle
    # Both partial and complete overlaps are included
    if complete_overlap:
        render_one_way()
    elif partial_overlap:
        render_another_way()
else:
    # No features intersect the rectangle
    dont_render()

This is very fast. However, I need to be able to distinguish between QgsRectangle objects that intersect features in the spatial index partially, and those that overlap completely. The QgsSpatialIndex.intersects method only tells me if there is any kind of intersection.

How can I efficiently determine whether it is a partial or complete overlap?

  • Out of curiosity, what is your aim? Why are you interested in partially overlapping bounding boxes? – user30184 Feb 8 '16 at 23:08
  • 1
    Conditional rendering of tiles. I am producing a cache of land tiles that are transparent over land, and opaque over the sea. These tiles' alpha channels then get inserted during the rendering of other tiles (my own tile server), as a fast way of masking land for oceanic phenomena. The best way to do that is to have oceanic tiles as 0 byte files, and full-land tiles as 1 byte files, so then I can just use os.stat to work out if a mask should be fully opaque or fully transparent without opening the tile image. So the only tiles I really need to render are on the coast: a partial overlap. – alphabetasoup Feb 8 '16 at 23:20
  • The first-pass intersection is already useful generally (github.com/nextgis/QTiles/pull/59), but I want this additional step to make rendering even faster for my own purposes. I need up to zoom level 18, so the only way to make this tractable with a static tile cache is to dramatically reduce the number of tiles that need to be rendered in QGIS. – alphabetasoup Feb 8 '16 at 23:21
  • That makes sense. I can't say right away how to do it without materializing the spatial index by writing the bounding boxes as polygons into a separate table. – user30184 Feb 9 '16 at 8:48
  • Can you get MinX, MinY, MaxX and MaxY out of the object of the spatial index? if you can you could make your own query: If both minx and maxx of the feature are between the minx and maxx of the index object, and the same for the y coordinates, the feature is totally inside. – user30184 Feb 9 '16 at 9:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.