Геопространственные агрегаты в Denali

Article
09/22/2011

Новой возможностью SQL Server 11 в части геопространственных расширений являются агрегатные функции CollectionAggregate, ConvexHullAggregate, EnvelopeAggregate и UnionAggregate для типов Geometry и Geography. Они работают так же, как и обычные SUM(), COUNT(), AVG() и т.д., но применительно к геопространственным типам данных. UnionAggregate строит суммарное объединение фигур в колонке, EnvelopeAggregate создает общий заключающий их прямоугольник (круг), ConvexHullAggregate натягивает на всех минимальную выпуклую оболочку, а CollectionAggregate добавляет их в результирующую GeometryCollection. Например, если стандартный метод STUnion() строит объединение двух объектов - @g1.STUnion(@g2), то функция UnionAggregate() пробегает вдоль колонки типа geometry или geography и строит объединение всех находящихся в ее ячейках геопространственных объектов - select UnionAggregate(g) from tbl. До появления геопространственных агрегатных функций подобные задачи приходилось решать традиционно через курсор, либо созданием CLRного агрегата, поэтому общественность давно высказывала пожелание иметь встроенную функциональность. Для SQL Server 2008 / R2 такие агрегаты на CLR были написаны в рамках открытого проекта SQL Server 2008 Spatial tools на Codeplex. Хочется надеяться, что и другие полезные вошедшие в него вещи, включая интерполяцию, аффинные преобразования и др., также станут штатными возможностями SQL Server.

Рассмотрим использование геопространственного агрегата на примере. Пусть имеется некоторое количество базовых станций c известными координатами и радиусами действия. Требуется определить, покрывает ли созданная ими сотовая сеть заданную территорию. В качестве территории я возьму построенный в посте Ориентация полигона на поверхности (см. тамошний Скрипт 1) внутримкадный многоугольник.

Чтобы не обижать никого из сотовых операторов, точки, обозначающие базовые станции, я набросаю случайным образом в описанный вокруг МКАД прямоугольник.

--Превращаем точки из таблицы #mkad в многоугольник.

--См. http://blogs.technet.com/b/isv\_team/archive/2011/09/11/3452373.aspx

