Как построить цветную топографическую карту в MATLAB

Сегодня мы построим топографическую карту (карту рельефа) для всего мира, а также для отдельных районов, используя средства MATLAB. Нам как океанологам, конечно, важна её батиметрическая часть, но строить мы будем и "сухопутную" часть тоже. Сегодня ограничимся только картой, где топография отлита цветом. Изолинии будем рисовать в другой раз.

Начнём с того, что используем родной файл с топографией, который встроен в MATLAB. 

load topo

В результате этой команды откроются несколько переменных, главная из которых topo, в которой и содержатся нужные нам значения высот

whos
  Name              Size              Bytes  Class     
  topo            180x360            518400  double              
  topolatlim        1x2                  16  double              
  topolegend        1x3                  24  double              
  topolonlim        1x2                  16  double              
  topomap1         64x3                1536  double              
  topomap2        128x3                3072  double              

В прошлый раз (советуем прочитать тот материал, прежде чем читать дальше) для построения карты мы использовали матрицы с широтами и долготам, равные по размерам матрице с данными. Здесь же нет даже векторов, из которых мы эти матрицы можем сделать. Но на самом деле, при равномерной координатной сетке мы можем обойтись и без них. Поэтому сейчас мы построим карту немного другим способом (отличным от того, которым пользовались ранее).

Если мы знаем размер матрицы с топографией и знаем пределы этой карты (по широте и долготе), то можно сказать, что наша координатная сетка нам известна. Функция georefcells по пределам широт и долгот, а также размеру матрицы topo создаст в матлабе референсную ячейку, которую впоследствии мы сможем использовать для построения:

topoR = georefcells(topolatlim,topolonlim,size(topo))

Теперь осталось дело за малым -  построить карту:

worldmap('world')
geoshow(topo,topoR,'DisplayType','texturemap')
топографическая карта MATLAB

Для того, чтобы цвета выглядели более привычно для географа, лучше воспользоваться функцией demcmap, которая создаст специальную цветовую политру (colormap), подходящую для отображения рельефа.

worldmap('world')
geoshow(topo,topoR,'DisplayType','texturemap')
demcmap(topo)
Цветовая палитра для топографической карты

Что ж, карта готова.  Но что, если мы хотим построить карту не для всего Мирового океана, а лишь для небольшого района?

worldmap([50 80],[-20 60])
geoshow(topo,topoR,'DisplayType','texturemap')
demcmap(topo)
Топография Европы

Как видно градусное разрешение данных topo плохо подходит для такого рода построений, а если мы построим подобную карту для отдельного моря, то совсем не увидим никакой полезной информации. Поэтому нам нужны данные с более детальным рельефом и на помощь нам придёт топографический массив данных ETOPO. Мы скачаем данные с минутным разрешением (ETOPO1) по этой ссылке. Поизучайте папку на сайте, в ней вы может найти красивые картинки с батиметрией, плакаты, файл с документацией по ETOPO, а также данные в разных форматах. Мы будем использовать файл etopo1_ice_c_f4.zip в папке data/ice_surface/cell_registered/binary/. Скачиваем и распаковываем.

Для работы с бинарным файлом ETOPO в MATLAB существует специальная функция etopo. У этой функции на выходе две переменной а на входе может быть использовано одно лишь название вашего файла:

[Z, refvec] = etopo('etopo1_ice_c_f4.flt');

В результате этого действия в рабочей области матлаба создастся достаточно объёмная матрица Z размером 10800x21600 и референсный вектор, состоящий из трёх значений. Первое значение определяет размер ячеек в долях градуса (1/refvec(1)), в нашем случае оно равно 60, т.е. дискретность сетки - 1/60 градуса. Второе значение вектора содержит значение широты самого северозападного угла, а третье - долготу того же самого угла. Как мы увидим далее, использование этого вектора - ещё один способ привязать вашу матрицу к координатной сетке при использовании функции geoshow. Сама же матрица Z содержит только значения нашей топографии. 

Поскольку нам не нужна такая большая матрица, то для ускорения процессов вырежем только нужный нам кусок. Делается это добавлением соответсвующих аргументов в функцию etopo:

