Перевод статьи, в которой при помощи GRASS GIS и PostGIS вычисляются плохо отрисованные в OpenStreetMap малонаселенные области.

Оригинал: OpenStreetMap Mapping Priority

[править] Приоритет отрисовки в OpenStreetMap

Работая в Французской Гвиане, автор узнал о проекте Mapazonia, целью которого является улучшение данных OpenStreetMap в регионе Амазонии. С начала участия в проекте автор научился использовать Менеджер задач OSM - инструмент, который позволяет разбить AOI (Areas Of Interest - область интереса) на отдельные тайлы задач. Будучи заинтересованным в отрисовке Гвианского плоскогорья, автор задался вопросом: возможно ли использовать подобный подход для этой области?

Как только тайлы задач были созданы (их насчиталось тысячи для такой большой области), возник вопрос: где наносить объекты в первую очередь?

Чтобы ответить на этот вопрос, нужно учесть несколько подходов:

• мы должны найти неотрисованные объекты там, где нет данных;
• пустые тайлы в удаленных областях, где живёт немного людей, имеют меньший приоритет;
• там, где исходные данные (спутниковые снимки) в низком разрешении, только объекты, покрывающие большие площади, могут быть оцифрованы (так называемый small scale маппинг);
• вносить данные нужно в первую очередь в центре AOI.

А теперь давайте построим схему приоритета отрисовки тайлов, используя эти идеи!

[править] Исходные данные

Начнем с получения исходных данных. В итоге мы должны получить схему, похожую на работу Мартина Райфера (Martin Raifer) по плотности узлов. Данные OSM доступны на сайте Geofabrik, и мы можем импортировать их в базу данных PostGIS, следуя инструкциям из руководства Регины Обе (Regina Obe). Тайлы задач соответствуют 12 уровню тайлов в проекции Web Mercator.

Создадим таблицу тайлов:

CREATE TABLE tile AS SELECT
  x, y
FROM
  generate_series(0, pow(2, 12)::INT - 1) AS x,
  generate_series(0, pow(2, 12)::INT - 1) AS y ;
COMMENT ON TABLE tile IS 'Map tiles' ;
COMMENT ON COLUMN tile.x IS 'Tile coordinate' ;
COMMENT ON COLUMN tile.y IS 'Tile coordinate' ;

и дополнительно построим их геометрии (таким образом мы будем использовать векторную и растровую модель представления данных):

CREATE FUNCTION ST_GeomFromTileXY(z int, x int, y int, OUT tile_geom geometry(Polygon, 4326)) AS $$
DECLARE
    side REAL;
    xmin REAL; ymin REAL; xmax REAL; ymax REAL;
    minlon REAL; maxlon REAL; maxlat REAL; minlat REAL;
BEGIN
    side := pow(2, z);
 
    xmin := x::REAL/side - .5;
    ymin := .5 - y::REAL/side;
    xmax := (x + 1)::REAL/side - .5;
    ymax := .5 - (y + 1)::REAL/side;
 
    minlon := 360.*xmin;
    maxlon := 360.*xmax;
    maxlat := 90. - 360.*ATAN(EXP(-ymin*2.*pi()))/pi();
    minlat := 90. - 360.*ATAN(EXP(-ymax*2.*pi()))/pi();
 
    tile_geom := ST_MakeEnvelope(minlon, minlat, maxlon, maxlat, 4326);
END; $$ LANGUAGE plpgsql;

Функция для получения X и Y координат тайла для заданной точки:

CREATE FUNCTION TileCoord(z int, point geometry(Point, 4326),
  OUT tile_x int, OUT tile_y int) AS $$
DECLARE
    tile_n REAL;
    sinphi REAL; x REAL; y REAL;
BEGIN
    tile_n := POWER(2, z)::REAL;
 
    x := (ST_X(point) + 180.)/360.;
    sinphi = SIN(ST_Y(point)*pi()/180.);
    y := .5 - LN((1. + sinphi)/(1. - sinphi))/(4.*pi());
 
    tile_x := CAST(TRUNC(tile_n*x + .5) AS int);
    tile_y := CAST(TRUNC(tile_n*y + .5) AS int);
