Форумы GIS-Lab.info

Добрый день. Помогите, пожалуйста, понять, в чём проблема.
Нужны растровые слои, где в значениях каналов будут расстояния в метрах до тех или иных объектов. Сначала до дороги. Основываюсь на алгоритме, описанном здесь:
http://www.mikemeredith.net/blog/1212_G ... n_QGIS.htm
Проект с СК WGS84/UTM54N.
ИЗ OSM беру слой дорог, обрезаю лишнее, пересохраняю слой в СК проекта - road.
Добавляю слой plots, чтобы ориентироваться на исследуемую площадь.
Растеризую слой road, вручную указывая границы искомого слоя - получаю road_raster.tif о значениями 1 по месту нахождения дороги, остальные пиксели со значением "0".
Далее иду Растр -> Анализ -> Карта близости. Исходный файл road_raster, целевой road_dist. Единицы измерения - GEO. Но прога выдаёт ошибку:

ERROR 6: GDALDriver::Create() ... no create method implemented for this format.
Traceback (most recent call last):
File "C:\OSGeo4W\bin\gdal_proximity.py", line 166, in
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
AttributeError: 'NoneType' object has no attribute 'SetGeoTransform'

Если убрать параметр "-of GMT", то процесс проходит, но снова ошибка:

Warning 1: Pixels not square, distances will be inaccurate.

Но растр road_dist добавляется. Как и положено, на слое дорог значения "0", увеличиваются по мере удаления.
На первый взгляд, всё адекватно, но на проверку оказывается, что значения сильно разнятся с фактическими.
Точки, которые линейкой в том же QGISе или SAS.Планете имеют одни значения удалённости от дороги, имеют другие значения, если ориентироваться на значения пикселя слоя road_dist. Наиболее близкие точки имеют разницу в 10-15 м, наиболее отдалённые - в 2000 м (фактические и "линеечные" значения меньше, чем на полученном растре).
Понимаю, что проблема, скорее всего, кроется в выборе проекций/систем координат. Но пробовал всё то же самое с World Equidistant Cylindrical, WGS84/Plate Carree, WGS84/World Mercator, Pseudo Mercator - бесполезно.
Поэтому очень прошу подсказать, в чём моя ошибка, или что дополнительно нужно указать для решения проблемы?
Спасибо!

Использовать нужно проекцию UTM, остальные из перечисленных сильно искажают расстояния.
Приведите вывод gdalinfo для исходного растра.

Значит, WGS 84 / UTM 54N должна подходить.

Вот вывод gdalinfo по итоговому растру с расстояниями (+ txt в приложении):

W>gdalinfo c:\road_dist.tif
Driver: GTiff/GeoTIFF
Files: c:\road_dist.tif
c:\road_dist.tif.aux.xml
Size is 3000, 3000
Coordinate System is:
PROJCS["WGS 84 / UTM zone 54N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",141],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32654"]]
Origin = (642675.000000000000000,5173418.000000000000000)
Pixel Size = (4.375000000000000,-3.301333333333334)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 642675.000, 5173418.000) (142d51'58.46"E, 46d41'56.75"N)
Lower Left ( 642675.000, 5163514.000) (142d51'47.43"E, 46d36'36.05"N)
Upper Right ( 655800.000, 5173418.000) (143d 2'16.12"E, 46d41'46.21"N)
Lower Right ( 655800.000, 5163514.000) (143d 2' 4.09"E, 46d36'25.54"N)
Center ( 649237.500, 5168466.000) (142d57' 1.53"E, 46d39'11.25"N)
Band 1 Block=3000x1 Type=Float32, ColorInterp=Gray
Min=0.000 Max=12542.172
Minimum=0.000, Maximum=12542.172, Mean=4546.033, StdDev=3178.091
Metadata:
STATISTICS_MAXIMUM=12542.171875
STATISTICS_MEAN=4546.0333195243
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=3178.0908246734

Tifoso писал(а): Pixel Size = (4.375000000000000,-3.301333333333334)

Тут надо разбираться, почему на этапе растеризации у вас получились не квадратные пиксели. Судя по предупреждению GDAL, алгоритм на такой случай не рассчитан. Можно также, если вы уверены, что с привязкой растра всё в порядке, сделать ресемплинг растра, чтобы у него было одинаковое разрешение по горизонтали и вертикали. Вроде бы в QGIS при сохранении растра можно было задать разрешение. Ну или в консоли gdalwarp с параметром -tr

Ariki писал(а):на этапе растеризации у вас получились не квадратные пиксели

Проблема прояснилась. Об этом предупреждал модуль Карта близости сообщением "Warning 1: Pixels not square, distances will be inaccurate".
Моя ошибка состояла в том, что при растеризации вектора, указав координаты крайних точек, я также задал параметр -ts вместо -tr, тем самым установил размер получаемого файла в пикселях вместо того, чтобы задать разрешение пикселя на единицу измерения СК (если я правильно уловил суть).
Большое спасибо, что натолкнули на решение!
Теперь все данные вполне адекватны.

Но всё равно, пока не удалю параметр -of GMT перед запуском Карты близости, так и выскакивает Error 6 (текст в первом посте).

А зачем вам GMT? Параметр -of задаёт выходной формат, в данном случае - какой-то экзотический вариант NetCDF, который, судя по всему, не поддерживается GDAL на запись (ну или в вашу сборку QGIS не включена поддержка). По умолчанию же создаётся растр в формате GeoTIFF.