Как я снова затаскивал в таблицу карту

Снова — потому что как-то больше 2-х лет назад я уже проделывал это упражнение. То был длительный многотрудный процесс:
blogs.technet.com/b/isv_team/archive/2010/01/18/3306462.aspx
blogs.technet.com/b/isv_team/archive/2010/01/21/3307201.aspx
blogs.technet.com/b/isv_team/archive/2010/01/23/3307719.aspx
blogs.technet.com/b/isv_team/archive/2010/01/24/3307804.aspx
С тех пор наука шагнула далеко вперед. В данном посте мы опять-таки загрузим в SQL Server карту нашей необъятной Родины, на этот раз гораздо проще и элегантней благодаря авторам карт, новым возможностям SQL Server и независимым разработчикам, которым всем большое спасибо. Нам понaдобятся:

1. Подходящая карта
2. SQL Server 2012 (можно Express)
3. Инструмент для конвертации карты из ее первоначального формата в SQL Server.

Большую часть работы возьмет на себя замечательная бесплатная утилита Shape2SQL, которую написал датский MVP Morten Ødegaard Nielsen, дай Бог ему здоровья, работающий GIS Lead Software Engineer в ESRI. Компания ESRI (Environmental Systems Research Institute) известна программными продуктами семейства ArcGIS, получившими широкое распространение в том числе в России. Векторный формат shape-файлов был представлен для ArcView GIS версии 2 в начале 90-х, и на сегодняшний день благодаря своей распространенности стал де-факто стандартом для обмена данными между геоинформационными системами. Как следует из названия, тула умеет грузить shape-файлы в SQL Server 2008 и выше. В этом посте используется 2012 из-за появившейся в нем поддержки геометрических/географических агрегатов — пригодится дальше. В 2008-м пришлось бы писать курсор.
Примечание 1. Помимо утилиты импорта .shp в SQL Server, в SqlSpatialTools v1.3.0 (348 kb) (build 3759)на сайте SharpGIS входит клиент визуализации карт SqlSpatial Query Tool на WPF и Silverlight, который гораздо богаче по своим возможностям, чем географическая панель результатов SSMS.
Примечание 2. Оба инструмента, как предупреждает автор, не предназначены для коммерческой эксплуатации, т.к. написаны им для собственных нужд или удовольствия.
Примечание 3. .NET-библиотека с открытым исходным кодом, умеющая парсить shape-файлы, находится здесь.

Остается найти shape-файл с административно-территориальным делением РФ. Карты регионов России в различных форматах, в том числе .shp, отыскались на сайте лаборатории GIS-Lab: этот слой распространяется по лицензии CC-BY, разрешающей любое использование, включая передачу, изменение, использование в коммерческих проектах, с условием упоминания авторства и наличием гиперссылки. Источниками данных для карт послужили Федеральная служба государственной регистрации, кадастра и картографии (Росреестр) — границы субъектов, OpenStreerMap Россия (OSM) — Сухопутная граница России, кроме Калининградской области + граница г.Москвы, VMap0 — береговая линия + границы Калининградской области и инициативная группа проекта GIS-Lab.info, за что им тоже выражается глубокая признательность.
Карты GIS-Lab приводятся в 3-х проекциях: WGS 1984 (GPS), Pulkovo-1942 и Albers-Siberia. Для отображения территории России лучше всего подходит третья:

Рис.1

чем первые две:

Рис.2

К сожалению, мне не удалось найти в sys.spatial_reference_systems равновеликую коническую проекцию Альберса с центральным меридианом = 1050 в.д. на основе сфероида Красовского, поэтому ее можно затащить в SQL Server как геометрию и использовать, например, в качестве средства визуализации отчетов. Первые две импортируются в тип географию (SRID=4326, 4284), и их можно использовать для картографических расчетов.

Для импорта карты нам понадобятся три файла:

Рис.3