declare @s nvarchar(max) = (select Format(p.Lat, 'N6') + ' ' + Format(p.Long, 'N6') + ' ' from #mkad order by km desc for xml path(''))

set @s += (select top 1 Format(p.Lat, 'N6') + ' ' + Format(p.Long, 'N6') from #mkad order by km desc)

declare @m geography = geography::GeomFromGml('<Polygon xmlns="http://www.opengis.net/gml">

</LinearRing>

</exterior>

</Polygon>', 4326)

--select @m

Скрипт 1

Теперь вокруг @m нужно описать минимальный содержащий его прямоугольник, стороны которого ориентированы вдоль координатных осей (axis-aligned bounding box). Это то, что делает в геометрии стандартный метод STEnvelope(). К сожалению, в географии такого метода нет, а упоминавшийся выше агрегат EnvelopeAggregate() в случае географии строит не прямоугольник, а круг типа CURVEPOLYGON с центром в EnvelopeCenter и радиусом EnvelopeAngle, которые мы разбирали в прошлом посте.

select geography::EnvelopeAggregate(@m).ToString()

select @m

union all

select geography::EnvelopeAggregate(@m).STCurveToLine()

---------------------------------------------------------------------------------------------------------------------------------------

CURVEPOLYGON (CIRCULARSTRING (37.834996307272576 55.869466790654435, 37.390710492930147 55.869466790654442, 37.392120687180309 55.620185967974976, 37.833586113022413 55.620185967974976, 37.834996307272576 55.869466790654435))

Скрипт 2

Для полигона, не включающего полюс, можно тупо обойти вершины, найти минимальную и ммаксимальную широту и долготу его вершин и из этих 4-х точек построить прямоугольник.

--Определяем координаты минимального описывающего прямоугольника

declare @широта_min float, @широта_max float, @долгота_min float, @долгота_max float

select @широта_min = min(p.Lat), @широта_max = max(p.Lat), @долгота_min = min(p.Long), @долгота_max = max(p.Long) from #mkad

declare @r geography = geography::STPolyFromText('POLYGON((' + Format(@долгота_max, 'N6') + ' ' + Format(@широта_max, 'N6') + ', ' +

Format(@долгота_min, 'N6') + ' ' + Format(@широта_max, 'N6') + ', ' +

Format(@долгота_min, 'N6') + ' ' + Format(@широта_min, 'N6') + ', ' +

Format(@долгота_max, 'N6') + ' ' + Format(@широта_min, 'N6') + ', ' +

Format(@долгота_max, 'N6') + ' ' + Format(@широта_max, 'N6') + '))', 4326)

select @m union all select @r

Скрипт 3

Как и в случае равномерной плоскости (geometry) сторонами POLYGONа в geography являются отрезки прямых, просто под прямой понимается большая окружность, полученная как результат сечения поверхности Земли плоскостью, проходящей через центр Земли. Для CURVEPOLYGON, сторонами которого выступают дуги, т.е. куски окружностей, нарисованных на поверхности Земли, построить прямоугольный конверт проблематичней, т.к. крайней точкой может оказаться не вершина, а точка где-нибудь на дуге. Хотя, в принципе, каждую дугу можно аппроксимировать ломаной с приемлемой точностью при помощи метода CurveToLineWithTolerance().

Набросаем в таблицу #cells @n точек со случайными координатами в пределах прямоугольника, где каждая точка будет центром окружности случайного радиуса от @радиус_min до @радиус_max. Я не специалист в сотовой связи, поэтому значения минимального и максимального радиуса покрытия базовой станции выбирал от балды.

if OBJECT_ID('tempdb..#cells', 'U') is not null drop table #cells

create table #cells (id int identity primary key, центр_широта float, центр_долгота float, радиус float)

declare @радиус_min float = 3000, @радиус_max float = 5000

declare @n int = 50, @i int = 0

while @i < @n begin

insert #cells (центр_широта, центр_долгота, радиус) values (@широта_min + rand() * (@широта_max - @широта_min),

@долгота_min + rand() * (@долгота_max - @долгота_min),

@радиус_min + rand() * (@радиус_max - @радиус_min))

set @i += 1

end

select * from #cells

Скрипт 4

По-хорошему, это не совсем равномерное распределение, т.к. в отличие от евклидовой плоскости, где метрика по X и по Y постоянна, на Земле, чем выше широта, тем меньше метров содержит каждый градус долготы. Однако для размеров МКАД, я думаю, можно пренебречь этой погрешностью. Добавим в таблицу #cells, собственно, круги:

alter table #cells add cell geography

update #cells set cell = geography::Point(центр_широта, центр_долгота, 4326).STBuffer(радиус)

select * from #cells

Скрипт 5

У, красота какая. Теперь нужно найти объединение всех этих кругов и посмотреть, накрывает ли собой их объединение внутримкадный многоугольник или же в нем остаются какие-либо дырки. Для объединения фигур в колонке cell используем появившуюся в Денали агрегатную функцию UnionAggregate():

--Пример геопространственного агрегата:

declare @coverage geography

select @coverage = geography::UnionAggregate(cell) from #cells

select @coverage union all select @m union all select @r

Скрипт 6

Бордовым цветом показана территория города в черте МКАД, синим – описывающий ее прямоугольник, голубым – результат объединения набросанных в него кругов (Скрипт 5). Визуально можно определить, что голубая фигура накрывает бордовую с брешами. Чтобы сделать это не на глазок, а точно, можно использовать метод STContains(), который, что особенно приятно, в Денали стал поддерживаться для географии. В 2008/R2 он был только в геометрии.

select @coverage.STContains(@m)

Скрипт 7

Мое сугубое пожелание к разработчикам геопространственных расширений SQL Server – добавить возможность задания порядка при расчете агрегата подобно тому, как это было сделано для sum(), count(), next value и др. Мы разбирали возможность задания порядка на оконных агрегатах, в частности, в посте Нарастающий итог в Денали. К сожалению, конструкция over(order by …) не поддерживается для геопространственных агрегатов. Это нелогично, т.к. GEOMETRYCOLLECTION предполагается упорядоченной – см. метод STGeometryN(), позволяющий вытащить N-й член коллекции. Однако при выполнении запроса с агрегатной функцией CollectionAggregate() порядок, в котором SQL Server будет перебирать записи в процессе выполнения select не гарантирован, от чего смысл этой агрегатной функции в известном смысле теряется.

Еще один пример. У меня была мысль применить CollectionAggregate() в посте Превращение последовательности точек в геометрическую фигуру. Можно видеть, что бинарные представления коллекции точек и LINESTRING из них очень схожи:

declare @l varbinary(max) = geometry::STLineFromText('LINESTRING(0 0, 1 1, 1 -1, -1 -1, -1 1)', 0).STAsBinary()

select @l

if OBJECT_ID('tempdb..#points', 'U') is not null drop table #points

create table #points (id int identity primary key, p geometry)

insert #points (p) values (geometry::Point(0, 0, 0)), (geometry::Point(1, 1, 0)), (geometry::Point(1, -1, 0)), (geometry::Point(-1, -1, 0)), (geometry::Point(-1, 1, 0))

declare @x varbinary(max) = (select Geometry::CollectionAggregate(p) from #points).STAsBinary()

select @x

select geometry::Point(0, 0, 0).STAsBinary()

---------------------------------------------------------------------------------------------------------------------------------------

0x01 02000000 05000000 00000000000000000000000000000000 000000000000F03F000000000000F03F 000000000000F03F000000000000F0BF 000000000000F0BF000000000000F0BF 000000000000F0BF000000000000F03F

0x01 07000000 05000000 010100000000000000000000000000000000000000 0101000000000000000000F03F000000000000F03F 0101000000000000000000F03F000000000000F0BF 0101000000000000000000F0BF000000000000F0BF 0101000000000000000000F0BF000000000000F03F

0x010100000000000000000000000000000000000000

Скрипт 8

Вначале идет байт-признак Little Endian, затем 4 байта – геопространственный тип 02000000 = LINESTRING, 07000000 = GEOMETRYCOLLECTION, , затем 4 байта 05000000 = кол-во объектов в коллекции / точек в LINESTRING, затем в LINESTRING идут координаты точек, а в коллекции – сами объекты, т.е. в данном случае – точки, у которых первый байт – это тоже Little Endian, затем 4 байта – геопространственный тип 01000000 = POINT и затем по 8 байт координаты точки.

Можно было бы заменить у коллекции геопространственный тип, обрезать у каждого объекта точки 5-байтный префикс, оставив только координаты, и получить вожделенный LINESTRING. К сожалению, от него пришлось отказаться, ибо, как я уже сказал, никто не обещал, что CollectionAggregate() будет соблюдать заданный порядок.

Менее быстрый, но более наглядный способ состоял бы в том, чтобы подкорректировать GML-представление – см. Продолжаем соединять точки в контур.

--Это есть:

select Geometry::CollectionAggregate(p).AsGml() from #points

--Нужно получить:

select geometry::STLineFromText('LINESTRING(0 0, 1 1, 1 -1, -1 -1, -1 1)', 0).AsGml()

---------------------------------------------------------------------------------------------------------------------------------------

<Point>

</Point>

<Point>

</Point>

<Point>

</Point>

<Point>

</Point>

<Point>

</Point>

</geometryMembers>

</MultiGeometry>

</LineString>

Скрипт 9

Я отказался от него сразу по двум причинам. Во-первых, все та же невозможность гарантировать порядок перечисления записей при вычислении агрегата. Во-вторых, я не знаю, как написать XQuery, который бы из верхнего XML делал нижний, т.к. функция string-join() в SQL Serverном варианте XQuery не поддерживается. В итоге пришлось строить нижний XML вручную, сконкатенировав строку координат при помощи известного хакерского приема через for xml path. Мы его видели в Скрипт 1. Как утверждается, SELECT … FOR XML строит XML (в данном случае plain text строку) в порядке, заданном ORDER BY.

Продолжение следует.

Алексей Шуленин

Геопространственные агрегаты в Denali

Additional resources