This is basically a question about running custom PostGIS functions inside the Django code. There is a number of related answers on this site, most close to my case is this one. It is suggested to use Func() or even GeoFunc() classes but there is no example for geospatial functions there. The latter ('GeoFunc') didn't even work for me throwing st_geofunc does not exist exception (Django 2.1.5).
The task that I have to complete is to filter LineStrings based on their Frechet Distance to the given geometry. Frechet Distance is supposed to be calculated using ST_FrechetDistance function provided by PostGIS.
In another project based on SQLAlchemy I complete the exact same task with the following function (it is working):
from geoalchemy2 import Geography, Geometry
from sqlalchemy import func, cast
def get_matched_segments(wkt: str, freche_threshold: float = 0.002):
matched_segments = db_session.query(RoadElement).filter(
func.ST_Dwithin(
RoadElement.geom,
cast(wkt, Geography),
10
)
).filter(
(func.ST_FrechetDistance(
cast(RoadElement.geom, Geometry),
cast(wkt, Geometry),
0.1
) < freche_threshold) |
# Frechet Distance is sensitive to geometry direction
(func.ST_FrechetDistance(
cast(RoadElement.geom, Geometry),
func.ST_Reverse(cast(wkt, Geometry)),
0.1
) < freche_threshold)
)
return matched_segments
As I said, the function above is working and I wanted to re-implement it in Django. I had to add additional SRS transformation of geometries because in SQLite-based projects LineStrings were in EPSG:4326 and in Django they are in EPSG:3857 initially. Here is what I came up with:
from django.db.models import Func, Value, Q, QuerySet, F
from django.contrib.gis.geos import GEOSGeometry
class HighwayOnlyMotor(models.Model):
geom = LineStringField(srid=3857)
def get_matched_segments(wkt: str, freche_threshold: float = 0.002) -> QuerySet:
linestring = GEOSGeometry(wkt, srid=4326)
transform_ls = linestring.transform(3857, clone=True)
linestring.reverse()
frechet_annotation = HighwayOnlyMotor.objects.filter(
geom__dwithin=(transform_ls, D(m=20))
).annotate(
fre_forward=Func(
Func(F('geom'), Value(4326), function='ST_Transform'),
Value(wkt),
Value(0.1),
function='ST_FrechetDistance'
),
fre_backward=Func(
Func(F('geom'), Value(4326), function='ST_Transform'),
Value(linestring.wkt),
Value(0.1),
function='ST_FrechetDistance'
)
)
matched_segments = frechet_annotation.filter(
Q(fre_forward__lte=freche_threshold) |
Q(fre_backward__lte=freche_threshold)
)
return matched_segments
It doesn't work, as the frechet_annotation QuerySet throws an exception:
django.db.utils.ProgrammingError: cannot cast type double precision to bytea
LINE 1: ...548 55.717805109,36.825235998 55.717761246)', 0.1)::bytea AS...
^
Seems that I incorrectly defined 'ST_FrechetDistance' calculation. How do I fix it?
UPDATE
Checked out the SQL that Django composed. It is overall correct but attempts to cast the result of FrecheDistance to bytea spoils it ST_FrechetDistance(...)::bytea. When I manually run the query without bytea cast, the SQL works. So the question is how to avoid this cast to bytea?
In your SQLAlchemy example, you are doing something that you didn't do in the GeoDjango one and that is to cast the WKT string to Geometry.
What happens here essentially is that you are trying to use a PostGIS function but instead of a Geometry, you are passing it a string.
Another problem that we would stumble upon after fixing the first one would be the following exception:
django.core.exceptions.FieldError: Cannot resolve expression type, unknown output_field
and that is why we need to create a custom database function based on GeoFunc. That poses some problems of its own though and we will need to consider the following:
Our DB Function will receive 2 Geometries as arguments.
That is a bit convoluted, but if we look at the code of GeoFunc we will see that the class inherits a mixin called: GeoFuncMixin which has the attribute geom_param_pos = (0,) and specifies the positions of the function arguments that will be geometries. (Yeaahhh frameworks are fun :P)
FloatField. Therefore our custom DB Function should look like this:
from django.contrib.gis.db.models.functions import GeoFunc
from django.db.models.fields import FloatField
class FrechetDistance(GeoFunc):
function='ST_FrechetDistance'
geom_param_pos = (0, 1,)
output_field = FloatField()
Now we can use this function in our query to calculate the ST_FrechetDistance.
We will also need to address the original issue of passing geometries to the function and not just WKT strings:
def get_matched_segments(wkt: str, freche_threshold: float = 0.002) -> QuerySet:
forward_linestring = GEOSGeometry(wkt, srid=4326)
backward_linestring = GEOSGeometry(wkt, srid=4326)
backward_linestring.reverse()
backward_linestring.srid = 4326 # On Django 2.1.5 `srid` is lost after `reverse()`
transform_ls = linestring.transform(3857, clone=True)
frechet_annotation = HighwayOnlyMotor.objects.filter(
geom__dwithin=(transform_ls, D(m=20))
).annotate(
fre_forward=FrechetDistance(
Func(F('geom'), Value(4326), function='ST_Transform'),
Value(forward_linestring),
Value(0.1)
),
fre_backward=FrechetDistance(
Func(F('geom'), Value(4326), function='ST_Transform'),
Value(backward_linestring),
Value(0.1)
)
)
matched_segments = frechet_annotation.filter(
Q(fre_forward__lte=freche_threshold) |
Q(fre_backward__lte=freche_threshold)
)
return matched_segments