ВИЗУАЛИЗИРАНЕ НА ОКЕАНОГРАФСКИ ДАННИ ЗА ЧЕРНО МОРЕ С MATLAB

Автор: Тодор Йорданов

Email: todor.iordanov@hippocampus-bg.org

Океанографските данни за Черно море са в различни формати (бинарни или ASCII) и всеки формат има различна структура. Например данните могат да са във вид на колонки:

Longitude       Latitude          Quantity (e.g. wind speed)”
30                    52                    22.8
30.5                 53                    22.9
.                       .                       .
.                       .                       .
.                       .                       .

Или във вид на матрица:

Longitude ———————————————————–>
L
A         22.8     23.9     27.1     …
T          21.4     24.6     30.9     …
I          .           .
T          .                       .
U         .                                  .
D
E
|
|
|
|
V

В това ръководство се разглеждат няколко различни начини за визуализиране на океанографски данни отнасящи се за Черно море. Естествено, изброените методи могат да бъдат приложени и за океанографски данни от други морета.

Удобен вариант за визуализиране на океанографски данни е използването на скриптове в Matlab. Предимството на Matlab е, че разполага с голям запас от функции за визуализиране и четене на различни формати данни. Недостатък е, че е платена програма и струва не малко пари.

Не е възможно да се разгледат всички видове данни и формати. Затова в това ръководство два вида данни са избрани за пример и различните скриптове са демонстрирани на базата на тях. Единият пример е на базата на данни за батиметрията на Черно море от проекта ETOPO (National Geophysical Data Center, National Oceanic and Atmospheric Administration), а другият е на базата на данни за повърхностен вятър от The Physical Oceanography Distributed Active Archive Center (PO.DAAC) при NASA.

Батиметрия

Първият пример е на базата на батиметрични данни във формат XYZ:

26.000000  48.000000    314
26.033333  48.000000    306
26.066667  48.000000    345
26.100000  48.000000    372
26.133333  48.000000    378
26.166667  48.000000    383
26.200000  48.000000    389
26.233333  48.000000    349
26.266667  48.000000    253
26.300000  48.000000    230
26.333333  48.000000    211
26.366667  48.000000    200
.
.
.

Да предположим, че данните се намират във файл с името BlackSea.xyz. В такъв случай командата на Matlab за прочитане на файла е:

%% Load data

load BlackSea.xyz

След като данните бяха прочетени се появява нова структура в Matlab със същото име както файла: BlackSea. Следващата ни цел е да преобразуваме данните от XYZ-формат във вид на матрица. Това може да се осъществи със следната функция:

% Transform the xyz data in matrix data

[bs_data, xgrid, ygrid] = xyz2Matrix(BlackSea);

