Как трасформировать растр и вектор по точкам?

[Boris]

Прошу помочь найти описание процесса, если он конечно есть, как растр и вектор можно пересчитать при наличии парных точек и некоторых предположений о характере трасформации.

[gamm]

Прошу помочь найти описание процесса, если он конечно есть, как растр и вектор можно пересчитать при наличии парных точек и некоторых предположений о характере трасформации.

непонятно, о чем речь: нужно натянуть растр на вектор, вектор на растр, или и то и другое на глобус?

общия принцип прост: считаете преобразование по точкам полиномом или МВА (тем более, что вид известен, подбираете по вкусу), делаете сетку (как показывает опыт, шага 20-30 пикселей достаточно), и с нее снимаете новые координаты что растра, что вектора. Я на Питоне баловался, растр 10,000х10,000, трансформация десяток секунд. Для вектора берете с сетки новые координаты точек билинейной интерполяцией, для растра ей же узнаете, из какого места оригинала брать каждый пиксель (если нужно, интерполируете). Да, сетки для Х и Y свои.

[Boris]

Точки у меня есть - пары точек XY -> X1Y1, плоские координаты. Есть предположение какая именно трансформация должна быть - у GDAL их, если не путаю, 5 или 6. Можно как то GDAL/OGR натравить на вектор и растр, что он их пересчитал?

[Игорь Белов]

Растр:

1. gdal_translate с -gcp для опорных точек
2. gdalwarp

Вектор:

• ogr2ogr с -gcp для опорных точек

[Boris]

Спасибо!
Поскольку для ogr2ogr эта часть описания написана с ошибками - в заголовке команды

Код: Выделить всё

[-gcp pixel line easting northing [elevation]]*

a в описании параметров

Код: Выделить всё

-gcp ungeoref_x ungeoref_y georef_x georef_y elevation

пришлось поломать голову, что бы это значило.
Обидно, что в параметр нельзя передать имя файла. В параметры проекции можно, сюда нельзя. Пришлось думать как "запихнуть" много строк в параметры утилиты командной строки. И в завершение потребовался эксперимент с тем в каком направлении ставятся точки трансформации. Все пакеты делятся на два типа: первые сперва ставят "неверную точку", вторые - "верную(базовую)". OGR2OGR - имеет порядок "исходная"->"базовая".
И не понятно, что означает единственное поле "elevation".
Теперь по существу:
Имелся файл расстановки точек трансформации, в одном ПО, которое хорошо ставит точки, и хорошо дает предпросмотр результата, только трансформировать умеет один слой за раз, и только аффинное или гельмета преобразование, такого вида:

Спойлер

Код: Выделить всё

Name:
Vector UTM37 to MSK50z2 Registration 1

Description:
Vector KU 50v41 UTM37 to MSK50z2 Registration 1

    SourceSource       Target       Target
 #  Type    X: (m)        Y: (m)       X: (m)       Y: (m)       Residuals: (m)
 1  Control 466182,08     6156232,49   466180,79    6156230,78   0,00
...

RMS error: (m)                             0 

Очень не хотелось изобретать свое, на коленке сделанное. Потому и искал как сделать через GDAL/OGR. И вроде получилось. Слой пересчитался аналогично тому как делалось в экранном ПО. С растром буду разбираться позже.
Порядок точек совпадает В файле точек с порядком точек в ogr2ogr - "откуда берем"-> "куда должно попасть".
Для создания командного файла (Win32>=WinXP) для трансформации пришлось создать специальный BAT-файл:

Спойлер

Код: Выделить всё

@if "%~1"==""       goto iHelp
@if "%~1"=="-h"     goto iHelp
@if "%~1"=="/?"     goto iHelp
@if "%~1"=="--help" goto iHelp
@setlocal ENABLEEXTENSIONS ENABLEDELAYEDEXPANSION
@echo off
@rem echo on
set iTmp1=%Temp%\%~n0-%RANDOM%-%RANDOM%.tmp
set iTmp2=%Temp%\%~n0-%RANDOM%-%RANDOM%.tmp
@set path=%path%;c:\gnuwin32\bin;
@echo %DATE% %TIME% Do Init 1>&2
:MakeInit
set infile=%~f2
set outfile=%~f1
set gcpfile=%~f3
@echo %DATE% %TIME% Step 1 1>&2
echo ogr2ogr >"%iTmp1%"
if "%~4"=="" @echo -order 1 >>"%iTmp1%"
if not "%~4"=="" @echo -order %~4 >>"%iTmp1%"
if "%~5"=="" @echo -f "ESRI Shapefile" >>"%iTmp1%"