Файлы .shp/.shx подаются на вход зачотной туле Shape2SQL, про которую говорилось выше. Указывается полное имя shape-файла, соединение с SQL Server, название таблицы, в которую будут импортированы данные из shape-файла, название геопространственной колонки в таблице, название колонки, которая будет выполнять роль первичного ключа. Если таблица не существовала, она создается. Импортировать данные в проекции Albers Siberia как географию не удается, т.к. SQL Server не знает такого SRIDa и неправильно преобразовывает координаты — Data projection or extent are outside the bounds of what is supported by SqlGeography type — поэтому импортируем их как геометрию. При желании можно тут же построить геопространственный индекс. Жмем кнопку Upload to Database, и таблица dbo.regions2010_albsib успешно создается и наполняется данными:

Рис.4

Можно сделать из нее select * from rfmap.dbo.regions2010_albsib и на закладке Spatial results в SSMS убедиться, что все ОК (см.Рис.1). Однако количество записей в таблице несколько превышает количество регионов России. Чтобы понять, в чем дело, нам понадобится третий файл — dbf, содержащий метаданные к карте. Как загрузить dbf в SQL Server, мы проходили в предыдущем посте. Создадим таблицу

CREATE TABLE dbo.regions2010_captions (
ID int IDENTITY(1,1) NOT NULL,
AREA float NULL,
PERIMETER float NULL,
region varchar(50) NULL
)

и воспользуемся мастером импорта/экспорта. В принципе, он автоматически связывает поля источника и назначения с одинаковыми именами, но лишний раз в этом убедиться можно, нажав

Рис.5

Если длина поля region заранее неизвестна, т.е. под рукой нет dbf-читалки, а парсить бинарным редактором заголовок лень, можно заложить максимальную длину, а после импорта урезать. Импорт dbf проходит успешно, и наряду с таблицей картографических объектов в БД SQL Server появляется таблица подписей к ним:

Рис.6

Записи в таблице картографических объектов regions2010_albsib связаны с таблицей подписей regions2010_captions исключительно физическим порядком следования, что не есть гуд, т.к. никто не обещал, что мастер импорта/экспорта переносит записи в том же порядке, как они лежали в dbf. Наилучшим выходом из этой ситуации было бы присутствие РК в dbf и в shp, что бы позволило их легко сопоставить и, кстати, избавило бы автора Shape2SQL от необходимости придумывать суррогатный ID. Но такие карты сданы. Обращает на себя внимание тот факт, что записей подписей 1506, а записей картографических объектов — 1505 (Рис.1), хотя в процессе загрузки на Рис.4 Shape2SQL тоже писал, что их 1506. Куда девалась еще одна? Здесь выясняется один серьезный незачет зачотной тулы, который состоит в том, что она не грузит пустые объекты. Она их пропускает, и записей в таблице картографических объектов может оказаться меньше, чем подписей, в результате чего соответствие по физическому порядку съезжает к чертовой бабушке, потому что как понять, какие по номеру записи из середины были выкинуты? Хорошо, что здесь присутствуют колонки area и perimeter, по которым можно вычислить запись, соответствующую пустому картографическому объекту в подписях:

Рис.7

После ее удаления из таблицы подписей delete from regions2010_captions where area = 0 появляется дырка. Хорошо, что в нашем случае это была предпоследняя запись, и номер последней можно просто поменять на 1505. Тем не менее честно пронумеруем подписи по-новой в порядке возрастания ID, как если бы удаленных записей было несколько, и они были бы разбросаны где-нибудь в середине.

alter table regions2010_captions add n int
;with cte as (select *, row_number() over (order by ID) n1 from regions2010_captions) update cte set n = n1
alter table regions2010_captions drop column ID

