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!

Угнали домены

30 июня рано утром мне пришла пара писем:

Prepayment amount 500 RUR for bill #3248738
was removed from your account at WebNames.ru.

и, как ни странно,

Уважаемый komzpa!
Заявка N 3248738 от 30.06.2015
Предмет заявки:
Заказ услуги «Виртуальный сервер (тариф «VDS-1″)» для домена fhack.pro
Сумма: 500.00 RUR

Сумма в размере 500.00 RUR на оплату заявки списана с Вашего
лицевого счёта.
Остаток на Вашем лицевом счёте составляет: 26.79 RUR.

Следом за ним прилетело письмо:

Уважаемый komzpa!

По Вашей заявке стартован виртуальный выделенный сервер fhack.pro (IP адрес 80.252.22.108) по
тарифному плану «Виртуальный сервер (тариф «VDS-1″)».

==============================================================
ЛОГИН И ПАРОЛЬ АДМИНИСТРАТОРА СЕРВЕРА
==============================================================
Логин администратора: 80.252.22.108
Пароль администратора: *******

================================

и как финальный аккорд,

Здравствуйте, komzpa!

Нами (ООО «Регтайм») был получен и выполнен Ваш запрос на изменение e-mail адреса к Вашему аккаунту. E-mail адрес был успешно изменён.

Если Вы не отправляли такой запрос, как можно быстрее обратитесь в службу технической поддержки по адресу support@webnames.ru

Ничтоже сумняшеся, я написал по предложенному адресу:

Втр Июн 30 11:24:15 2015, me@komzpa.net писал:
> Здравствуйте,  я не изменял адрес, похоже, угнали аккаунт — срочно
> откатите
> последние действия
Добрый день.
К сожалению, это возможно только по письменному заявлению.
Форму такого заявления Вы можете найти по адресу:
https://www.webnames.ru/help/forgot_password
Спасибо за обращение.С уважением, Михаил Егоров
Группа по поддержке клиентов и партнеров
ООО «Регтайм» (Webnames.ru)
Email: support@webnames.ru

Круто, да?
Сначала «пишите как можно скорее», а потом «мы без письменного заявления ничего делать не будем». Напоминаю, живу я в Минске.
Обнаружилось, что заявление можно завести в ближайший офис Экстмедии. Написал, занёс. Приятный молодой человек сделал копию моих паспартов (того, который во Whois, и того, который теперишний), и обещал переслать всё как можно скорее. «А домены за это время с аккаунта не унесут?» — «…унесут.»
icon-car.pngKML-LogoFullscreen-LogoQR-code-logoGeoJSON-LogoGeoRSS-Logo
Гостиница Виктория Олимп

Карта загружается. Пожалуйста, подождите.

Гостиница Виктория Олимп 53.935418, 27.494026
Вчера — ура! — пришло письмо.

Здравствуйте, Пролесковский Дорофей Аркадьевич!

Сегодня нами (ООО «Регтайм») было рассмотрено Ваше письмо о восстановлении доступа к аккаунту komzpa.
Сообщаем Вам, что ссылка на сброс пароля была отправлена на Ваш новый контактный e-mail адрес komzpa@gmail.com
Срок действия этой ссылки — час. Если Вы не успеете ей воспользоваться, Вы можете запросить её повторно, указав Ваш контактный e-mail адрес в разделе «забыли пароль» нашего портала Webnames.ru


С уважением, Павел Глазунов
Группа по поддержке клиентов и партнеров
ООО «Регтайм» (Webnames.ru)
Email: support@webnames.ru
Как-то так.
Теперь у меня есть ненужная VPS на месяц. Надо кому? :)
Мораль: денег на счету у Webnames для автопродления держать нельзя, шанс, что угонят и бросят высок. А вот сделать что-то фатальное, вроде унести домен, не получится — злоумышленнику топать с вашим паспортом в офис будет куда сложнее.
По этому поводу посмотрел на продающиеся домены и сделал постгресячий using.vodka. Ждём коммитов от товарищей :)

Стройка полугодия

Прошлый год в деревне закончился вестью о том, что в старом доме в деревне упала крыша, утянув за собой стену.

Вот что бывает, если в доме не жить десятилетиями, и копать огород вместо ведения хозяйства.

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

Решили попробовать обкатать постройку такого дома — раз уж на участке должно что-то стоять, чтобы не забрали :)

icon-car.pngKML-LogoFullscreen-LogoQR-code-logoGeoJSON-LogoGeoRSS-Logo
Ялча, дом 57

Карта загружается. Пожалуйста, подождите.

Ялча, дом 57 53.780503, 28.269464

Началась стройка с истерики родителей — «там же аварийное здание! его надо срочно сносить!» Таким образом, срочно снесли старый дом:

(Чем хорош Петя — всё пишет на видео, эффект присутствия на стройке полный :)

Пока дом сносили, оформляли проект. Гимп и планшет — лучший способ скетчить такие вещи.

building-plan-v0

 

Дом диаметром 6 метров, с крыльцом, лесенкой на второй этаж, второй этаж для кошачьей лежанки, лайт-кухня и туалет.

Потом подбили бюджет. Потом подбили ещё раз. Оказалось хорошо, но много :)

После фичаката у дома не осталось открытой веранды, собственного туалета и водопровода, а также утепления — всё равно домик летний.

Накинули к бюджету ещё чуть-чуть и зафиксировали как верхнюю планку. Если что пойдёт не так, то будем ещё резать фичи, но дом достроим до законченного состояния :)

Таймлайн проекта, с учётом расходов, решили выставить в ровно полгода. То есть уже через пять месяцев можно будет тыкать пальцем и говорить, что можно было сделать ещё лучше :)

Сейчас вот оформили эскизный проект.

Ну и в лучших традициях: так как проект для меня — фановый, а для Пети — рекламно-тренировочный, все планы, хауту и чертежи будут складироваться на гитхаб, куда-то в окрестности https://github.com/RudenkoPeter. Пулл-реквесты приветствуются ;)

Все интересные вопросы можно задавать мне лично :)