0

I want to add features to a gpkg using python, but before I add them, I want to test if they intersect with the gpkg_layer. For performance, I want to first filter the gpkg_layer using the bounding box. The filter using MbrIntersects() returns all the features in the gpkg_layer, whereas the filter using ST_Intersects returns the correct amount. (At start, the code didn't convert the geometry to bounding box using Envelope, but the behaviour was the same.)

Here is the code with MbrIntersects()

    sql = """\
        SELECT {cols}
        FROM btk_evaluation
        WHERE MbrIntersects(Envelope(GeomFromText('{geometry}', 2169)),geom);
    """.format(cols=cols, geometry=geometry)

Here is the code with ST_Intersects()

    sql = """\
        SELECT {cols}
        FROM btk_evaluation
        WHERE ST_Intersects(Envelope(GeomFromText('{geometry}', 2169)),geom);
    """.format(cols=cols, geometry=geometry)

As you can see, it's exactly the same code, except for the functions used. What am I missing?

2
  • Perhaps the bounding boxes of some features really intersect but the features themselves don't. Have you checked that?
    – user30184
    Feb 13 at 19:11
  • They are at completely different places in a large area.
    – ThiPa
    Feb 13 at 20:27

1 Answer 1

1

There is indeed something wrong with MBRIntersects. It always seems to return -1. ST_EnvIntersects, which is an alias for MBRIntersects, does seem to work though.

Test script:

import geopandas as gpd
import shapely

p1 = shapely.Polygon([[0, 20], [0, 30], [10, 30], [0, 20]])
p2 = shapely.Polygon([[0, 0], [0, 10], [10, 10], [0, 0]])
p3 = shapely.Polygon([[5, 0], [10, 1], [10, 0], [5, 0]])
gdf = gpd.GeoDataFrame(data={"desc": ["p1", "p2", "p3"]}, geometry=[p1, p2, p3])
path = "data.gpkg"
gdf.to_file(path)

geometry = p3.wkt
sql = f"""
    SELECT desc
        ,ST_EnvIntersects(GeomFromText('{geometry}'), d.geom) AS env_intersects 
        ,MBRIntersects(GeomFromText('{geometry}'), d.geom) AS mbr_intersects
        ,ST_Intersects(GeomFromText('{geometry}'), d.geom) AS st_intersects
    FROM data d;
"""

print(gpd.read_file(path, sql=sql, engine="pyogrio"))

Output:

  desc  env_intersects  mbr_intersects  st_intersects
0   p1               0              -1              0
1   p2               1              -1              0
2   p3               1              -1              1

Note though that performance-wise these all aren't great solutions, as they won't use the spatial index of the gpkg. There are many options, but a very easy one that will offer good performance is to use the mask parameter of geopandas.read_file. Under the hood, it will first use the spatial index to find geometries with intersecting MBRs, and then will check for those if the geometries really intersect, which is what you are looking for I think.

import geopandas as gpd
import shapely

p1 = shapely.Polygon([[0, 20], [0, 30], [10, 30], [0, 20]])
p2 = shapely.Polygon([[0, 0], [0, 10], [10, 10], [0, 0]])
p3 = shapely.Polygon([[5, 0], [10, 1], [10, 0], [5, 0]])
gdf = gpd.GeoDataFrame(data={"desc": ["p1", "p2", "p3"]}, geometry=[p1, p2, p3])
path = "data.gpkg"
gdf.to_file(path)

print(gpd.read_file(path, mask=p3))

Output:

  desc                                           geometry
0   p3  POLYGON ((5.00000 0.00000, 10.00000 1.00000, 1...
2
  • Thanks, this works indeed. I will look into 'mask', but I was actually trying to avoid gpd (i don't remember the reason tho). In spatialite, performance can also be increased with the index (gaia-gis.it/gaia-sins/spatialite-tutorial-2.3.1.html#query2). Wasn't yet brave enough to look into it.
    – ThiPa
    Feb 13 at 22:22
  • 1
    geopandas.read_file with the sql parameter will also use spatialite under the hood, so it all boils down to the same thing. The page you link to uses the table and column names for the spatial index as can be expected for "spatialite" files. They are very similar to "Geopackage", but there are some differences. The answer to the following stackoverflow post shows a sample for a Geopackage file.
    – Pieter
    Feb 13 at 22:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.