[Z, refvec] = etopo('etopo1_ice_c_f4.flt',1, [50 80],[-20 60]);

Аргумент, который у нас равен '1' - это samplefactor, если он равен единице, то значит мы на выходи получим исходное разрешение. Если нам не нужна большая подробность, то можно это число увеличить. Отсальные аргументы -  векторы с пределами карты. 

Теперь строим карту следующим образом:

worldmap([50 80],[-20 60])
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)
Топографическая карта Европы в MATLAB

Обратите внимание, что если мы при построении захотим быстро уменьшить район (к примеру, Белое море сюда полностью попадает и мы хотим просто "подрезать" карту со всех сторон, изменив лишь аргументы в worldmap):

worldmap([63.8 67.1],[32 45])
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap(Z)

то может получиться не очень хороший вариант. Например пропали градации по глубине:

Батиметрия Белого моря

Произошло это из-за того, что цветовая палитра была построена для диапазона глубин из исходного большого района. Поэтому нам нужно либо ещё раз запустить функцию etopo с новыми координатами, либо перезапустить demcmap с указанием диапазона значений (в метрах), в пределах которых изменяется топография в нашем районе. Мы выберем второй вариант:

worldmap([63.8 67.1],[32 45])
geoshow(Z, refvec, 'DisplayType', 'texturemap');
demcmap([-300 600])
Построение батиметрии моря в MATLAB

Таким образом, мы можем более тонко настраивать нашу цветовую шкалу.

На получившейся карте слишком мало градаций света, увеличить их число можно при помощи дополнительного параметра в demcmap, также добавим шкалу цветов на карту.

demcmap([-300 600],150)
colorbar
colorbar для батиметрии и топографической карты MATLAB

До этого момента мы использовали для построения тип 'texturemap', однако можно использовать и 'surface', который пригодится, если вы захотите изменить угол просмотра (функция view) вашей карты (иными словами, сделать её трёхмерной). Например, для нашей карты сделаем азимут 10 и высоту 35 (view(10,35)). Однако одного view будет недостаточно. Для того, чтобы эффект трёхмерности был заметен нужно увеличить масштам по оси z. Делается, это при помощи daspectm, поиграйтесь с разными значениями, мы остановились на 200:

worldmap([63.8 67.1],[32 45])
geoshow(Z, refvec, 'DisplayType', 'surface');
demcmap([-300 600])
daspectm('m',200)
view(10,35)
Трёхмерная карта рельефа

Можем добавить эффект освещения, запустив функцию

camlight
Трёхмерный рельеф Белого моря с батиметрией и топографией в MATLAB

Давайте посмотрим, как будет выглядеть карта с эффектом освещения, но только со стандартным азимутом (90):

worldmap([63.8 67.1],[32 45])
geoshow(Z, refvec, 'DisplayType', 'surface');
demcmap([-300 600])
daspectm('m',200)
camlight
colorbar
Эффект освещения на карте в MATLAB

Если вам не хочется скачивать полный архив ETOPO который занимает почти 1 гб, то на сайте есть сервис для вырезания данных для отдельного района: https://maps.ngdc.noaa.gov/viewers/wcs-client/

Скачать можно в разных форматах, давайте выберем знакомый нам netCDF формат и попробуем построить батиметрическую карту, но используя способ, подобный тому, который мы использовали в прошлой заметке. Файл скачался под названием etopo1.nc. Долго объяснять не будем, приведём код и результат (всё в общем то похоже, но в качестве готового рецепта решили привести код здесь)

Z = ncread('etopo1.nc','Band1');
lat = ncread('etopo1.nc','lat');
lon = ncread('etopo1.nc','lon');
[lon,lat]=ndgrid(lon,lat);

worldmap([63.8 67.1],[32 45])
geoshow(lat,lon, Z, 'DisplayType', 'texturemap');
demcmap([-300 600],120)
colorbar
Батиметрия Белого моря

Ну что ж, кажется у вас теперь достаточно навыков, чтобы построить топографическую карту не только по данным ETOPO, но и по своим собственным.  

 telegram

Чтобы не пропустить новые материалы, подпишитесь на канал в Telegram: https://t.me/koldunovaleksey