Приоритет отрисовки в OpenStreetMap

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Эта страница является черновиком статьи.


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

Приоритет отрисовки в OpenStreetMap

Приоритет отрисовки в OpenStreetMap

Приоритет отрисовки в OpenStreetMap


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

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

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

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

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

Исходные данные

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

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

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 в этой проекции и перейдем к данной location: 

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


Затем создадим region на 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


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

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


или в текст (перед тем как импортировать в 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.

Используя location в проекции Web Mercator с region на уровне 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 построим mask, сдвинем и масштабируем значения к диапазону [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

Взвешенное расстояние до центра 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

Разрешение спутниковых снимков 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

Источник — https://wiki.gis-lab.info/index.php?title=Приоритет_отрисовки_в_OpenStreetMap&oldid=25289