Приоритет отрисовки в OpenStreetMap: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Нет описания правки
 
(не показаны 4 промежуточные версии этого же участника)
Строка 1: Строка 1:
{{Статья|Черновик}}
{{Статья|Опубликована|osm-mapping-priority}}


{{Аннотация|Перевод статьи, в которой при помощи GRASS GIS и PostGIS вычисляются плохо отрисованные в OpenStreetMap малонаселенные области.}}
{{Аннотация|Перевод статьи, в которой при помощи GRASS GIS и PostGIS вычисляются плохо отрисованные в OpenStreetMap малонаселенные области.}}
Строка 7: Строка 7:
== Приоритет отрисовки в OpenStreetMap ==
== Приоритет отрисовки в OpenStreetMap ==


[[Файл:Osmmp_mapping_priority_map.jpeg|700|Приоритет отрисовки в OpenStreetMap]]
[[Файл:Osmmp_mapping_priority_map.jpeg|700px|Приоритет отрисовки в OpenStreetMap]]




Работая в Французской Гвиане, автор узнал о проекте [http://mapazonia.org/ Mapazonia], целью которого является улучшение данных  [http://www.openstreetmap.org/ OpenStreetMap] в регионе [http://wiki.openstreetmap.org/wiki/File:Amazonrivermap.svg Амазонии]. С начала участия в проекте автор научился использовать [http://tasks.hotosm.org/?sort_by=priority&direction=asc&search=mapazonia Менеджер задач OSM] - инструмент, который позволяет разбить AOI (Areas Of Interest - область интереса) на отдельные участки по задачам. Будучи заинтересованным в отрисовке [http://ru.wikipedia.org/wiki/Гвианское_плоскогорье Гвианского плоскогорья], автор задался вопросом: возможно ли использовать подобный подход для этой области?
Работая в Французской Гвиане, автор узнал о проекте [http://mapazonia.org/ Mapazonia], целью которого является улучшение данных  [http://www.openstreetmap.org/ OpenStreetMap] в регионе [http://wiki.openstreetmap.org/wiki/File:Amazonrivermap.svg Амазонии]. С начала участия в проекте автор научился использовать [http://tasks.hotosm.org/?sort_by=priority&direction=asc&search=mapazonia Менеджер задач OSM] - инструмент, который позволяет разбить AOI (Areas Of Interest - область интереса) на отдельные тайлы задач. Будучи заинтересованным в отрисовке [http://ru.wikipedia.org/wiki/Гвианское_плоскогорье Гвианского плоскогорья], автор задался вопросом: возможно ли использовать подобный подход для этой области?


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


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


Начнем с получения исходных данных. В итоге мы должны получить схему, похожую на [http://www.openstreetmap.org/user/tyr_asd/diary/22363 работу Мартина Райфера (Martin Raifer) по плотности узлов]. Данные OSM доступны на сайте [http://download.geofabrik.de/south-america-latest.osm.bz2 Geofabrik], и мы можем импортировать их в базу данных [http://postgis.net/ PostGIS], следуя инструкциям из [http://www.bostongis.com/?content_name=loading_osm_postgis#229 руководства Регины Обе (Regina Obe)]. Тайлы задач соответствуют [http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Zoom_levels 12 уровню] тайлов в [http://wiki.openstreetmap.org/wiki/EPSG:3857 проекции веб-картографии].
Начнем с получения исходных данных. В итоге мы должны получить схему, похожую на [http://www.openstreetmap.org/user/tyr_asd/diary/22363 работу Мартина Райфера (Martin Raifer) по плотности узлов]. Данные OSM доступны на сайте [http://download.geofabrik.de/south-america-latest.osm.bz2 Geofabrik], и мы можем импортировать их в базу данных [http://postgis.net/ PostGIS], следуя инструкциям из [http://www.bostongis.com/?content_name=loading_osm_postgis#229 руководства Регины Обе (Regina Obe)]. Тайлы задач соответствуют [http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Zoom_levels 12 уровню] тайлов в [http://wiki.openstreetmap.org/wiki/EPSG:3857 проекции Web Mercator].


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




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


<syntaxhighlight lang="plsql">
<syntaxhighlight lang="plsql">
Строка 130: Строка 130:
== Малонаселенные области ==
== Малонаселенные области ==


По всей видимости плотность населения выше вдоль побережья и главных рек, уменьшаясь вверх по течению реки. Мы будем моделировать '''удаленность''' используя '''гидрологическую дистанцию''', которая определяется [http://ru.wikipedia.org/wiki/Поверхностный_дренаж поверхностным дренажем]. Вычисления производятся в модуле [http://grasswiki.osgeo.org/wiki/R.stream.*#r.stream.distance r.stream.distance] системы [http://grass.osgeo.org/ GRASS GIS]. В гидрологическом моделировании обычно используется схема '''направлений потоков''', полученная из  [http://en.wikipedia.org/wiki/Digital_elevation_model цифровой модели местности] (Digital Elevation Model - DEM):  
Очевидно, что плотность населения выше вдоль побережья и главных рек, уменьшаясь вверх по течению реки. Мы будем моделировать '''удаленность''', используя '''гидрологическую дистанцию''', которая определяется [http://ru.wikipedia.org/wiki/Поверхностный_дренаж поверхностным дренажём]. Вычисления производятся в модуле [http://grasswiki.osgeo.org/wiki/R.stream.*#r.stream.distance r.stream.distance] системы [http://grass.osgeo.org/ GRASS GIS]. В гидрологическом моделировании обычно используется схема '''направлений потоков''', полученная из  [http://en.wikipedia.org/wiki/Digital_elevation_model цифровой модели местности] (Digital Elevation Model - DEM):  


[[Файл:Osmmp_hydrology_map.jpeg|700|Направления гидрологических потоков]]
[[Файл:Osmmp_hydrology_map.jpeg|700px|Направления гидрологических потоков]]


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


USGS предоставляет набор данных [http://lta.cr.usgs.gov/HYDRO1K GTOPO30 HYDRO1k], который можно скачать с сайта [http://earthexplorer.usgs.gov/ Earth Explorer] (ID GT30H1KSA для Южной Америки). Набор данных представлен в азимутальной равновеликой проекции Ламберта (LAEA), которой соответствуют следующие параметры [http://proj4.org/ PROJ.4]:
USGS предоставляет набор данных [http://lta.cr.usgs.gov/HYDRO1K GTOPO30 HYDRO1k], который можно скачать с сайта [http://earthexplorer.usgs.gov/ Earth Explorer] (ID GT30H1KSA для Южной Америки). Набор данных представлен в азимутальной равновеликой проекции Ламберта (LAEA), которой соответствуют следующие параметры [http://proj4.org/ PROJ.4]:
Строка 141: Строка 140:
<pre>+proj=laea +ellps=sphere +lon_0=-60 +lat_0=-15 +x_0=0.0 +y_0=0.0 +units=m</pre>
<pre>+proj=laea +ellps=sphere +lon_0=-60 +lat_0=-15 +x_0=0.0 +y_0=0.0 +units=m</pre>


Исходя из этих параметров создадим ''location'' с использованием [http://grass.osgeo.org/grass71/manuals/g.proj.html g.proj] и ''mapset'':
Исходя из этих параметров создадим ''location'' с использованием модулей [http://grass.osgeo.org/grass71/manuals/g.proj.html g.proj] и ''mapset'':


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 149: Строка 148:




Набор данных содержит несколько слоев (смотрите [http://lta.cr.usgs.gov/HYDRO1KReadMe файл README]), но нам в основном интересен импорт высот и направлений течения:
Набор данных содержит несколько слоёв (см. [http://lta.cr.usgs.gov/HYDRO1KReadMe файл README]), но нам в основном интересен импорт растра высот и направлений течения:


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 178: Строка 177:
<br />
<br />


[[Файл:Osmmp hydrology distance map.jpeg|700|Расстояние до русла]]
[[Файл:Osmmp hydrology distance map.jpeg|700px|Расстояние до русла]]


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


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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 192: Строка 190:




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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 200: Строка 198:




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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 206: Строка 204:
r.resamp.stats in=dir_dist_1k out=dir_dist meth=median
r.resamp.stats in=dir_dist_1k out=dir_dist meth=median
</syntaxhighlight>
</syntaxhighlight>


== Экспорт ==
== Экспорт ==
Строка 218: Строка 215:




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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 233: Строка 230:




Обратите внимание, что NULL или MASK ячейки не будут экспортированы (см. [http://grass.osgeo.org/grass71/manuals/r.out.xyz.html r.out.xyz documentation]) что может привести процедуру вставки к неожиданному результату.
Обратите внимание, что NULL или MASK ячейки не будут экспортированы (см. [http://grass.osgeo.org/grass71/manuals/r.out.xyz.html r.out.xyz documentation]), что может привести процедуру импорта в базу данных к неожиданному результату.


<syntaxhighlight lang="plsql">
<syntaxhighlight lang="plsql">
Строка 253: Строка 250:
Чтобы увеличить приоритет тайлов, которые находятся в центре нашей AOI, построим сетку весов. Вес в данном случае - это расстояние по прямой до границ AOI.
Чтобы увеличить приоритет тайлов, которые находятся в центре нашей AOI, построим сетку весов. Вес в данном случае - это расстояние по прямой до границ AOI.


Используя location в проекции Web Mercator с region на уровне 12 импортируйте (используя [http://grass.osgeo.org/grass71/manuals/v.proj.html v.proj] при необходимости) и растеризуйте вашу AOI:
Используя локацию в проекции Web Mercator с регионом на 12-м тайловом уровне, импортируем (используя [http://grass.osgeo.org/grass71/manuals/v.proj.html v.proj] при необходимости) и растеризуем нашу область интересов (AOI):


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 263: Строка 260:


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
r.grow.distance in=aoi     dist=dist_outer_aoi
r.grow.distance in=aoi dist=dist_outer_aoi
</syntaxhighlight>
</syntaxhighlight>


Строка 282: Строка 279:




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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 293: Строка 290:




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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 300: Строка 297:




[[Файл:Osmmp weighted distance map.jpeg|700|Взвешенное расстояние до центра AOI]]
[[Файл:Osmmp weighted distance map.jpeg|700px|Взвешенное расстояние до центра AOI]]


''Взвешенное расстояние до центра AOI''
''Взвешенное расстояние до центра AOI''
Строка 307: Строка 304:
== Разрешение снимков ==
== Разрешение снимков ==


Информация о разрешении снимков может быть получена, например, как в инструменте [http://mvexel.dev.openstreetmap.org/bingimageanalyzer/?lat=0.882&lon=-71.132&zoom=8 анализа покрытия спутниковых снимков Bing] (смотрите статью [http://wiki.openstreetmap.org/wiki/Bing_Maps/Coverage Coverage] на OSM Wiki).
Информация о разрешении снимков может быть получена, например, в [http://mvexel.dev.openstreetmap.org/bingimageanalyzer/?lat=0.882&lon=-71.132&zoom=8 инструменте анализа покрытия спутниковых снимков Bing] (см. статью [http://wiki.openstreetmap.org/wiki/Bing_Maps/Coverage Coverage] на OSM Wiki).




Строка 321: Строка 318:




и функция PostgreSQL, вычисляющая индекс [http://wiki.openstreetmap.org/wiki/QuadTiles quadkey] для заданного уровня (zoom level) и координат тайла (x, y):
Функция PostgreSQL, вычисляющая индекс [http://wiki.openstreetmap.org/wiki/QuadTiles quadkey] для заданного уровня (zoom level) и координат тайла (x, y):


<syntaxhighlight lang="plsql">
<syntaxhighlight lang="plsql">
Строка 341: Строка 338:




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


<syntaxhighlight lang="plsql">
<syntaxhighlight lang="plsql">
Строка 352: Строка 349:




[[Файл:Osmmp bing resolution map.jpeg|700|Разрешение спутниковых снимков Bing]]
[[Файл:Osmmp bing resolution map.jpeg|700px|Разрешение спутниковых снимков Bing]]


''Разрешение спутниковых снимков Bing''
''Разрешение спутниковых снимков Bing''
Строка 359: Строка 356:
== Создание индекса приоритетности ==
== Создание индекса приоритетности ==


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


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




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


Для того, чтобы использовать схему приоритетов в [http://qgis.org/ QGIS], добавьте столбец status в таблицу tile и создайте в плагине [http://docs.qgis.org/2.0/en/docs/training_manual/databases/db_manager.html DB Manager] динамически обновляемый слой TODO на основе SQL запроса:  
Для того, чтобы использовать схему приоритетов в [http://qgis.org/ QGIS], добавьте столбец ''status'' в таблицу ''tile'' и создайте в плагине [http://docs.qgis.org/2.0/en/docs/training_manual/databases/db_manager.html DB Manager] динамически обновляемый слой TODO на основе SQL-запроса:  


<syntaxhighlight lang="plsql">
<syntaxhighlight lang="plsql">
Строка 381: Строка 378:




[[Файл:Osmmp high priority tiles.jpeg|700|Тайлы задач с наивысшим приоритетом]]
[[Файл:Osmmp high priority tiles.jpeg|700px|Тайлы задач с наивысшим приоритетом]]


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




Перед тем, как начать наносить объекты в [http://wiki.openstreetmap.org/wiki/JOSM JOSM], либо сохраните выбранный тайл в GPX (отметьте опцию FORCE_GPX_TRACK и необходимость использования  [http://postgis.net/docs/ST_ExteriorRing.html ST_ExteriorRing]), или вызовите GDAL из [http://docs.qgis.org/2.8/en/docs/training_manual/create_vector_data/actions.html layer action]:
Перед тем, как начать наносить объекты в [http://wiki.openstreetmap.org/wiki/JOSM JOSM], либо сохраните выбранный тайл в файл формата GPX (отметьте опцию FORCE_GPX_TRACK и необходимость использования  [http://postgis.net/docs/ST_ExteriorRing.html ST_ExteriorRing]), или вызовите GDAL из [http://docs.qgis.org/2.8/en/docs/training_manual/create_vector_data/actions.html layer action]:


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
Строка 393: Строка 390:




Когда все интересующие Вас объекты, входящие в границы тайла задания, добавлены - обновите атрибут его статуса в базе данных (status = DONE) и обновите экран QGIS (нажав клавишу F5).
Когда все интересующие вас объекты, входящие в границы тайла задания, добавлены, обновите атрибут его статуса в базе данных (status = DONE) и обновите экран QGIS (нажав клавишу F5).
 
Используя данную схему приоритетов, автор замечал, что иногда [http://wiki.openstreetmap.org/wiki/DigitalGlobe снимки Mapbox] в более высоком разрешении, чем снимки Bing. Было бы полезно выяснить, предоставляют ли они схему разрешения снимков. Кроме того, две плитки последовательного ранга могут быть пространственно сильно удалены друг от друга, и было бы интересно найти способ их использования кластерами. Наконец, не следует ли включить приоритет непосредственно в Менеджер задач OSM?


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


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


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


== Источники ==


== Источник ==
Оригинал статьи: [http://www.mapaou-web.fr/2016/02/18/osm-mapping-priority/ OpenStreetMap Mapping Priority]


Автор статьи: [http://www.openstreetmap.org/user/adrienandrem adrienandrem]
Автор статьи: [http://www.openstreetmap.org/user/adrienandrem adrienandrem]


Автор перевода: [http://www.openstreetmap.org/user/Sibri Sibri]
Автор перевода: [http://www.openstreetmap.org/user/Sibri Sibri]

Текущая версия от 18:10, 17 июня 2017

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/osm-mapping-priority.html


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

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

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

Приоритет отрисовки в 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

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