Сильно надеюсь, что таблицу regions2010_albsib Нильсен тоже получает в порядке физического сканирования объектов. Тогда ее по полю ID можно связать с таблицей regions2010_captions (поле n). Количество картографических объектов превышает фактическое количество регионов, потому что, как мы видим из подписей, один регион может состоять более, чем из одного полигона: Москва и Зеленоград, прибрежная территория и острова и т.д. Например, select s.geom from regions2010_captions c join regions2010_albsib s on s.ID = c.n where c.region = 'Чукотский автономный округ':

Рис.8

Значит, их нужно объединить, сагрегировав по признаку одного и того же region:

if exists (select 1 from sys.sequences where name = 'Seq') drop sequence Seq
create sequence Seq as int start with 1 increment by 1

select * from regions2010

В этом коде иллюстрируются две новые возможности SQL Server 2012: последовательности и геопространственные агрегаты.

Рис.9

Судя по расцветке похоже на правду. Значит, со связыванием подписей и полигонов мы не ошиблись.
Несмотря на то, что полигоны, относящиеся к одному региону, закрашиваются одним и тем же цветом, это отдельные объекты, однотипную геометрическую коллекцию которых составляет регион (мультиполигон). В этом можно убедиться, выполнив запрос select ID, Region, geom.STGeometryType(), geom.STNumGeometries() from regions2010:

Рис.10

Для их развертывания из коллекции в таблицу для выбранного региона я написал простую функцию:

if object_id('dbo.ufnGeometries', 'TF') is not null drop function dbo.ufnGeometries
go
create function dbo.ufnGeometries(@g geometry) returns @t table (id int identity(1, 1), g geometry) as begin
declare @i int = 1
while @i <= @g.STNumGeometries() begin
insert @t (g) values (@g.STGeometryN(@i))
set @i += 1
end
return
end

которая позволяет из агрегатной, по сути, таблицы regions2010 вновь вернуться к первоначальной regions2010_albsib:

select r.ID, r.Region, f.ID, f.g from regions2010 r cross apply dbo.ufnGeometries(geom) f

Рис.11

На самом деле, это не совсем она. Обратите внимание, что детальных полигонов уже 1503, а не 1505. Я объясняю это тем, что имелось две пары пересекающихся полигонов или три пересекающихся полигона в пределах одного региона, которые объединились в один. Давайте проверим:

with cte as (select s.ID, s.geom, c.region from regions2010_albsib s join regions2010_captions c on s.ID = c.n)
select cte1.ID, cte1.geom, cte2.ID, cte2.geom, cte1.geom.STIntersection(cte2.geom).STGeometryType() from cte cte1 join cte cte2 on cte1.ID < cte2.ID and cte1.region = cte2.region
where cte1.geom.STIntersects(cte2.geom) = 1

Действительно, всего имеются семь таких пересечений, из них 5 точечных, которые при объединении не схлопываются в один полигон, и две пары смежных полигонов, т.е. имеющих общую границу в виде ломаной. Каждая из этих пар объединилась в один полигон, таким образом, их количество сократилось на 2:

Рис.12

В завершении рассмотрим практическое действие метода Reduce(). Обратите, пожалуйста, внимание на то, что отрисовка карты regions2010 (Рис.9) происходит с задержкой 2-3 сек. Это вызвано значительным количеством точек — порядка тысяч или даже десятков тысяч — задающих контур каждого региона. Для карты того размера, что будет показываться в отчете, такая скрупулезность является лишней, неоправданно увеличивающей время рендеринга. Для огрубления используется метод Reduce(), действующий по алгоритму Дугласа-Пекера. Чем больше значение параметра, тем сильнее огрубление (мы разбирали толерантные методы здесь). Например, select ID, Region, geom.STNumPoints() [Points Total], geom.Reduce(100).STNumPoints() Reduce100, geom.Reduce(500).STNumPoints() Reduce500, geom.Reduce(1000).STNumPoints() Reduce1000 from regions2010

Рис.13

Таким образом, select ID, Region, geom.Reduce(1000) from regions2010 отрисовывается не в пример быстрее при сохранении видимого для данного масштаба карты качества:

Рис.14

Автор: alexejs