Резултатът от тази функция е записан в структурите: bs_data, xgrid и ygrid. bs_data съдържа батиметричните данни, а xgrid и ygrid съдържат съответните координати. Допълнителна (и същевременно опционална) стъпка за визуализирането на данните е прочитането на подходяща цветова таблица. Предварително създадена цветова таблица може да бъде свалена тук. Ако искате да създадете собствена цветова таблица може да използвате или вградената в Matlab помощна програма colormapeditor или някоя от свободно достъпните програми (например http://www.mathworks.com/matlabcentral/fileexchange/5043-colormapsnc). Новосъздадената цветова таблица трябва да бъде запаметена за да може да се използва и в бъдеще. Описание за това как създадената цветова таблица може да бъде запаметена можете да намерите на страницата на Matlab централата (http://www.mathworks.com/matlabcentral/newsreader/view_thread/97964). За да продължим с визуализирането трябва да заредим предварително запаметената цветова таблица. В този случай това става така:

% Load colormap

load ‘topocolormap.mat’

Вече можем да видим и първите резултати:

Black Sea bathymetryТази графика може да бъде създаде със следния скрипт:

% Contour plot of the data

figure;

contourf(xgrid, ygrid, bs_data, 321, ‘LineStyle’, ‘none’);

colormap(topocolormap)

xlim([26 43])

caxis([-4000 4000])

hold on

[C, h] = contour(xgrid, ygrid, bs_data, ‘-k’, ‘LineWidth’, 2);

colormap(topocolormap)

text_handle = clabel(C,h);

set(text_handle, ‘FontSize’, 14);

set(h,’ShowText’,’on’,’TextStep’,get(h,’LevelStep’)*1)

xlim([26 43])

caxis([-4000 4000])

hold off

Можете спокойно да експериментирате с параметрите на функциите за визуализиране и да създадете графика пригодена на личните нужди. В основата на това визуализиране е функцията contourf. Опростена употреба на тази функция може да изглежда така:

figure;

contourf(xgrid, ygrid, bs_data);

xlim([26 43])

caxis([-4000 4000])

colorbar

Резултатът от този код е следната графика:

Black Sea bathymetryТази графика е опростена в сравнение с предишната и липсва информация (например надморската височина).

Втори метод за визуализиране на данните е с функцията surf. Един примерен резултат от тази функция може да изглежда така:

Black Sea bathymetryФункцията surf оперира в три-димензионалното пространство. При нея могат да се използват параметри за осветление и сянка, с помощта на които 3D ефектът се засилва (както е видно от горната графика). Кодът за генериране на графиката е следният:

% Plot the barhymetry

h = surf(xgrid, ygrid, bs_data, ‘LineStyle’, ‘none’);

%shading interp

%lightangle(0, 30)

material dull

set(h,’FaceLighting’,’phong’,’FaceColor’,’interp’,…

‘AmbientStrength’,0.4)

light(‘Position’,[1 0 0],’Style’,’infinite’);

view(0, 90)

xlim([26 43])

caxis([-4000 4000])

colormap(topocolormap)

colorbar

hold on

[C, h1] = contour3(xgrid, ygrid, bs_data);

set(h1,’LineWidth’, 2)

hold off

Следваща стъпка, която бихме искали да изпълним е да запаметим графиката във формата на файл. Добро качество на запаметената графика може да се постигне ако се използва форматът PDF. Следните команди ориентират и запаметяват графиката:

orient landscape

print -dpdf -r300 bs_bathymetry.pdf

Файлът bs_bathymetry.pdf може да бъде едитиран с Photoshop например и трансформиран в png-графика.

Повърхностен вятър

За вторият пример ще използваме данни за повърхностният вятър над Черно море. Допълнителна трудност при тези данни е, че са в HDF-формат. За щастие Matlab разполага с функции за четене на този формат данни. Прочитането на данните се извършва със следните команди:

% Read HDF surface wind data

des_avg_wind_vel_v = double(hdfread(‘QS_XWGRD3_2009179.hdf’, ‘des_avg_wind_vel_v’));

des_avg_wind_vel_u = double(hdfread(‘QS_XWGRD3_2009179.hdf’, ‘des_avg_wind_vel_u’));

След изпълнението на тези команди ще се появят две нови матрици в работното прострнаство на Matlab: des_avg_wind_vel_u и des_avg_wind_vel_v.

Тъй като тези данни не са само за Черно море, а за целия Свтовен океан, първо е необходимо да екстрахираме само тези данни, които се отнасят за Черно море. Това се осъществявя със следния скрипт:

% Generate axis data

x = 0:0.25:359.75;

y = -89.75:0.25:90;

z = sqrt((des_avg_wind_vel_u.^2 + des_avg_wind_vel_v.^2));

 

% Extract Black sea data

start_x = 26;

end_x = 43;

start_y = 40;

end_y = 48;

step_x = x(2)-x(1);

step_y = y(2)-y(1);

start_x_idx = find(x > start_x – step_x & x < start_x + step_x);

end_x_idx = find(x > end_x – step_x & x < end_x + step_x);

start_y_idx = find(y > start_y – step_y & y < start_y + step_y);

end_y_idx = find(y > end_y – step_y & y < end_y + step_y);

bs_wind_speed = z(start_y_idx:end_y_idx, start_x_idx:end_x_idx);

bs_wind_vel_u = des_avg_wind_vel_u(start_y_idx:end_y_idx, start_x_idx:end_x_idx);

bs_wind_vel_v = des_avg_wind_vel_v(start_y_idx:end_y_idx, start_x_idx:end_x_idx);

x1 = x(start_x_idx:end_x_idx);

y1 = y(start_y_idx:end_y_idx);

Екстрахираните данни можем да визуализираме със следните команди:

% Plot wind speed and direction

n_isolines = 100; % number of isolines for the contour plot

arrow_length = 2;

arrow_color = ‘-w’; % w=white, r=red, b=blue, g=green

arrow_thickness = 1.5;

 

figure; contourf(x1, y1, bs_wind_speed, n_isolines, ‘LineStyle’, ‘none’)

hold on

[mx, my] = meshgrid(x1, y1);

quiver(mx, my, bs_wind_vel_u, bs_wind_vel_v, arrow_length, arrow_color, ‘LineWidth’, arrow_thickness);

Резултатът от този скрипт е показан на следващата графика:

Black Sea surface windНедостатъкът на тези данни е, че няма стойности за всяка точка от матрицата. Това е видно на горната графика. Примерно решение на проблема е да се интерполират данните които липсват. Преди да интерполираме данните трябва да ги подготвим за интерполация като първо ги трансформираме от матричен формат в XYZ-формат и после изтрием всички точки, които нямат валидни стойности (в този случай това са всички точки със стойност 0).

% Transform the matrix data to xyz data

bs_wind_speed_xyz = matrix2XYZ(bs_wind_speed, x1, y1);

 

% Remove the zero data points

tmp_xyz = removeZeros(bs_wind_speed_xyz);

% Interpolate the data

[xi, yi] = meshgrid(x1, y1);

zi = griddata(tmp_xyz(:, 1), tmp_xyz(:, 2), tmp_xyz(:, 3), xi, yi, ‘cubic’);

% Plot the interpolated data

n_isolines = 100; % number of isolines for the contour plot

arrow_length = 2;

arrow_color = ‘-k’; % w=white, r=red, b=blue, g=green, k=black

arrow_thickness = 1.5;

figure;

contourf(xi, yi, zi, n_isolines, ‘LineStyle’, ‘none’);

xlabel(‘Longitude’), ylabel(‘Latitude’), zlabel(‘Depth in Feet’)

hold on

quiver(xi, yi, bs_wind_vel_u, bs_wind_vel_v, arrow_length, arrow_color, ‘LineWidth’, arrow_thickness);

Резултатът от горните операции е показан на следващата графика:

Black Sea surface windИнтерполацията запълни липсващите данни но се загуби формата на Черно море. За да решим този проблем можем да интерполираме данните още веднъж като използваме същият грид както при батиметричните данни и после използваме батиметричните данни за да маскираме частта от повърхностния вятър, която не е над водната повърхност.

% Interpolate the wind data on the grid of the bathymetry data

XI = xgrid;

YI = ygrid;

ZI = griddata(tmp_xyz(:, 1), tmp_xyz(:, 2), tmp_xyz(:, 3), XI, YI, ‘v4’);

 

% Create masked data

masked_z = createDataMask(ZI, bs_data);

 

% Plot the interpolated masked data

n_isolines = 100; % number of isolines for the contour plot

arrow_length = 1.5;

arrow_color = ‘-k’; % w=white, r=red, b=blue, g=green, k=black

arrow_thickness = 1.5;

figure;

contourf(XI, YI, masked_z, n_isolines, ‘LineStyle’, ‘none’)

caxis([0 700])

colorbar

hold on

quiver(xi, yi, bs_wind_vel_u, bs_wind_vel_v, arrow_length, arrow_color, ‘LineWidth’, arrow_thickness);

Резултатът от горните операции е показан на следващата графика:

Black Sea surface windЕстествено и тази графика може да бъде запаметена като PDF-файл както беше описано в частта за батиметричните данни.

С това завършваме краткото ръководство за визуализиране на океанографски данни с Matlab. Не всички възможности на Matlab са разгледани, но това би отнело твърде много време. Надявам се въпреки това то да помогне на начинаещите океанографи.

Leave a Reply

Your email address will not be published. Required fields are marked *