echo "%outfile%" >>"%iTmp1%"
echo "%infile%" >>"%iTmp1%"

@echo %DATE% %TIME% Step 2 1>&2
:LOOP01
@if "%~6"=="" goto ENDLOOP01

@echo "%~6" >>"%iTmp1%"
SHIFT /6
goto LOOP01
:ENDLOOP01
@echo %DATE% %TIME% Step 3 1>&2
gawk -F " " "{print \"-gcp \" $0;}" "%gcpfile%" >>"%iTmp1%"

@echo %DATE% %TIME% Step 4 1>&2
gawk  -F " " "{sub(/[ ]+$/,\"\",$0);sub(/^[ ]+/,\"\",$0);};(NR==1){t=$0;next;};{t=t \" \" $0;};END{print t;};" "%iTmp1%" >"%iTmp2%"
:ExitIt
@echo %DATE% %TIME% Exit... 1>&2
@rem type "%iTmp1%"
@rem echo -------------------------------------------
@type "%iTmp2%"
@rem @if exist "%iTmp1%" del /q/f "%iTmp1%"
@rem @if exist "%iTmp2%" del /q/f "%iTmp2%"
@echo Generated cmd file "%iTmp2%". >1&2
@endlocal
@exit /b 0
:iHelp
@echo Usage of: %~f0
@echo %~nx0 "Output OGR file" "Input OGR file" "GCP-list file" [Order] [OGR format] [OGR2OGR other argumets]
@echo Generated cmd file with multi OGR2OGR -gcp params. Gets -gcp from "GCP-list file".
@echo Arguments : 
@echo "Output OGR file"        Is String for OGR2OGR "dst_datasource_name"
@echo "Input OGR file"         Is String for OGR2OGR "src_datasource_name"
@echo "GCP-list file"          Is List of records for OGR2OGR "-gcp " option. Lines Are "ungeoref_x ungeoref_y georef_x georef_y [elevation]". No quotes. Space as delimeter.
@echo [Order]                  Is Digit for OGR2OGR "-order n" option. Value Is from 1 to 3, defaul = 1;
@echo [OGR format]             Is String for OGR2OGR "-f format_name" option. Value Is output file format name (default is ESRI Shapefile);
@echo [OGR2OGR other argumets] Is Digit for OGR2OGR "-order n" option. Value Is from 1 to 3, defaul = 1;

@exit /b -1

Кроме ogr2ogr >= 1.10.0 (не проверял, делал в 2.1.3), так же используется GAWK из GnuWin32.

[gamm]

Для создания командного файла (Win32>=WinXP) для трансформации пришлось создать специальный BAT-файл

чем городить такие турусы на колесах, не проще ли было затаскивать растры/вектора в R/Python (если они не по сотне Гб), и там спокойно всме делать. Если уж очень gdal для трансформации нужен, он там тоже привинчен, но и без него хватает разных методов. Размер программы будет меньше, ИМХО.

[Игорь Белов]

чем городить такие турусы на колесах, не проще ли было затаскивать растры/вектора в R/Python (если они не по сотне Гб), и там спокойно всме делать. Если уж очень gdal для трансформации нужен, он там тоже привинчен, но и без него хватает разных методов. Размер программы будет меньше, ИМХО.

Если поставить MSYS с той же OSGeo, турусы на колёсах превратятся в несколько строчек шелл-скрипта. Пусть, например, в папке лежат растры вместе с файлами точек привязки, сделанными в QGIS Georeferencer:

Код: Выделить всё

#!/bin/bash

for fpnt in *.points

do

    fimg=${fpnt%.*} ; f=${fimg%.*}

    cat ${fpnt} | awk 'BEGIN { FS = "," } (NR > 1 && $NF == 1) { printf " -gcp %.2f %.2f %.2f %.2f", $3-0.5, -$4-0.5, $1, $2 }' > tmp1.gcp

    xargs gdal_translate -of GTiff -a_srs EPSG:28406 -co COMPRESS=LZW ${fimg} tmp2.tif < tmp1.gcp

    gdalwarp -t_srs EPSG:3857 -dstnodata 0 -tps -co COMPRESS=LZW tmp2.tif ${f}_modified.tif

done

Файлы опорных точек других прогамм так же непринуждённо парсятся командой AWK.