END; $$ LANGUAGE plpgsql;

Для вычисления плотности узлов сначала создадим временную таблицу, в которой будет содержаться информация о координатах тайла, в который попадает каждая точка:

CREATE TEMP TABLE planet_osm_point_dump AS
SELECT id, TYPE, (coord).tile_x, (coord).tile_y, geom
  FROM
    (SELECT
      ROW_NUMBER() OVER (ORDER BY TYPE, id, (DUMP).PATH) AS id,
      TYPE,
      TileCoord(12, (DUMP).geom) AS coord,
      (DUMP).geom AS geom
    FROM
        (SELECT 0 AS TYPE, osm_id AS id, ST_DumpPoints(way) AS DUMP
        FROM planet_osm_point
      UNION ALL
        SELECT 1 AS TYPE, osm_id AS id, ST_DumpPoints(way) AS DUMP
        FROM planet_osm_line
      UNION ALL
        SELECT 2 AS TYPE, osm_id AS id, ST_DumpPoints(way) AS DUMP
        FROM planet_osm_polygon
        ) AS point
    ) AS point_tiled;

и затем посчитаем число узлов в каждом тайле нашей AOI:

ALTER TABLE tile ADD COLUMN node_count INT DEFAULT 0;
 
UPDATE tile t SET node_count = COUNT.num
FROM (
  SELECT x, y, COUNT(*) AS num
  FROM planet_osm_point_dump
  GROUP BY x, y ORDER BY x, y) AS COUNT
WHERE 
  t.x = COUNT.x AND t.y = COUNT.y;

[править] Малонаселенные области

Очевидно, что плотность населения выше вдоль побережья и главных рек, уменьшаясь вверх по течению реки. Мы будем моделировать удаленность, используя гидрологическую дистанцию, которая определяется поверхностным дренажём. Вычисления производятся в модуле r.stream.distance системы GRASS GIS. В гидрологическом моделировании обычно используется схема направлений потоков, полученная из цифровой модели местности (Digital Elevation Model - DEM):

Направления гидрологических потоков

USGS предоставляет набор данных GTOPO30 HYDRO1k, который можно скачать с сайта Earth Explorer (ID GT30H1KSA для Южной Америки). Набор данных представлен в азимутальной равновеликой проекции Ламберта (LAEA), которой соответствуют следующие параметры PROJ.4:

+proj=laea +ellps=sphere +lon_0=-60 +lat_0=-15 +x_0=0.0 +y_0=0.0 +units=m

Исходя из этих параметров создадим location с использованием модулей g.proj и mapset:

g.proj -c proj4="+proj=laea +ellps=sphere +lon_0=-60 +lat_0=-15 +x_0=0.0 +y_0=0.0 +units=m" location=south_america
g.mapset -c mapset=hydro location=south_america

Набор данных содержит несколько слоёв (см. файл README), но нам в основном интересен импорт растра высот и направлений течения:

r.in.gdal    in=gt30h1ksa/sa_dem.bil out=elev
r.in.gdal -e in=gt30h1ksa/sa_fd.bil  out=dir_arc

Преобразуем коды направлений в формат GRASS GIS, используя таблицу сопоставления в файле arc2grass.lut:

255 = 0
128 = 1
64  = 2
32  = 3
16  = 4
8   = 5
4   = 6
2   = 7
1   = 8

r.reclass in=dir_arc out=dir rules=arc2grass.lut
r.stream.distance -o dir=dir elev=elev method=up dist=dir_dist

Расстояние до русла

Для того, чтобы назначить расстояния тайлам, мы перепроецируем сетку из проекции LAEA в проекцию Web Mercator. На 15-м уровне длина стороны тайла примерно равна 1223 м, что наиболее близко к разрешению нашей сетки расстояний. Сначала создадим локацию (прим. переводчика: "location" в терминологии GRASS GIS) в этой проекции и перейдем к данной location:

