TIGER 2015 and OpenStreetMap difference, The SQL way

MapBox тут запустил очередную версию TileReduce. Одной из заявленных фишечек стало «вычисление разницы между TIGER и OpenStreetMap за 11 минут на лаптопе разработчика». Мне показалось, что это тривиально делается на постгресе, Вова Агафонкин подначил меня проверить.

Добываем TIGER

1. Качаем TIGER с census.gov:

% wget -r ftp://ftp2.census.gov/geo/tiger/TIGER2015/ROADS/
 Total wall clock time: 4h 0m 32s
 Downloaded: 3231 files, 4,1G in 3h 29m 31s (342 KB/s)

1. % cd ftp2.census.gov/geo/tiger/TIGER2015/ROADS
   % ls *.zip | parallel unzip -o

2. Приготовим табличку под тигера:

   shp2pgsql -s 4269:900913 -p tl_2015_01001_roads.shp tiger | psql gis gis -q -X

3. Убьём лишние констрейнты, чтобы параллельный импорт был быстрее:

   % psql gis gis
   alter table tiger DROP CONSTRAINT tiger_pkey;
   alter table tiger ALTER COLUMN geom TYPE geometry;

4. Импортнём всё в заготовленную таблицу в много потоков:

   % ls *.shp | parallel --eta "shp2pgsql -s 4269 -S -D -a {} tiger | psql gis gis -q -X || shp2pgsql -s 4269 -D -a {} tiger | psql gis gis -q -X"

5. Перепроецируем всё в одну проекцию, выбросим лишнее:

   [local] gis@gis=> create table tiger_merc as (select gid, ST_Simplify(ST_Transform(geom, 900913),3) as geom from tiger where fullname is not null  order by geom);
    SELECT 6254781
    Time: 163278,633 ms

Добываем OSM

1. Качаем NA с geofabrik:

   axel http://download.geofabrik.de/north-america-latest.osm.pbf Downloaded 6,9 Gigabyte in 21:46 seconds. (5572,74 KB/s)

2. Обновляем дамп осма:

   kom@thinkcat ~ % osmupdate --drop-author --drop-version --drop-relations north-america-latest.osm.pbf north-america-latest-new.osm.pbf 
   631,18s user 11,10s system 92% cpu 11:36,61 total

   # для тех, кто любит смотреть на прогрессбары - для просмотра прогресса потокового прожёвывания дампов можно пользоваться утилитой pv
   kom@thinkcat ~ % pv north-america-latest-new.osm.pbf | osmconvert --drop-relations --out-pbf - > north-america-latest-new1.osm.pbf

3. Пишем стиль для osm2pgsql, узнаём, что он не умеет импортировать только линии без точек, лайкаем баг: https://github.com/openstreetmap/osm2pgsql/issues/272
4. Забираем конфиги furry-sansa для импорта в базу и импортируем:

   kom@thinkcat furry-sansa % python furry.py config/roads.conf importAll 
   .....
   indexes on highway_osm_line created in 1371s
    Completed highway_osm_line
    Stopped table: highway_osm_ways in 5514s
    node cache: stored: 491756655(60.75%), storage efficiency: 75.04% (dense blocks: 424942, sparse nodes: 110109832), hit rate: 66.98%
   Osm2pgsql took 9027s overall

5. Компьютер на импорте начинает резко педалить, что читать вконтактик становится неудобно — переключаем планировщик ввода-вывода:

   echo cfq | sudo tee /sys/block/sda/queue/scheduler

6. Не забываем вернуть шедулер перед дальнейшей работой:

   echo deadline | sudo tee /sys/block/sda/queue/scheduler

   Посчитаем разницу

1. Для начала просто сделаем это в один поток, прицепимся к базе qgis и посмотрим, что всё хорошо:

2. create table tiger_not_in_osm as (select (ST_Dump(ST_Safe_Difference(t.geom, (select ST_Buffer(ST_Collect(way), 5, 1) from highway_osm_line l where ST_DWithin(l.way, t.geom, 5))))).geom as geom from (select geom from tiger_merc t join lateral (select way from highway_osm_line l where t.geom && l.way and ((t.geom <-> l.way) < 5) order by t.geom <-> l.way limit 1) z on true where z.way is null or (ST_HausdorffDistance(t.geom, z.way) > 5)) t);

   Для незнакомых с SQL / PostGIS:

   create table tiger_not_in_osm as (
       select (
                  ST_Dump( -- разбираем мультигеометрии на простые одиночные кусочки
                      ST_Safe_Difference( -- считаем разницу, безопасный вариант с https://github.com/wgnet/globalmap/tree/master/code/postgis_wrappers
                          t.geom,
                          ( -- подзапрос, выбираем все лежащие рядом линии и раздуваем на 5 метров
                              select ST_Buffer(ST_Collect(way), 5, 1) -- последний аргумент позволяет уменьшить детализацию кружочков на стыках линий
                              from highway_osm_line l
                              where ST_DWithin(l.way, t.geom, 5) -- из запроса может вернуться NULL, эта ситуация обрабатывается ST_Safe_Difference, но не обрабатывается ST_Difference - ST_Difference(geom, NULL) = NULL
                          )
                      )
                  )
              ).geom as geom
       from (select geom
             from tiger_merc t -- выбираем геометрии из тайгера
                 join lateral (select way -- выбираем геометрию из осма для каждой геометрии из тайгера
                               from highway_osm_line l -- OSM - по факту тайгер старых годов, многие геометрии просто совпадут
                               where
                                   t.geom && l.way -- только если они пересекаются bbox'ами
                                   and ((t.geom <-> l.way) < 5) -- только если между центрами bbox'ов меньше 5 метров
                               order by t.geom <-> l.way -- сортируем по расстоянию между центрами
                               limit 1) z on true
             where
                 z.way is null or -- нас интересуют только те объекты тайгера, для которых в осме ничего не нашлось
                 (ST_HausdorffDistance(t.geom, z.way) > 5) -- или нашедшееся больше, чем на 5 метров отличается
            ) t
   );

   Без распараллеливания запрос выполняется 20 с половиной минут — всего в два раза дольше, чем на TileReduce.

