如何在 Django 中计算 Frechet 距离？

Posted 2023-02-24

【中文标题】如何在 Django 中计算 Frechet 距离？

【英文标题】：How to calculate Frechet Distance in Django?

【发布时间】：2019-10-21 20:54:44

【问题描述】：

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

我必须完成的任务是根据 LineStrings 到给定几何的 Frechet 距离来过滤它们。 Frechet 距离应该使用 PostGIS 提供的ST_FrechetDistance 函数计算。

在另一个基于 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 的项目中 LineStrings 位于 EPSG:4326 中，而在 Django 中它们最初位于 EPSG:3857 中。这是我想出的：

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

它不起作用，因为frechet_annotation QuerySet 抛出异常：

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。当我在没有bytea 演员的情况下手动运行查询时，SQL 工作。所以问题是如何避免这种转换为bytea？

【问题讨论】：

我会使用 shell 尝试逐步分解问题：首先仅使用 Func(F('geom'), Value(4326), function='ST_Transform') 进行注释以查看输出，然后仅使用 Value(wkt) 进行注释并检查每次结果是否是你所期望的。不知道为什么它需要向 bytea 投射一些东西

【参考方案1】：

在您的 SQLAlchemy 示例中，您正在做一些在 GeoDjango 中没有做的事情，即将WKT 字符串转换为Geometry。 这里发生的事情本质上是您尝试使用 PostGIS 函数，但您传递的不是 Geometry，而是一个字符串。

修复第一个问题后我们会偶然发现的另一个问题是以下异常：

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

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

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

这有点令人费解，但是如果我们查看GeoFunc 的代码，我们会看到该类继承了一个名为：GeoFuncMixin 的mixin，它具有geom_param_pos = (0,) 属性并指定了函数参数的位置将是几何形状。 （是的，框架很有趣：P）

我们的函数将输出FloatField。

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

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()

现在我们可以在查询中使用这个函数来计算 ST_FrechetDistance。我们还需要解决将几何传递给函数的原始问题，而不仅仅是 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   

【讨论】：

感谢GeoFunc的用法说明！我在您的答案中编辑了代码 - 调用 reverse() 方法后，srid 数据丢失的错误。

@SS_Rebelious 不错不错！

@SS_Rebelious 实际上这是一个非常好的用例，我已将其添加到链接的 QA（顺便说一句，感谢您在问题中提及）。