如何沿着某一直线绘制特定物理场（高程、负荷场）剖面（Arcgis\GMT\MATLAB）

目录

一、如何在DEM中绘制特定测线的剖面

二、如何在其他物理场中绘制特定测线的剖面（以负荷场为例）

在GIS中，我们通常需要绘制特定剖面线上的值的变化情况，剖面可以清晰反映出某个物理场具体沿着特定方向的变化情况，比如以下的情况就非常的易于理解（图片来自于Pan et al, JGR）。

---

接下来我们就针对于这样的需求来分析我们可以采用什么方法实现？我将尽可能采用多的方法和软件工具达到这一目的。首先我们先从直观的高程入手。

一、如何在DEM中绘制特定测线的剖面

【数据准备】我是直接利用GMT5工具将全球DEM提取出中国地区的DEM。具体的下载地址如下：

https://www.ncei.noaa.gov/products/etopo-global-relief-model

在gmt中，采用以下的命令得到了中国大陆地区的DEM：resamp_china.grd

gmt grdcut ETOPO1_Bed_g_gdal.grd -R%R% -Gchina.grd
gmt grdsample china.grd -Gresamp_china.grd -I0.1/0.1

同理，我还利用本人自带的负荷场数据，得到了中国大陆地区的负荷场数据：1.grd,用于第二节的使用。

xyz2grd  green_orgin.txt -R%R% -I0.125  -G1.grd
gmt grdsample 1.grd -G1.grd -I0.1/0.1

【数据转换】由于gmt得到的grd格式的文件无法直接打开，且都是乱码，因此需要在Global Mapper中进行数据格式的转换，具体的操作如下：打开软件，选择对应的grd文件，载入数据，然后右键点击，选择Layer，Export数据至xyz格式，得到我们需要的数据源。

【gmt直接绘制测线剖面】本文的测线剖面均采用下图所示

绘图的代码：采用GMT6平台，具体的参考代码见gmt中文社区。

REM!/usr/bin/env bash
REM
REM 绘制地形起伏剖面图
gmt begin profile png
    rem 边框参数
    gmt set MAP_FRAME_PEN 0.2p,black
    rem 刻度相关参数
    gmt set MAP_TICK_PEN 0.5p,black
    gmt set MAP_TICK_LENGTH 2P
    gmt set LABEL_FONT 2p
    gmt set FONT_LABEL 8p
    gmt set LABEL_FONT_SIZE 3p
    gmt set ANNOT_FONT_SIZE_PRIMARY 3p
    gmt set MAP_FRAME_WIDTH 5p
    gmt set MAP_FRAME_TYPE=plain
    gmt set FORMAT_ANNOT_PRIMARY 0.5p,Times-Roman,black
    REM 绘制地形图
    gmt makecpt -Cterra -T-100/8000/100 -Z -D

    gmt basemap -JM12c -R75/106/18/42 -Ba5f5a -BWSen
    gmt grdimage @earth_relief_01m -I+d
    gmt coast -N2 -A1000 -Df -Lg79/22+c40+w500k+f+u

    REM 绘制自定义指北针
    echo 103 40 0 | gmt plot -Skcompass/0.5 -W0.6p
    echo 103 39 N | gmt text -F+f10p,4
    REM 选取测线AB
    echo 90 30 A > tmp1
    echo 99 40 B >> tmp1
    REM 绘制测线AB
    gmt plot tmp1 -W1p,red
    REM 标注AB
    gmt text tmp1 -F+f15p -D0c/0.2c

    REM 沿测线AB绘制地形高度
    gmt basemap -R0/12/0/7000 -JX12c/3c -Bya2000+l"Elevation (m)" -Bxa5f5+l"Distance"+u"\260" -BWSrt -Y12c
    REM 标注AB位置
    echo 0 7000 A | gmt text -F+f10p+jBC -N
    echo 10 7000 B | gmt text -F+f10p+jBC -N
    REM 沿测线提取地形高度
    gmt project -C90/30 -E99/40 -G0.1 | gmt grdtrack -G@earth_relief_02m > tmp2
    REM 将海平面以下填充为淡蓝色
    echo 0 0 > tmp3
    echo 10 0 >> tmp3
    gmt plot tmp3 -Wblack -Gblue -L+y-8000
    REM 将地形填充为灰色
    gmt plot tmp2 -i2,3 -Wblack -Ggray -L+y-8000
gmt end

需要注意的地方有：

1.如何选取测线，A（90，30），B（99，40）

REM 选取测线AB
    echo 90 30 A > tmp1
    echo 99 40 B >> tmp1

2.沿测线提取测线高度

REM 沿测线提取地形高度
gmt project -C90/30 -E99/40 -G0.1 | gmt grdtrack -G@earth_relief_02m > tmp2

我们观察一下tmp2文件，其中的第一列是经度，第二列是纬度，第三列是剖面的编号（对应剖面的x轴）第四列是高程（对应剖面的y轴）。

我们在matlab中可以读取tmp2文件，绘制出剖面线。

A = load('tmp2')
x = A(:,3);
zdata = A(:,4);
plot(x,zdata,'r')
set(gca,'gridlinestyle',':','LineWidth',0.5,'GridAlpha',0.5)
xlabel('Distance', 'Fontname', 'Times New Roman', 'Fontsize',13);
ylabel('Elevation/m', 'Fontname', 'Times New Roman', 'Fontsize',13);
grid on

接下来我们看看Arcgis中是如何实现的。在arcgis中读入转换导出的dem.tif数据【在Global Mapper中将grd导出为tif格式，便于Arcgis处理】。

【在Arcgis中创建剖面】由于没有找到加载外部线的功能，因此需要我们自己添加一个剖面线的shp文件。注意在创建折线的时候，我们需要使用【绝对X\Y，此处输入在gmt中设置的起始点坐标A（90，30）B（99，40）】

完成剖面，下一步选择Arcgis的3D Analyst工具，

关于arcgis剖面线，可参考：https://blog.csdn.net/mrib/article/details/120289686

最后在MATLAB中实现了DEM绘制，绘制剖面，并插入子图显示剖面高程变化

二、如何在其他物理场中绘制特定测线的剖面（以负荷场为例）

利用上述导出的中国大陆地区的地表负荷场数据，我们可以利用matlab绘制如下的结果

其余的方法如DEM的一致。