Меню
Главная
Случайная статья
Настройки
|
Здесь описывается способ получения карты рельефа с помощью GMT, описанный de:User:Uwe Dedering вот здесь, с некоторыми дополнениями и вариациями.
- Инструменты
Ниже предполагается что всё это выполняется под Windows 7, хотя переносимости в принципе ничего не мешает.
Команды запускаются в папке C:\programs\GMT5\bin, а рабочие файлы все складываются в C:\programs\GMT5.
Вариант 1, на основе данных SRTM- 1. Получение базовой карты
Сначала определяемся с участком карты, для которого будем получать рельеф. Допустим, это будет позиционная карта, тогда у нас уже есть два диапазона, по широте и долготе, в градусах. Также есть размеры SVG либо PNG контурной карты, в пикселях — нужно перевести в дюймы по формуле: inches = pixels / 150.0 * 2.54. Координаты границ и размеры изображения подставляем в следующий командный файл:
cd bin
set GMTBIN=C:\programs\GMT5\bin
set WORKDIR=C:\programs\GMT5
set GSBIN=C:\PROGRA~1\gs\gs9.04/bin
set COORDSCUT=119.53/135.1/47.68/57.78
rem width = xmaxsvg / 150.0 * 2.54
set PAPERX=15.0706666667
rem height = ymaxsvg / 150.0 * 2.54
set PAPERY=16.2390666667
%GMTBIN%\gmtset GRID_PEN_PRIMARY = 0.1p,#404040
%GMTBIN%\pscoast.exe -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -Na -B1g1 -Ia/0.25p,#0978AB -W0.25,#0978AB -Df -P --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 > %WORKDIR%\Oblast_bounds.eps
if errorlevel 0 goto End
%GSBIN%\gswin32c.exe -dSAFER -dBATCH -dNOPAUSE -dGraphicsAlphaBits=4 -sDEVICE=pngalpha -dEPSCrop -r150 -sOutputFile=%WORKDIR%/Oblast_bounds.png %WORKDIR%\Oblast_bounds.eps
cd ..
В результате получаем изображение карты в равноугольной цилиндрической проекции (-JX), с градусной сеткой (-B1g1), изображением рек (-Ia) и водоёмов, государственных границ (-Na). Такая карта нужна для того чтобы проверить что указанные координаты правильные. Кроме того, по этой карте можно понять насколько адекватна имеющаяся контурная карта — совмещаем в Inkscape эти две карты друг с другом слоями с прозрачностью и смотрим что получилось.
- 2. Исходные данные
По изображению базовой карты мы видим сколько нам понадобится «квадратиков». Данные берём отсюда, каждый квадратик это отдельный файл с именем вида N42E044.hgt.zip. Скачиваем нужные файлы и разархивируем в C:\programs\GMT5.
Затем каждый квадратик переводим из hgt в grd:
xyz2grd ../N42E044.hgt -R44/45/42/43 -I3c -N-32768 -ZTLhw -G../n042e044.grd
- 3. Обработка данных
Для дальнейших вычислений объединяем квадратики в большие площади. Объёмы данных больше 7-9 квадратиков могут не поместиться в память при дальнейших вычислениях (зависит от объёма памяти, доступного для программ). Обычно я объединяю квадратики в «горизонтальные» блоки и дальше работаю уже с ними:
grdpaste ../n042e044.grd ../n042e045.grd -G../temp.grd
grdpaste ../temp.grd ../n042e046.grd -G../n042.grd
Переводим данные блоков в текстовый формат.
grd2xyz ../n042.grd > ../n042.txt
Теперь используем утилиту surface, которая выполняет сплайновую интерполяцию для точек, пропущенных в исходных данных. На этом этапе происходит наибольший расход оперативной памяти, поэтому мы и обрабатываем не всю карту целиком, а по блокам. Поскольку на границе блоков данные не интерполируются, то теоретически это может привести к тому что границы блоков станут различимы. На практике этого не наблюдается, видимо за счёт того что плотность точек на входе гораздо выше пиксельной сетки на выходе.
surface ../n042.txt -R44/47/42/43 -I3c -T0.35 -G../n042_cor.grd
Полученные данные нужно обрезать по границам интересующей нас области.
grdcut.exe ../n042_cor.grd -R44.8/46.7/42.4/44.1 -G../n042_land.grd
Получаем shaded relief с помощью grdgradient.
grdgradient.exe ../n042_land.grd -A315 -Ne0.2 -G../n042_land_schatten.grd
Дальше уже можем объеденить полученные файлы в один общий блок:
grdpaste.exe ../n042_land.grd ../n043_land.grd -G../Oblast_land.grd
grdpaste.exe ../n042_land_schatten.grd ../n043_land_schatten.grd -G../Oblast_land_schatten.grd
В результате у нас получились файлы Oblast_land.grd и Oblast_land_schatten.grd, которые используются при рендеринге.
- 4. Рендеринг
Для рендеринга используются утилиты pscoast (даёт обрезку по береговым линиям) и grdimage (рендерит рельеф), обе они выдают результат в PostScript. Для правильной раскраски используется файл wiki-land-verlauf2.cpt с содержимым отсюда.
%GMTBIN%\pscoast.exe -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -Gc -Df -P -K --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 > %WORKDIR%\Oblast_land.eps
%GMTBIN%\grdimage.exe ../Oblast_land.grd -I../Oblast_land_schatten.grd -C%WORKDIR%\wiki-land-verlauf2.cpt -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -P -K -O --PAPER_MEDIA=Custom_14.0cx14.0c -X0 -Y0 >> %WORKDIR%/Oblast_land.eps
%GMTBIN%\pscoast.exe -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -Q -Df -P -O -K >> %WORKDIR%\Oblast_land.eps
%GMTBIN%\pscoast.exe -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -O -W0.25,#0978AB -Df -P --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 >> %WORKDIR%\Oblast_land.eps
Здесь ключ -W задаёт рисование береговых линий.
Наконец, для получения PNG из PostScript используется GS:
%GSBIN%\gswin32c.exe -dSAFER -dBATCH -dNOPAUSE -dGraphicsAlphaBits=4 -sDEVICE=pngalpha -dEPSCrop -r150 -sOutputFile=%WORKDIR%/Oblast_relief.png %WORKDIR%\Oblast_land.eps
При использовании нескольких вызовов утилит GMT, отдающих результат в виде PostScript, нужно указывать ключи: -K если не нужно завершать PS-файл (будет продолжение), -O если не нужно генерировать начало PS-файла (если оно уже сгенерировано до этого).
Вариант 2, на основе данных ETOPO
Исходные данные: ETOPO1 Global Relief Model — сетка рельефа и батиметрии с разрешением в 1 угловую минуту. Хорошо подходит для карт размером от нескольких градусов по каждой стороне.
Пример создания карты Курильских островов:
set PATH=C:\programs\GMT5\bin;%PATH%
set GSBIN=C:\PROGRA~1\gs\gs9.04/bin
set COORDSCUT=144.7492/157.3007/42.9694/51.3837
rem width = xmaxsvg / 150.0 * 2.54
set PAPERX=14.9352
rem height = ymaxsvg / 150.0 * 2.54
set PAPERY=16.0189333333333
grdcut.exe ETOPO1_Bed_g_gmt4.grd -R%COORDSCUT% -Gh_cor_cut.grd
grdgradient h_cor_cut.grd -Ne0.3 -A315 -M -Ghi.grd
if not errorlevel 0 goto End
grdimage h_cor_cut.grd -Ihi.grd -Cwiki-water-verlauf2.cpt -P -R%COORDSCUT% -JX%PAPERX%cd/%PAPERY%cd --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 -K > map.eps
if not errorlevel 0 goto End
pscoast.exe -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -Gc -P -Df --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 -O -K >> map.eps
if not errorlevel 0 goto End
grdimage h_cor_cut.grd -Ihi.grd -Cmount.cpt -P -R%COORDSCUT% -JX%PAPERX%cd/%PAPERY%cd --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 -O -K >> map.eps
if not errorlevel 0 goto End
pscoast.exe -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -Q -P -Df --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 -O -K >> map.eps
if not errorlevel 0 goto End
pscoast.exe -JX%PAPERX%cd/%PAPERY%cd -R%COORDSCUT% -Na -Ia/0.25p,#0978AB -W0.25,#0978AB -P -Df --PAPER_MEDIA=Custom_%PAPERX%cx%PAPERY%c -X0 -Y0 -O >> map.eps
if not errorlevel 0 goto End
%GSBIN%\gswin32c.exe -dSAFER -dBATCH -dNOPAUSE -dGraphicsAlphaBits=4 -sDEVICE=pngalpha -dEPSCrop -r150 -sOutputFile=Oblast_etopo.png map.eps
if not errorlevel 0 goto End
del h.grd
del hi.grd
:End
Постобработка
После того как получено изображение рельефа, остаётся использовать его для создания готовой карты — обычно для этого нужно только наложить сверху контурную карту области. Это можно сделать с помощью слоёв в Inkscape.
Если на контурной карте нет рек и водоёмов, можно получить их через GMT — см. выше пример с получением базовой карты, ключи: -Ia, -W, -S. Реки и водоёмы должны лежать выше всего остального, поэтому лучше сделать с ними отдельное изображение, которое потом наложить отдельным слоем в Inkscape.
|
|
