TIGER 2015 and OpenStreetMap difference, The SQL way

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

osm-tiger

Добываем 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!

TIGER 2015 and OpenStreetMap difference, The SQL way: 2 комментария

    1. Параллельные сканы — это прекрасно, ждём :)
      Но в посте они, в принципе, сэмулированы четырьмя параллельными селектами, так что по идее должны быть просто синтаксическим сахаром. Или нет?

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *