Как рассчитать угол, под которым две линии пересекаются в PostGIS?

Я хочу рассчитать угол между двумя линиями, где они пересекаются в PostGIS.

Начальная точка для расчета угла в PostGIS, кажется, ST_Azimuth - но она принимает точки в качестве входных данных. Моей первой мыслью было взять конечные точки пересекающихся линий и выполнить для них азимутальные вычисления. Это недостаточно хорошо, потому что большинство линейных объектов не прямые, и меня интересует угол на пересечении. Итак, я придумал вложенную операцию, которая проходит следующие шаги:

Определите все пересечения между двумя таблицами линейных объектов.
Создайте очень маленький буфер вокруг точки пересечения
Определите точки, где линейные объекты пересекают внешнюю часть буфера (беря первую точку, если их больше одной - меня действительно интересует только, близок ли угол к 0, 90 или 180 градусам)
Рассчитайте ST_Azimuth для этих двух точек.

Полный SQL довольно долго размещать здесь, но я привел его здесь, если вам интересно. (Кстати, есть ли лучший способ, чем перенести все поля, идущие вниз по операторам WITH?)

Результаты выглядят неправильно, поэтому я явно делаю что-то не так:

выходной пример 1 выходной пример 2

РЕДАКТИРОВАТЬ Я переделал расчеты в EPSG: 3785, и результаты немного отличаются, но все еще не правильно:

выход в 3785 # 1 выход в 3785 # 2

Мой вопрос, где недостатки в этом процессе. Я неправильно понимаю, что делает ST_Azimuth? Есть ли проблема с CRS? Что-то еще вообще? Или, может быть, есть гораздо более простой способ сделать это?

postgis geometry mvexel
источник

Каким был оригинальный CRS? Вычисления угла должны выполняться с помощью конформной проекции, а не с непроецированной широтой / долготой (SRID = 4326).

Майк Т

Это был EPSG: изначально было 4326 координат, я включил ST_Translate, чтобы быть на 100% уверенным, что вся обработка будет выполняться в одном и том же CRS. Я попробую конформную проекцию, спасибо.

mvexel

Я переделал расчеты EPSG: 3785, и это действительно имеет значение - я исправлю вопрос, чтобы показать новые результаты - но результат все еще не отражает фактический угол.

mvexel

Ответы:

У меня было прозрение. Это довольно обыденно. Я оставил одну важную информацию для PostGIS, чтобы вычислить правильный угол.

Я рассчитывал угол между двумя точками, пересекающими внешнюю поверхность небольшого буфера. Чтобы вычислить угол пересечения, мне нужно рассчитать оба угла между обеими точками на внешней стороне буфера и точкой пересечения двух линейных объектов и вычесть их.

Я обновил полный SQL , но вот что заметно:

SELECT
    ...
    abs
    (
        round
        (
            degrees
            (
            ST_Azimuth
            (
                points.point2,
                points.intersection
            )
            -
            ST_Azimuth
            (
                points.point1,
                points.intersection
            )
        )::decimal % 180.0
        ,2
    )
)
AS angle
...
FROM
points

mvexel
источник

Я думал об угле буферизованной точки относительно пересечения, но у меня нет времени вдаваться в детали. Другой аспект - угловые единицы. Вам нужно умножить результат в радианах из ST_Azimuth на 180.0 / pi (), чтобы получить результаты в градусах.

Майк Т

Да, спасибо, я использую функцию PostgreSQL Степени () для этого.

mvexel

Кливер. (Я даже не знаю , что было градусы не функционируют до сих пор.) Было бы хорошо , чтобы обернуть всю эту логику в вызове функции, но у меня возникают трудности в концептуализации , как это будет работать, то есть ST_IntersectionAngle(...?

Майк Т

Я был на самом деле удивлен, что это не функция PostGIS. Спасибо за ваш отзыв об этом.

mvexel

Мне недавно пришлось рассчитывать то же самое, но я выбрал более простой и, вероятно, более быстрый подход.

Чтобы найти дополнительные точки для вычисления азимута, я просто проверяю пермириаду длины за пересечением (или после того, в редком случае, когда это происходит в самом начале линии), используя ST_Line_Locate_Point и ST_Line_Interpolate_Point :

abs(degrees( 
  ST_Azimuth (
    intersection, 
    ST_Line_Interpolate_Point(
      line1, 
      abs(ST_Line_Locate_Point(line1, intersection) - 0.0001)
    )
  )
  -
  ST_Azimuth (
    intersection, 
    ST_Line_Interpolate_Point(
      line2, 
      abs(ST_Line_Locate_Point(line2, intersection) - 0.0001)
    )
  )
))

Пермириада была произвольной, и для более последовательных результатов было бы лучше использовать абсолютное смещение. Например, чтобы заранее проверить 20 м, вы должны изменить 0,0001 на 20/ST_Length(line1)и 20/ST_Length(line2)соответственно.

lynxlynxlynx
источник