Partager via

Ориентация полигона на поверхности

  • Article

Задача. Найти расстояние от Москвы до мест проведения конференции разработчиков DevCon в 2011 и 1996 г.

 

Чтобы было интересней, под расстоянием от Москвы будем понимать не традиционное расстояние от нулевого километра 55.755831 с.ш., 37.617673 в.д. (это тривиально), а, скажем, кратчайшее расстояние от МКАД, координаты каждого километра которого указаны, например, здесь.

 

use tempdb

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

create table #mkad (km int primary key, p geography)

 

insert #mkad (km, p) values

( 1, geography::Point(55.774558, 37.842762, 4326)),

( 2, geography::Point(55.76522, 37.842789, 4326)),

( 3, geography::Point(55.755723, 37.842627, 4326)),

( 4, geography::Point(55.747399, 37.841828, 4326)),

( 5, geography::Point(55.739103, 37.841217, 4326)),

( 6, geography::Point(55.730482, 37.840175, 4326)),

( 7, geography::Point(55.721939, 37.83916, 4326)),

( 8, geography::Point(55.712203, 37.837121, 4326)),

( 9, geography::Point(55.703048, 37.83262, 4326)),

( 10, geography::Point(55.694287, 37.829512, 4326)),

( 11, geography::Point(55.68529, 37.831353, 4326)),

( 12, geography::Point(55.675945, 37.834605, 4326)),

( 13, geography::Point(55.667752, 37.837597, 4326)),

( 14, geography::Point(55.658667, 37.839348, 4326)),

( 15, geography::Point(55.650053, 37.833842, 4326)),

( 16, geography::Point(55.643713, 37.824787, 4326)),

( 17, geography::Point(55.637347, 37.814564, 4326)),

( 18, geography::Point(55.62913, 37.802473, 4326)),

( 19, geography::Point(55.623758, 37.794235, 4326)),

( 20, geography::Point(55.617713, 37.781928, 4326)),

( 21, geography::Point(55.611755, 37.771139, 4326)),

( 22, geography::Point(55.604956, 37.758725, 4326)),

( 23, geography::Point(55.599677, 37.747945, 4326)),

( 24, geography::Point(55.594143, 37.734785, 4326)),

( 25, geography::Point(55.589234, 37.723062, 4326)),

( 26, geography::Point(55.583983, 37.709425, 4326)),

( 27, geography::Point(55.578834, 37.696256, 4326)),

( 28, geography::Point(55.574019, 37.683167, 4326)),

( 29, geography::Point(55.571999, 37.668911, 4326)),

( 30, geography::Point(55.573093, 37.647765, 4326)),

( 31, geography::Point(55.573928, 37.633419, 4326)),

( 32, geography::Point(55.574732, 37.616719, 4326)),

( 33, geography::Point(55.575816, 37.60107, 4326)),

( 34, geography::Point(55.5778, 37.586536, 4326)),

( 35, geography::Point(55.581271, 37.571938, 4326)),

( 36, geography::Point(55.585143, 37.555732, 4326)),

( 37, geography::Point(55.587509, 37.545132, 4326)),

( 38, geography::Point(55.5922, 37.526366, 4326)),

( 39, geography::Point(55.594728, 37.516108, 4326)),

( 40, geography::Point(55.60249, 37.502274, 4326)),

( 41, geography::Point(55.609685, 37.49391, 4326)),

( 42, geography::Point(55.617424, 37.484846, 4326)),

( 43, geography::Point(55.625801, 37.474668, 4326)),

( 44, geography::Point(55.630207, 37.469925, 4326)),

( 45, geography::Point(55.641041, 37.456864, 4326)),

( 46, geography::Point(55.648794, 37.448195, 4326)),

( 47, geography::Point(55.654675, 37.441125, 4326)),

( 48, geography::Point(55.660424, 37.434424, 4326)),

( 49, geography::Point(55.670701, 37.42598, 4326)),

( 50, geography::Point(55.67994, 37.418712, 4326)),

( 51, geography::Point(55.686873, 37.414868, 4326)),

( 52, geography::Point(55.695697, 37.407528, 4326)),

( 53, geography::Point(55.702805, 37.397952, 4326)),

( 54, geography::Point(55.709657, 37.388969, 4326)),

( 55, geography::Point(55.718273, 37.383283, 4326)),

( 56, geography::Point(55.728581, 37.378369, 4326)),

( 57, geography::Point(55.735201, 37.374991, 4326)),

( 58, geography::Point(55.744789, 37.370248, 4326)),

( 59, geography::Point(55.75435, 37.369188, 4326)),

( 60, geography::Point(55.762936, 37.369053, 4326)),

( 61, geography::Point(55.771444, 37.369619, 4326)),