g.proj -c epsg=3857 location=mercator
g.mapset -c mapset=priority location=mercator

Затем создадим регион (прим. переводчика: "region" в терминологии GRASS GIS) на 15-м уровне и перепроецируем:

g.region w=-20037508.342789244 s=-20037508.342789244 e=20037508.342789244 n=20037508.342789244 rows=32768 cols=32768
r.proj in=dir_dist loc=south_america mapset=hydro out=dir_dist_1k meth=near

Теперь установим для региона разрешение тайлов и выполним ресэмплинг сетки:

g.region w=-20037508.3427892 s=-20037508.3430388 e=20037508.3427892 n=20037508.3430388 rows=4096 cols=4096
r.resamp.stats in=dir_dist_1k out=dir_dist meth=median

[править] Экспорт

В конечном итоге мы получили данные, которые можем экспортировать в растровый файл:

r.mapcalc "dir_dist = if(isnull(dir_dist), -1, int(dir_dist))"
r.out.gdal -t  in=dir_dist out=hydro_dist.tif type=UInt32 nodata=4294967295

или в текстовый файл XYZ (перед импортом в PostGIS):

r.mapcalc "y = row() - 1"
r.mapcalc "x = col() - 1"
r.out.xyz in=y        sep=space out=hydro_dist-y.xyz
r.out.xyz in=x        sep=space out=hydro_dist-x.xyz
r.out.xyz in=dir_dist sep=space out=hydro_dist.xyz

paste hydro_dist-x.xyz hydro_dist-y.xyz hydro_dist.xyz | awk '!/-1$/ {print $3, $6, $9}' > hydro_dist.ssv

Обратите внимание, что NULL или MASK ячейки не будут экспортированы (см. r.out.xyz documentation), что может привести процедуру импорта в базу данных к неожиданному результату.

CREATE TEMP TABLE tile_attr (
z INTEGER DEFAULT 12, x INTEGER, y INTEGER,
VALUE REAL);

COPY tile_attr (x, y, VALUE) FROM 'hydro_dist.ssv' CSV DELIMITER ' ' ;
UPDATE tile t SET hydro_dist = a.VALUE
FROM tile_attr a
WHERE t.LEVEL = a.z AND t.x = a.x AND t.y = a.y ;

[править] Расстояние до AOI

Чтобы увеличить приоритет тайлов, которые находятся в центре нашей AOI, построим сетку весов. Вес в данном случае - это расстояние по прямой до границ AOI.

Используя локацию в проекции Web Mercator с регионом на 12-м тайловом уровне, импортируем (используя v.proj при необходимости) и растеризуем нашу область интересов (AOI):

v.to.rast in=aoi type=line,area out=aoi use=val value=1

Рассчитаем расстояние до AOI для внешних ячеек:

r.grow.distance in=aoi dist=dist_outer_aoi

Выполним то же самое для внутренних ячеек:

r.mapcalc "aoi_inv = if(isnull(aoi), 1, null())"
r.grow.distance in=aoi_inv dist=dist_inner_aoi

Комбинируем сетки расстояний в инвертированных значениях для внутренних тайлов:

r.mapcalc "dist_aoi = if(isnull(aoi), dist_outer_aoi, -dist_inner_aoi)"

В случае использования нескольких AOI построим растровую маску, сдвинем и масштабируем значения к диапазону [0, 1]:

r.mapcalc "roi = aoi_1 ||| aoi_2"
r.grow in=roi out=roi old=1 new=1
r.mask roi
eval $(r.univar -g dist_aoi | grep "\(min\|max\)")
r.mapcalc "dist_aoi = (dist_aoi - $min)/($max - $min)" --o

Вычислим средневзвешенное значение для каждого расстояния:

r.mapcalc "aois_dist = (dist_aoi_1 + dist_aoi_2)/2."

Взвешенное расстояние до центра AOI

