Центр и радиус дуги в географии
В прошлом посте я писал, что стандартный способ OGC задания дуг по трем точкам представляется мне неестетственным и непрактичным. Вместо 6 координат дугу окружности можно задать 5-ю величинами: координаты центра, радиус, начальный и конечный угол, откладываемые в положительном направлении. В этом посте мы посмотрим, как определить центр окружности, которой принадлежит задааемая 3-мя точками дуга. Радиус определяется элементарно, как расстояние от центра до любой из этих точек. Мы проделаем это упражнение для географии, в геометрии его можно выполнить схожим образом. Более того, благодаря простоте евклидовой метрики задачу нахождения центра в геометрии можно вообще решить практически аналитическим путем.
Выберем некоторые три точки на поверхности Земли и очертим через них дугу окружности.
declare @Astana geography = geography::Point(51.18, 71.40, 4326)
declare @Moscow geography = geography::Point(55.76, 37.62, 4326)
declare @Minsk geography = geography::Point(53.93, 27.63, 4326)
declare @duga geography = geography::STGeomFromText('CIRCULARSTRING(' + Format(@Astana.Long, '#.##') + ' ' + Format(@Astana.Lat, '#.##') + ', ' +
Format(@Moscow.Long, '#.##') + ' ' + Format(@Moscow.Lat, '#.##') + ', ' +
Format(@Minsk.Long, '#.##') + ' ' + Format(@Minsk.Lat , '#.##') + ')', 4326)
select @Astana.STBuffer(25000) union all
select @Moscow.STBuffer(25000) union all
select @Minsk.STBuffer(25000) union all
select @duga.STCurveToLine()
Рис.1
SSMS в CTP3 все еще не умеет рисовать CIRCULARSTRING (см. Рис.4 предыдущего поста), поэтому чтобы нарисовать дугу, ее надо предварительно превратить в ломаную с помощью метода STCurveToLine(). Чтобы точки были заметны, их нужно отобразить маленькими кружочками. Для рисования честного круга появился метод BufferWithCurves(), однако что проку к нему сейчас прибегать, если CURVEPOLYGON она тоже в СТР3 отобразить не сможет. Поэтому аппроксимируем кружочки многоугольничками при помощи старого доброго метода STBuffer(). Все равно эти кружки используются сугубо для визуализации, ни в каких геометрических операциях не участвуют, и на точности вычислений эта аппроксимация не отразится.
Нарисовали дугу, теперь соединяем точки дуги попарно хордами. Для этого выдираем из бинарного представления точек координаты и формируем из них LINESTRINGи о двух точках. Подробнее структура WKB - well-known binary - бинарное представление геопространственных данных разбиралась в посте Превращение последовательности точек в геометрическую фигуру. Хорды выводятся в дополнение к Рис.1:
declare @Astana_Moscow geography = geography::STLineFromWKB(0x01 + 0x02000000 + 0x02000000 + substring(@duga.STAsBinary(), 1 + 4 + 4 + 1, 8 + 8 + 8 + 8), 4326)
declare @Moscow_Minsk geography = geography::STLineFromWKB(0x01 + 0x02000000 + 0x02000000 + substring(@duga.STAsBinary(), 1 + 4 + 4 + 8 + 8 + 1, 8 + 8 + 8 + 8), 4326)
select @Astana.STBuffer(25000) union all
select @Moscow.STBuffer(25000) union all
select @Minsk.STBuffer(25000) union all
select @duga.STCurveToLine() union all
select @Astana_Moscow union all select @Moscow_Minsk
Рис.2
Пересечение перпендикуляров к серединам хорд даст искомый центр окружности. Вспоминаем начальный курс геометрии - как при помощи циркуля и линейки провести перпендикуляр через середину отрезка. Нарисуем окружности радиусом больше половины (я возьму 0.8) длины отрезка с центрами в его вершинах. Окружности будем создавать при помощи метода BufferWithCurves(), а не STBuffer(), чтобы это получились настоящие окружности, а не полигональное приближение, т.к. нас будет интересовать их пересечение, и точность здесь важна.
declare @r1 float = @Astana.STDistance(@Moscow) * 0.8
declare @Astana_circ geography = @Astana.BufferWithCurves(@r1).RingN(1)
declare @Moscow_circ1 geography = @Moscow.BufferWithCurves(@r1).RingN(1)
declare @r2 float = @Minsk.STDistance(@Moscow) * 0.8
declare @Minsk_circ geography = @Minsk.BufferWithCurves(@r2).RingN(1)
declare @Moscow_circ2 geography = @Moscow.BufferWithCurves(@r2).RingN(1)
select @Astana.STBuffer(25000) union all
select @Moscow.STBuffer(25000) union all
select @Minsk.STBuffer(25000) union all
select @duga.STCurveToLine() union all
select @Astana_Moscow union all select @Moscow_Minsk union all
select @Astana_circ.STCurveToLine() union all select @Moscow_circ1.STCurveToLine() union all
select @Minsk_circ.STCurveToLine() union all select @Moscow_circ2.STCurveToLine()
Рис.3
Если внимательно присмотреться к SSMS, видно, что окружность вокруг Москвы, парную с Минской, она все-таки рисует, но таким похабным цветом, что разглядеть ее практически невозможно. Цвета фигур на закладке Spatial results выбираются автоматически, и повлиять на них нельзя. Во всяком случае, я не знаю, как это сделать.
Найдем попарное пересечения окружностей и соединим точки пересечения отрезками. Результатом пересечения двух окружностей является пара точек - MULTIPOINT. Например,
declare @x_Astana_Moscow geography = @Astana_circ.STIntersection(@Moscow_circ1)
select @x_Astana_Moscow.STAsText()
---------------------------------------------------------------------------------------------------------------------------------------------
MULTIPOINT ((51.475251686921311 42.088886593299925), (62.988335441188283 66.977227021261527))
Из этого требуется изготовить отрезок LINESTRING(51.475251686921311 42.088886593299925, 62.988335441188283 66.977227021261527). Сделать это можно через преобразование бинарного формата, как уже проделывалось ранее, например, на Рис.2:
select geography::STGeomFromText('LINESTRING(51.475251686921311 42.088886593299925, 62.988335441188283 66.977227021261527)', 4326).STAsBinary()
select @x_Astana_Moscow.STAsBinary()
---------------------------------------------------------------------------------------------------------------------------------------------
0x01 02000000 02000000 12591A0CD5BC4940 55A3C9A2600B4540, D2B096C6817E4F40 7D4534E38ABE5040
0x01 04000000 02000000 01 01000000 12591A0CD5BC4940 55A3C9A2600B4540, 01 01000000 D2B096C6817E4F40 7D4534E38ABE5040
Требуется из нижней строчки получить верхнюю. Согласно стандарту бинарного представления OGC в нижней строчке первым идет байт-признак Little Endian, затем 4 байта, обозначающие, что это тип MULTIPOINT, затем 4 байта, обозначающие кол-во точек в коллекции, и затем сами эти точки как географические объекты в своем бинарном представлении. Первый байт - признак кодировки Little Endian (в общей бинарной последовательности это будет уже 10-й байт), 4 байта - тип объекта POINT, затем пара 8-байтовых координат этой точки. Далее - вторая точка коллекции. Из этого требуется собрать верхнюю строчку: первый байт, как водится, признак кодировки Little Endian, далее 4 байта - признак типа LINESTRING, далее 4 байта - кол-во точек в ломаной, далее пары 8-байтовых координат точек. Они уже имеются в MULTIPOINT строчкой ниже. Чтобы добраться до первой пары координат, надо пропустить 1 + 4 + 4 + 1 + 4 байта и стать на следующий байт (+1). Поскольку точек всего две, вторая пара будет просто последние 16 байт. По-моему, предельно ясно. Если кому-то еще не понятно - читатйте Превращение последовательности точек в геометрическую фигуру.
declare @x_Astana_Moscow geography = @Astana_circ.STIntersection(@Moscow_circ1)
declare @x_Moscow_Minsk geography = @Minsk_circ.STIntersection(@Moscow_circ2)
declare @i_Astana_Moscow geography = geography::STGeomFromWKB(0x01 + 0x02000000 + 0x02000000 + substring(@x_Astana_Moscow.STAsBinary(), 1 + 4 + 4 + 1 + 4 + 1, 8 + 8) + cast(right(@x_Astana_Moscow.STAsBinary(), 8 + 8) as varbinary), 4326)
declare @i_Moscow_Minsk geography = geography::STGeomFromWKB(0x01 + 0x02000000 + 0x02000000 + substring(@x_Moscow_Minsk.STAsBinary(), 1 + 4 + 4 + 1 + 4 + 1, 8 + 8) + cast(right(@x_Moscow_Minsk.STAsBinary(), 8 + 8) as varbinary), 4326)
select @Astana.STBuffer(25000) union all
select @Moscow.STBuffer(25000) union all
select @Minsk.STBuffer(25000) union all
select @duga.STCurveToLine() union all
select @Astana_Moscow union all select @Moscow_Minsk union all
select @Astana_circ.STCurveToLine() union all select @Moscow_circ1.STCurveToLine() union all
select @Minsk_circ.STCurveToLine() union all select @Moscow_circ2.STCurveToLine() union all
select @i_Astana_Moscow union all select @i_Moscow_Minsk
Рис.4
Ну вот. Теперь блеклым цветом ей взбрело отрисовать отрезок @i_Astana_Moscow, соединяющий точки пересечения больших окружностей. Надеюсь, что читателям все-таки удастся его разглядеть. Чтобы получить центр окружности, проходящей через точки @Astana, @Moscow, @Minsk, следует продлить отрезки @i_Astana_Moscow и @i_Moscow_Minsk и найти их точки пересечения. В географии это даже проще, чем в геометрии. Прямой, проходящей через точки А и Б, на поверхности а ля сферы считается след от сечения поверхности плоскостью, проходящей через ее центр и А и Б. Чтобы нарисовать проходящий через точки А и Б экватор, достаточно принять во внимание, что на нем будут также находиться их антиподы, т.е. это будет LINESTRING, задаваемая, как А, Б, антипод А, антипод Б, А. Без антиподов эта схема, к сожалению, не работает. Если задать LINESTRING как А, Б, А, SSMS выдает, что фигура not valid при выполнении операции пересечения прямых. Видимо, SQL Server не догадывается, что звено Б-А надо пустить в обход земного шара, а пускает обратно внахлест.
Под антиподом точки А понимается протвоположная ей точка на обратной стороне Земли, там, где протыкает глобус спица, проходящая через А и центр Земли. Координаты антипода А с координатами @x в.д., @у с.ш. будут 180 - @х з.д., @у ю.ш. и наоборот, т.е. долготу @x ∈[-180, 180] надо обратить в - sign(@x) * (180 - abs(@x)), а широту @у ∈[-90, 90]- просто взять с противоположным знаком. На эту тему я даже не поленился написать функцию:
if exists(select 1 from sys.objects where name = 'udfAntipodalPoint' and type = 'FN' and schema_id = SCHEMA_ID('dbo')) drop function dbo.udfAntipodalPoint
go
create function dbo.udfAntipodalPoint(@p geography) returns geography as
begin
declare @antip geography
if @p.InstanceOf('POINT') = 1 set @antip = geography::Point(-@p.Lat, -sign(@p.Long) * (180 - abs(@p.Long)), @p.STSrid)
return @antip
end
Досадно, что несмотря на появление оператора THROW в Denali обработка ошибок в пользовательских функциях осталась на прежнем уровне, то есть никакая. Операторы RAISERROR, THROW, конструкцию BEGIN TRY … BEGIN CATCH в UDF использовать нельзя. По уму функцию udfAntipodalPoint следовало писать как-то так:
create function udfAntipodalPoint(@p geography) returns geography as
begin
if @p.InstanceOf('POINT') = 0 throw 66666, 'Произошла ошибка. Мне на входе требуется точка', 1
...
Но - Msg 443, Level 16, State 14... Invalid use of a side-effecting operator 'THROW' within a function. Складывается впечатление, что люди, которые когда-то реализовали пользовательские функции в SQL Server 2000, все давно поувольнялись и с тех пор этим механизмом никто всерьез не занимался.
Сделаем при помощи функции udfAntipodalPoint прямые из отрезков @i_Astana_Moscow, @i_Moscow_Minsk прямые двумя способами для разнообразия: через текстовое и бинарное представления.
declare @l_Astana_Moscow geography = geography::STGeomFromText('LINESTRING(' + Format(@i_Astana_Moscow.STPointN(1).Long, 'N13') + ' ' +
Format(@i_Astana_Moscow.STPointN(1).Lat , 'N13') + ',' +
Format(@i_Astana_Moscow.STPointN(2).Long, 'N13') + ' ' +
Format(@i_Astana_Moscow.STPointN(2).Lat , 'N13') + ',' +
Format(dbo.udfAntipodalPoint(@i_Astana_Moscow.STPointN(1)).Long, 'N13') + ' ' +
Format(dbo.udfAntipodalPoint(@i_Astana_Moscow.STPointN(1)).Lat , 'N13') + ',' +
Format(dbo.udfAntipodalPoint(@i_Astana_Moscow.STPointN(2)).Long, 'N13') + ' ' +
Format(dbo.udfAntipodalPoint(@i_Astana_Moscow.STPointN(2)).Lat , 'N13') + ',' +
Format(@i_Astana_Moscow.STPointN(1).Long, 'N13') + ' ' +
Format(@i_Astana_Moscow.STPointN(1).Lat , 'N13') + ')',
4326)
--Берем за основу имеющийся отрезок @i_Moscow_Minsk. Увеличиваем кол-во точек в LINESTRING до 5-ти (0x05000000)
--и добавляем координаты антиподных и начальной точек
declare @l_Moscow_Minsk geography = geography::STGeomFromWKB(cast(stuff(@i_Moscow_Minsk.STAsBinary(), 1 + 4 + 1, 4, 0x05000000) as varbinary(max)) +
cast(right(dbo.udfAntipodalPoint(@i_Moscow_Minsk.STPointN(1)).STAsBinary(), 8 + 8) as varbinary(max)) +
cast(right(dbo.udfAntipodalPoint(@i_Moscow_Minsk.STPointN(2)).STAsBinary(), 8 + 8) as varbinary(max)) +
substring(@i_Moscow_Minsk.STAsBinary(), 1 + 4 + 4 + 1, 8 + 8),
4326)
select @duga.STCurveToLine() union all
select @Astana_Moscow union all select @Moscow_Minsk union all
select @Astana_circ.STCurveToLine() union all select @Moscow_circ1.STCurveToLine() union all
select @Minsk_circ.STCurveToLine() union all select @Moscow_circ2.STCurveToLine() union all
select @l_Astana_Moscow union all select @l_Astana_Moscow union all select @l_Moscow_Minsk
Рис.5
Мне пришлось прочертить прямую @l_Astana_Moscow два раза, т.к. SSMS опять вознамерилась первый раз рисовать ее тусклым цветом. Так хоть что-то, я надеюсь, заметно.
Результатом пересечения двух прямых, в том числе параллельных, является пара точек.. Из них по смыслу выбираем ту, которая ближе. Она и будет центром проходящей через Астану, Москву и Минск окружности, а расстояние до любой из трех точек - ее радиусом.
declare @кандидатывцентры geography = @l_Astana_Moscow.STIntersection(@l_Moscow_Minsk)
declare @центр geography = iif (@Moscow.STDistance(@кандидатывцентры.STPointN(1)) < @Moscow.STDistance(@кандидатывцентры.STPointN(2)), @кандидатывцентры.STPointN(1), @кандидатывцентры.STPointN(2))
declare @радиус float = @Moscow.STDistance(@центр)
select @duga.STBuffer(100000) union all
select @Astana_circ.STCurveToLine() union all select @Moscow_circ1.STCurveToLine() union all
select @Minsk_circ.STCurveToLine() union all select @Moscow_circ2.STCurveToLine() union all
select @l_Astana_Moscow union all select @l_Astana_Moscow union all select @l_Moscow_Minsk union all
select @центр.STBuffer(200000) union all select @центр.BufferWithCurves(@радиус).RingN(1).STCurveToLine()
union all select @центр.BufferWithCurves(@радиус).RingN(1).STCurveToLine()
Рис.6
Собственно, все. Мы определили центр и радиус окружности, которой принадлежит заданная дуга из трех точек. Самым сложным в этом упражнении оказалось побороть цвета SSMS. Так, искомую окружность я прорисовываю два раза не только на радостях, что ее нашел, а потому что по закону подлости SSMS выбрала самый незаметный цвет, чтобы отобразить венец трудов. При помощи метода STBuffer() пожирнее были выделены исходная дуга и найденный центр. Обратите внимание, что SQL Server Denali поддерживает операции на полном глобусе в отличие от своих предшественников 2008/R2, которые умели работать только с полушарием. При сравнении расстояния до дуги от двух точек пересечения прямых, чтобы определить, кто из них в действительности получается ее центром, использовалась функция iif(), которая, к слову сказать, тоже является маленьким приятным новшеством T-SQL в Denali.
Обработка исключительных ситуаций, например, когда дуга оказывается отрезком прямой, предоставляется читателям в качестве самостоятельного упражнения.
Алексей Шуленин