( 62, geography::Point(55.779722, 37.369853, 4326)),

( 63, geography::Point(55.789542, 37.372943, 4326)),

( 64, geography::Point(55.79723, 37.379824, 4326)),

( 65, geography::Point(55.805796, 37.386876, 4326)),

( 66, geography::Point(55.814629, 37.390397, 4326)),

( 67, geography::Point(55.823606, 37.393236, 4326)),

( 68, geography::Point(55.83251, 37.395275, 4326)),

( 69, geography::Point(55.840376, 37.394709, 4326)),

( 70, geography::Point(55.850141, 37.393056, 4326)),

( 71, geography::Point(55.858801, 37.397314, 4326)),

( 72, geography::Point(55.867051, 37.405588, 4326)),

( 73, geography::Point(55.872703, 37.416601, 4326)),

( 74, geography::Point(55.877041, 37.429429, 4326)),

( 75, geography::Point(55.881091, 37.443596, 4326)),

( 76, geography::Point(55.882828, 37.459065, 4326)),

( 77, geography::Point(55.884625, 37.473096, 4326)),

( 78, geography::Point(55.888897, 37.48861, 4326)),

( 79, geography::Point(55.894232, 37.5016, 4326)),

( 80, geography::Point(55.899578, 37.513206, 4326)),

( 81, geography::Point(55.90526, 37.527597, 4326)),

( 82, geography::Point(55.907687, 37.543443, 4326)),

( 83, geography::Point(55.909388, 37.559577, 4326)),

( 84, geography::Point(55.910907, 37.575531, 4326)),

( 85, geography::Point(55.909257, 37.590344, 4326)),

( 86, geography::Point(55.905472, 37.604637, 4326)),

( 87, geography::Point(55.901637, 37.619603, 4326)),

( 88, geography::Point(55.898533, 37.635961, 4326)),

( 89, geography::Point(55.896973, 37.647648, 4326)),

( 90, geography::Point(55.895449, 37.667878, 4326)),

( 91, geography::Point(55.894868, 37.681721, 4326)),

( 92, geography::Point(55.893884, 37.698807, 4326)),

( 93, geography::Point(55.889094, 37.712363, 4326)),

( 94, geography::Point(55.883555, 37.723636, 4326)),

( 95, geography::Point(55.877501, 37.735791, 4326)),

( 96, geography::Point(55.874698, 37.741261, 4326)),

( 97, geography::Point(55.862464, 37.764519, 4326)),

( 98, geography::Point(55.861979, 37.765992, 4326)),

( 99, geography::Point(55.850257, 37.788216, 4326)),

(100, geography::Point(55.850383, 37.788522, 4326)),

(101, geography::Point(55.844167, 37.800586, 4326)),

(102, geography::Point(55.832707, 37.822819, 4326)),

(103, geography::Point(55.828789, 37.829754, 4326)),

(104, geography::Point(55.821072, 37.837148, 4326)),

(105, geography::Point(55.811599, 37.838926, 4326)),

(106, geography::Point(55.802781, 37.840004, 4326)),

(107, geography::Point(55.793991, 37.840965, 4326)),

(108, geography::Point(55.785017, 37.841576, 4326))

Скрипт 1

 

Сделаем из них многоугольник. На протяжении трех предыдущих постов мы учились разными способами соединять точки в контуры и делать из них фигуры. Kirill написал комментарий, что "я пока только фантазировать могу, как это на практике использовать". На самом деле задача самая что ни на есть практическая, т.к. на практике любая территория чаще всего задается в виде многоугольника. Вернее сказать, в виде набора координат вершин, соединив которые в нужном порядке, можно получить контур территории, а закрасив внутреннюю область контура - саму территорию. В качестве примера сошлюсь на один из своих ранних постов - Как я затаскивал в таблицу карту. Короче, мы уже так долго до этого занимались соединением точек, что я не буду, с вашего позволения, подробно останавливаться на процессе. Используем, например, метод Рис.5 из позапрошлого поста.

 

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

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

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

                                    <exterior>

                                       <LinearRing>

                                         <posList>' + @s + '</posList>

                                       </LinearRing>

                                                      </exterior>

                                                      </Polygon>', 4326)

select @m

Скрипт 2

 

clip_image002

Рис.1

 

По-моему, похоже. Для пропорционального отображения в соответствии с выбранным SRID в комбобоксе Select projection следует выбрать Меркатор.

