GeoDjango does not create correct PostGIS index for 3D geometry fields

[Casper van der Wel]

Ref https://code.djangoproject.com/ticket/30267

Hi,

I stumbled again on this (one year old) ticket and wanted to follow up on the suggestion of Tim Graham to ask your opinion about this issue.

When I create a model with a 3D geometryfield, an index USING gist (geom gist_geometry_ops_nd) is created automatically by geodjango in my PostGIS database. However, the spatial lookups 'intersects' and 'bboverlaps' do not use this index, which I can confirm by explaining the queries. Below some code to reproduce this issue.

I am not sure if this is a bug or if I am reporting a non-intended use case for 3D geometryfields. It appears to me as a pretty straightforward use case.

The model:

from django.contrib.gis.db import models
from django.db import connections
from django.db.models.query import QuerySet

class QuerySetExplain(QuerySet):
def explain(self):
cursor = connections[self.db].cursor()
query, params = self.query.sql_with_params()
cursor.execute('EXPLAIN %s' % query, params)
return '\n'.join(r[0] for r in cursor.fetchall())

class PolyModel(models.Model):
objects = QuerySetExplain.as_manager()
geometry = models.PolygonField(srid=28992, dim=3)
After creating and running migrations, I explain the following query:

>>> from django.contrib.gis.geos import Polygon
>>> from my_app.models import PolyModel
>>> bbox = Polygon.from_bbox((0, 0, 1, 1))
>>> qs = PolyModel.objects.filter(geometry__intersects=bbox)
>>> print(qs.explain())

Seq Scan on my_app_polymodel (cost=0.00..2964.00 rows=1 width=241)
Filter: ((geometry && '0103000020E61000000100000005000000000000000000000000000000000000000000000000000000000000000000F03F000000000000F03F000000000000F03F000000000000F03F000000000000000000000000000000000000000000000000'::geometry) AND _st_intersects(geometry, '0103000020E61000000100000005000000000000000000000000000000000000000000000000000000000000000000F03F000000000000F03F000000000000F03F000000000000F03F000000000000000000000000000000000000000000000000'::geometry))
This says that is used a sequential scan. It did not use the index
The solution for me is creating the indices "USING gist (geometry)" with raw SQL.
Looking at PostGIS manuals, it seems you have to use the operator &&& to actually use an ND index. I can confirm that in the sql shell.

Versions: I am using django 1.11.20, PostGRES 9.3 with PostGIS 2.4. I also reproduced the issue in raw SQL using PostGRES 10.

After glancing the django source code, it appears to me that it persists in master. But I did not run any tests.

A quickfix would be always creating 2 indexes on 3D geometry fields. But that does seem a waste of resources. Using &&& instead of && as the 'bboverlaps' lookup would probably a better fix, but I am not such a PostGIS expert that I can oversee the changes necessary in the django ORM.

Thanks in advance for looking at this bug report.