Форумы GIS-Lab.info

Всем снова добрый день!
Задача состоит в следующем:
Даны результаты расчета двух моделей, каждая на поверхности Земли имеет регулярную сетку. Шаг сетки модели 2 меньше, чем у модели 1.
Требуется получить (проинтерполировать) значения модели 1 в узлы сетки модели 2. Конечно же с минимальными искажениями. Проекции у моделей тоже разные. Но даже если перепроецировать одну в другую, сетки будут лежать друг относительно друга как показано на рисунке.
Итак, значения в зеленых узлах нужно проинтерполировать в красные узлы. Существует ли способ лучше, чем перевод одной из сеток в вектор и использование инструментов интерполяции растра в векторные точки?
Изначально формат файлов GRIB, netCDF.

В горном софте есть такая процедура – «Опустить точки на поверхность». Построить по регулярной сетке MODEL1 поверхность. «Опустить точки» модели MODEL2 на эту поверхность, т.е. присвоить им значение атрибута Z.

Конечно, задача стандартная. По сути, вам нужно одно изображение трансформировать в пространственный домен другого. Пространственный домен это не только проекция, это:
1. Проекция - более общо скажем "система координат"
2. Границы сетки (xmin, ymin, xmax, ymax)
3. Размер ячейки (пространственное разрешение)

Так мы однозначно определяем (если опустить поворот, которого в 99% случаев нет) пространство нашей "сетки".

Решение на python для python gdal
https://github.com/silenteddie/RasterAdjuster

Просто через консольный gdal это делается тоже просто, gdal_warp:
gdalwarp -te <xmin ymin xmax ymax> -s_srs <проекция второго растра> -t_srs <проекция первого растра> -tr <размер ячейки по X, размер ячейки по Y> <второй растр> <трансформированный растр>

Таким образом приведем второй растр к первому.
xmin ymin xmax ymax и размеры ячейки должны совпадать с таковыми у первого растра. Можно посмотреть через gdalinfo.

Типом интерполяции можно управлять через ключ -r, там много, посмотрите здесь:
http://www.gdal.org/gdalwarp.html

если это модель, то лучше интерполировать с точек, поскольку интерполяция растра в растр - штука довольно специфическая и в данном случае не очень пригодная. Исходную сетку лучше перевести в XYZ, поменять проекцию, и использовать интерполятор на целевую сетку в SAGA, там методов хватает.

Да, я тоже думаю о конвертации, потому что у меня дало желаемый результат, но работало очень долго.
А при gdal_translate в xyz не будет потерь точности?

Вторая попытка сделать то же самое (будем называть это "алгоритм 1" - из растра в растр. Ошибки привожу в картинках.
Растр явно не записался, хотя привязка правильная, хоть и не совсем совпадают, как в первом тесте.
"алгоритм 2" - предложенный gamm будет опробован завтра, когда я пойму как записывать SET в cmd.

А попробуйте добавить в команду параметры:
-wo SAMPLE_GRID=YES -wo SAMPLE_STEPS=1000

Эдуард Казаков писал(а): ↑

18 июл 2018, 09:40

-wo SAMPLE_GRID=YES -wo SAMPLE_STEPS=1000

Привело к аналогичному результату.
https://lists.osgeo.org/pipermail/gdal- ... 37659.html тут кажется об этом.
Видимо, нужно обрезать большой растр сначала. вот только как получить маску обрезки.