В этом году DevCon проходил в пансионате Покровское, что подле деревни Часцы (55.657042, 36.819745). 15 лет назад местом проведения конференции был Центральный институт повышения квалификации в славном городе Обнинске недалеко от памятника Курчатову (55.10419, 36.61894). Уставшие после официальной программы конференции, а также неофициальной части разработчики не смогли в ночи опознать выдающегося ученого и путем логических умозаключений пришли к мысли, что это, видимо, товарищ Обнин. После чего отправились купаться на реку Протву, по имени которой получило название вино протвейн. Бр-р. Обратите внимание, что в WKT-методах сначала идет долгота, затем широта, в отличие от конструктора Point, где все наоборот. Это сделано для удобства тренировки памяти.

 

declare @devcon2011 geography = geography::STPointFromText('POINT(36.819745 55.657042)', 4326)

declare @devcon1996 geography = geography::STPointFromText('POINT(36.61894 55.10419 )', 4326)

Скрипт 3

 

Для определения расстояния в SQL Server используется OGC-метод STDistance(). В Denali стало возможным его не только отмерить, но и проложить. Для построения кратчайшей геодезической между фигурами в SQL Server Denali можно использовать метод ShortestLineTo():

select @devcon1996.ShortestLineTo(@m).ToString()

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

LINESTRING EMPTY

Скрипт 4

 

Но что это? Как пустое множество? Разве он не в состоянии проложить кратчайшее расстояние от Обнинска до Москвы? Может, глюк в СТР3? Надо срочно написать в connect! Хотя, приводимые в документации примеры, вроде, работают. Тогда в чем затык в нашем случае?

 

Плоскость (Geometry) бесконечна, из-за этого на ней не существует двусмысленности, что считать многоугольником (Polygon) в случае замкнутой границы. Естественно, внутреннюю часть, ограниченную контуром, потому что все, что снаружи, - бесконечно. Поверхность Земли (Geography) конечна, и любой замкнутый контур делит ее на две конечные области. Какую из них считать внутренностью многоугольника? Воткнем в глобус спицу так, чтобы она проходила через центр Земли и была опоясана замкнутым контуром. Наденем стрелку на тот конец спицы, что торчит из северного полушария. Посмотрим на замкнутый контур с высоты спицыной стрелки. Если обходить по нему против часов, то, что останется по левую руку, и будет считаться внутренней областью многоугольника. Теперь вернемся к нашему многоугольнику под названием МКАД. Нумерация километров МКАД ведется по часовой стрелке, поэтому в том многоугольнике, который мы построили на Рис.1, внутренней областью будет весь остальной глобус, внешний по отношению к закрашенной территории. SSMS рисует неправильно. По уму следовало бы территорию внутри МКАД оставить белой, а закрасить деревню Часцы, город Обнинск, Париж, Нью-Йорк и остальное замкадье. Поэтому, например, получается пустой кратчайшая линия от точки @devcon1996 до многоугольника @m - куда ее тянуть, если эта точка и так уже находится внутри этого многоугольника? Надо было в Скрипте 2 сделать select … from #mkad order by km desc, чтобы точки обходились в обратном порядке. Исправить ситуацию можно при помощи метода ReorientObject(), который выворачивает многоугольник наизнанку, меняя местами его внутреннюю и внешнюю области.

 

set @m = @m.ReorientObject()

select @m.RingN(1).STUnion(@devcon1996.ShortestLineTo(@m)).STUnion(@devcon2011.ShortestLineTo(@m))

Скрипт 5

 

clip_image002[1]

Рис.2

 

Теперь все правильно, и мы видим на карте контур МКАД вместе с кратчайшими линиями до него от пансионата Покровское и ЦИПК г. Обнинск.

Прикрываясь хилой отмазкой о невозможности различить на сфере внешнюю и внутреннюю границы, гадский OGC не предусмотрел в спецификации для Geography аналогов плоских методов STExteriorRing() и STInteriorRingN(), а заодно и STBoundary(). Правильно, чего мелочиться. Из всего многообразия имеется только RingN(), да и тот не стандартный, а наш. Спасибо разработчикам SQL Server. Без него на Земле было бы совсем тоскливо. В отличие от евклидовой плоскости. В нашем случае граница состоит из одного контура, который и возвращает RingN(1).

Как легко догадаться, длиной кратчайшего отрезка ShortestLineTo() является традиционный STDistance(). В этой метрике расстояние дается в метрах.

 

select @devcon1996.ShortestLineTo(@m).STLength(), @devcon1996.STDistance(@m), @devcon2011.ShortestLineTo(@m).STLength(), @devcon2011.STDistance(@m)

Скрипт 6

 

clip_image004

Рис.3

 

Как мы видим, DevСon 15-летней давности значительно отстоит от нынешнего не только во времени, но и по расстоянию.

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

 

 

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

Comments

  • Anonymous
    January 01, 2003
    По-видимому, я был все-таки неправ на тему формализации внутренней области полигона. За полигон тупо считается тот кусок, у которого площадь меньше. Надо посмотреть, как она себя проявит в случае границы = экватор.