[править] Разрешение снимков

Информация о разрешении снимков может быть получена, например, в инструменте анализа покрытия спутниковых снимков Bing (см. статью Coverage на OSM Wiki).

Ниже пример на Python:

PATTERN = "http://t{srv}.domain.tld/tile/{key}.jpg"
HEADER = "Content-Length"; server = randint(0, 7)
url = PATTERN.format(srv=server, key=quadkey)
meta = urlopen(url).info()
size = meta.getheaders(HEADER)[0]

Функция PostgreSQL, вычисляющая индекс quadkey для заданного уровня (zoom level) и координат тайла (x, y):

CREATE FUNCTION QuadKey(z int, x int, y int, OUT quadkey character varying(16)) AS $$
DECLARE digit int; mask int; BEGIN
    quadkey := '';
 
    FOR i IN REVERSE z..1 LOOP
        digit := 0;
        mask  := 1 << (i - 1);
 
        IF (x & mask) != 0 THEN digit := digit + 1; END IF;
        IF (y & mask) != 0 THEN digit := digit + 2; END IF;
 
        quadkey := quadkey || digit;
    END LOOP;
END; $$ LANGUAGE plpgsql;

В областях с низким разрешением снимков Bing тайлы на высоких масштабных уровнях (zoom level) в основном отсутствуют. Давайте оценим качество покрытия снимками высокого разрешения путём подсчета ненулевых тайлов на 14-м уровне:

SELECT
  substring(quadkey FOR 12) AS key,
  COUNT(quadkey)/16 AS resolution_index
FROM tile_bing WHERE z = 14 AND size > 0
GROUP BY substring(quadkey FOR 12) ORDER BY key ;

Разрешение спутниковых снимков Bing

[править] Создание индекса приоритетности

Используя полученные ранее измерения для тайлов, мы можем создать индекс приоритетности в соответствии с нашими пожеланиями. Например, формула

priority = ((1 - density) + dist_hydro + (1 - dist_aoi) + (1 - resolution))/4

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

Для того, чтобы использовать схему приоритетов в QGIS, добавьте столбец status в таблицу tile и создайте в плагине DB Manager динамически обновляемый слой TODO на основе SQL-запроса:

SELECT
  ROW_NUMBER() OVER (ORDER BY priority DESC) AS id,
  priority, x, y, ST_ExteriorRing(geom) AS geom
FROM tile WHERE status <> 'DONE'
ORDER BY priority DESC
LIMIT 32 ;

Данный слой показывает следующие 32 тайла задач для отрисовки:

Тайлы задач с наивысшим приоритетом

Перед тем, как начать наносить объекты в JOSM, либо сохраните выбранный тайл в файл формата GPX (отметьте опцию FORCE_GPX_TRACK и необходимость использования ST_ExteriorRing), или вызовите GDAL из layer action:

ogr2ogr -f GPX -sql "SELECT x || ' ' || y AS name, ST_ExteriorRing(geom) FROM tile WHERE x = '[% "x" %]' AND y = '[% "y" %]'" [% "x" %]_[% "y" %].gpx PG:dbname=osm -lco FORCE_GPX_TRACK=YES -nlt LINESTRING

Когда все интересующие вас объекты, входящие в границы тайла задания, добавлены, обновите атрибут его статуса в базе данных (status = DONE) и обновите экран QGIS (нажав клавишу F5).

Используя данную схему приоритетов, автор замечал, что иногда снимки Mapbox в более высоком разрешении, чем снимки Bing. Было бы полезно выяснить, предоставляют ли они схему разрешения снимков. Кроме того, две плитки последовательного ранга могут быть пространственно сильно удалены друг от друга, и было бы интересно найти способ их использования кластерами. Наконец, не следует ли включить приоритет непосредственно в Менеджер задач OSM?

Успешного Вам маппинга!

[править] Источники

Оригинал статьи: OpenStreetMap Mapping Priority

Автор статьи: adrienandrem

Автор перевода: Sibri