如何计算Django的Frechet距离？

Question

这基本上是关于在Django代码中运行自定义PostGIS函数的问题。这个网站上有很多相关的答案，最接近我的例子是this one。建议使用Func()甚至GeoFunc()类，但没有地理空间函数的例子。后者（'GeoFunc'）甚至不适合我抛出st_geofunc does not exist异常（Django 2.1.5）。在

我要完成的任务是根据给定几何体的Frechet距离过滤LineStrings。Frechet距离应该使用PostGIS提供的^{}函数来计算。在

在另一个基于SQLAlchemy的项目中，我使用以下函数完成了完全相同的任务（它正在工作）：

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

正如我所说，上面的函数正在工作，我想在Django中重新实现它。我不得不添加额外的几何图形的SRS转换，因为在基于SQLite的项目中，linestring是EPSG:4326和他们在Django最初是3857。我想到的是：

^{pr2}$

它不起作用，因为frechet_annotation查询集抛出一个异常：

django.db.utils.ProgrammingError: cannot cast type double precision to bytea
LINE 1: ...548 55.717805109,36.825235998 55.717761246)', 0.1)::bytea AS...
                                                             ^

似乎我错误地定义了“ST_FrechetDistance”计算。我怎么修？在

---

更新

签出了Django编写的SQL。总体上是正确的，但是尝试将FrecheDistance的结果转换为bytea会破坏ST_FrechetDistance(...)::bytea。当我在没有byteacast的情况下手动运行查询时，SQL可以工作。所以问题是如何避免这种转换到bytea？在

Answer 1

在SQLAlchemy示例中，您正在做一些在GeoDjango one中没有做的事情，即将WKT字符串转换为Geometry。
这里发生的基本情况是您试图使用一个PostGIS函数，但是您传递给它的不是一个几何体，而是一个字符串。在

在解决第一个问题之后，我们会偶然遇到的另一个问题是以下例外：

django.core.exceptions.FieldError: Cannot resolve expression type, unknown output_field

这就是为什么我们需要创建一个基于GeoFunc的自定义数据库函数。不过，这也带来了一些问题，我们需要考虑以下几点：

• 我们的DB函数将接收2个几何图形作为参数。在

  这有点复杂，但是如果我们看一下^{}的代码，我们会看到这个类继承了一个名为：^{}的mixin，它具有geom_param_pos = (0,)属性，并指定了将成为几何体的函数参数的位置。（是啊，框架很有趣：P）

• 我们的函数将输出FloatField。在

因此，我们的自定义DB函数应该如下所示：

^{pr2}$

现在我们可以在查询中使用此函数来计算^{}。
我们还需要解决将几何图形传递给函数的原始问题，而不仅仅是WKT字符串：

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