新浪博客

GMT/MatlabToolbox的高级举例

2023-03-31 10:05阅读：

http://blog.sina.cn/dpool/blog/u/2798219670


GMT/Matlab Toolbox的高级举例

GMT 的 Matlab 接口提供了在Matlab 中调用 GMT 命令的功能。通过该接口，GMT 的所有模块命令都可以在 Matlab 脚本中嵌入执行。GMT 命令生成的结果（grid 格网数据、table 表格数据、CPT 颜色表、文本文件、图片等）都可以作为 Matlab 变量进行运算；Matlab 中的矩阵变量也可以直接作为GMT的输入。
如何设置调用路径和一些简单的例子，可参考GMT中文手册中GMT/Matlab Toolbox相关内容（https://docs.gmt-china.org/6.4/api/matlab/）。
本例子最初是由pygmt官方给出(https://www.pygmt.org/latest/gallery/lines/vector_styles.html#sphx-glr-gallery-lines-vector-styles-py)。
我把它转化成了MATLAB版本，主要遇到的问题是pstext需要传递的是一个多种数据类型的元胞数组（cell array）。Wessel and Luis （2017）的文章（Wessel, P., and J. F. Luis (2017), The GMT/MATLAB Toolbox, Geochem.Geophys. Geosyst., 18, 811–823,doi:10.1002/2016GC006723.）对此说明并不是很清楚，也没有例子。于是我在github上提了问题（https://github.com/GenericMappingTools/gmtmex/issues/32）得到了开发者Wessel的回复，才得以解决。

具体代码如下

outfn = 'ex2HL_arrows'; % must single quotaion

%% basemap
gmt(['pscoast -R70/150/0/60 -JM15c -Ba -A4000 -N1 -W0.25p,black -K > ' outfn '.ps'])

%% plot 12 Cartesian vectors with different lengths
x = linspace(109, 109, 12)';
y = linspace(45, 55, 12)';
direction = zeros(size(x));
length = linspace(0.5, 1.5, 12)';
% Cartesian vectors (v) with red pen and fill (+g, +p), vector head at
% end (+e), and 40 degree angle (+a) with no indentation for vector head (+h)
style = 'v0.2c+e+a40+gred+h0+p1p,red';
vector = [x y direction length];
% gmt('psxy', ['-R20/60/-180/0 -JM12c -W1p,red -S' style ' > ' outfn '.ps'], vector);
gmt(['psxy -R -J -Ba -W1p,red -S' style ' -K -O >> ' outfn '.ps'], vector);
txt = gmt('record', [113, 43], 'CARTESIAN');
gmt(['pstext -R -J -F+f13p,Helvetica-Bold,red -K -O >> ' outfn '.ps'], txt);

%% plot 7 math angle arcs with different radii
num = 7;
x = linspace(95, 95, num)'; % x coordinates of the center
y = linspace(33, 33, num)'; % y coordinates of the center
radius = 1.8 - 0.2 * (0 : num - 1)';% radius
startdir = linspace(90, 90, num)'; % start direction in degrees
stopdir = 180 + 40 * (0 : num - 1)'; % stop direction in degrees
%# data for circular vectors
data = [x y radius startdir stopdir];
arcstyle = 'm0.5c+ea'; % Circular vector (m) with an arrow at end
gmt(['psxy -R -J -Gred3 -W1.5p,black -S' arcstyle ' -K -O >> ' outfn '.ps'], data)
txt = gmt('record', [95, 42], 'CIRCULAR');
gmt(['pstext -R -J -F+f13p,Helvetica-Bold -K -O >> ' outfn '.ps'], txt);

%% plot geographic vectors using endpoints
PEK = [116.41, 39.9]; % Peking
XAN = [108.9, 34.3]; % Xi'an
XJG = [87.62, 43.79]; % XinJiang
KMG = [102.71, 25.04]; % KunMing
city = [PEK; XAN; XJG; KMG];
for i = 1 : size(city, 1)
gmt(['psxy -R -J -Sc0.1c -K -O >> ' outfn '.ps'] , city(i, :));
txt = gmt('record', city(i, :) , ['city-' num2str(i)]);
gmt(['pstext -R -J -F+f9p,Times-Roman -K -O >> ' outfn '.ps'], txt);
end
% `=` means geographic vectors.
% With the modifier '+s', the input data should contain coordinates of start
% and end points
style = '=0.5c+s+e+a30+gblue+h0.5+p1p,blue';
data = [[PEK XAN]; [PEK XJG]; [PEK KMG]];
gmt(['psxy -R -J -S' style ' -W1.0p,blue -K -O >> ' outfn '.ps'], data);
% txt = {[110, 30], 'CIRCULAR'}; % not works
txt = gmt('record', [125, 40], 'GEOGRAPHIC');
gmt(['pstext -R -J -F+f13p,Helvetica-Bold,blue -O >> ' outfn '.ps'], txt);
%gmt(['pstext -R -J -O >> ' outfn '.ps'], txt);

gmt(['psconvert -A -Tg -E600 ' outfn '.ps']);
gmt(['psconvert -A -Tgf ' outfn '.ps']);
gmt('destroy');
open([outfn '.pdf']);


成图效果

GMT/MatlabToolbox的高级举例

我的更多文章

下载客户端阅读体验更佳

APP专享