3. Параллелим на ParPSQL, клиенте Postgres для Slightly Big Data:

   kom@thinkcat par_psql % cat tiger.sql
   drop table if exists tiger_not_in_osm_1;
   drop table if exists tiger_not_in_osm_2;
   drop table if exists tiger_not_in_osm_3;
   drop table if exists tiger_not_in_osm_4;
   drop table if exists tiger_not_in_osm_parallel;
   create unlogged table tiger_not_in_osm_1 as (with tiger_id_range as (select min(gid) as min_gid, max(gid) as max_gid from tiger_merc) select (ST_Dump(ST_Safe_Difference(t.geom, (select ST_Buffer(ST_Collect(way), 10, 1) from highway_osm_line l where ST_DWithin(l.way, t.geom, 10))))).geom as geom from (select geom from (select geom from tiger_merc, tiger_id_range where gid between min_gid+0*(max_gid - min_gid)/4 and min_gid+1*(max_gid - min_gid)/4) t join lateral (select way from highway_osm_line l where t.geom && l.way and ((t.geom <-> l.way) < 10) order by t.geom <-> l.way limit 1) z on true where (z.way is null or (ST_HausdorffDistance(t.geom, z.way) > 10))) t ); --&
   create unlogged table tiger_not_in_osm_2 as (with tiger_id_range as (select min(gid) as min_gid, max(gid) as max_gid from tiger_merc) select (ST_Dump(ST_Safe_Difference(t.geom, (select ST_Buffer(ST_Collect(way), 10, 1) from highway_osm_line l where ST_DWithin(l.way, t.geom, 10))))).geom as geom from (select geom from (select geom from tiger_merc, tiger_id_range where gid between min_gid+1*(max_gid - min_gid)/4 and min_gid+2*(max_gid - min_gid)/4) t join lateral (select way from highway_osm_line l where t.geom && l.way and ((t.geom <-> l.way) < 10) order by t.geom <-> l.way limit 1) z on true where (z.way is null or (ST_HausdorffDistance(t.geom, z.way) > 10))) t ); --&
   create unlogged table tiger_not_in_osm_3 as (with tiger_id_range as (select min(gid) as min_gid, max(gid) as max_gid from tiger_merc) select (ST_Dump(ST_Safe_Difference(t.geom, (select ST_Buffer(ST_Collect(way), 10, 1) from highway_osm_line l where ST_DWithin(l.way, t.geom, 10))))).geom as geom from (select geom from (select geom from tiger_merc, tiger_id_range where gid between min_gid+2*(max_gid - min_gid)/4 and min_gid+3*(max_gid - min_gid)/4) t join lateral (select way from highway_osm_line l where t.geom && l.way and ((t.geom <-> l.way) < 10) order by t.geom <-> l.way limit 1) z on true where (z.way is null or (ST_HausdorffDistance(t.geom, z.way) > 10))) t ); --&
   create unlogged table tiger_not_in_osm_4 as (with tiger_id_range as (select min(gid) as min_gid, max(gid) as max_gid from tiger_merc) select (ST_Dump(ST_Safe_Difference(t.geom, (select ST_Buffer(ST_Collect(way), 10, 1) from highway_osm_line l where ST_DWithin(l.way, t.geom, 10))))).geom as geom from (select geom from (select geom from tiger_merc, tiger_id_range where gid between min_gid+3*(max_gid - min_gid)/4 and min_gid+4*(max_gid - min_gid)/4) t join lateral (select way from highway_osm_line l where t.geom && l.way and ((t.geom <-> l.way) < 10) order by t.geom <-> l.way limit 1) z on true where (z.way is null or (ST_HausdorffDistance(t.geom, z.way) > 10))) t ); --&
   create table tiger_not_in_osm_parallel as (select * from tiger_not_in_osm_1 union all select * from tiger_not_in_osm_2 union all select * from tiger_not_in_osm_3 union all select * from tiger_not_in_osm_4);%

4. Запускаем:

5. ./par_psql gis gis --file=tiger.sql

6. Через 10 с половиной минут наслаждаемся сгенерированной таблицей дорог, которые есть в TIGER, но нет в OSM. PROFIT!