gmt-template/gmtxy-image2.sh
2023-03-19 19:41:55 +08:00

308 lines
11 KiB
Bash
Executable File
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#!/bin/bash
# 包含GMT自带脚本文件 其中包含了一些有用的功能 比如获取网格文件的范围
. gmt_shell_functions.sh
# 包含dispOption脚本
. dispOptions.sh
. generic_func.sh
# 定义一个函数 执行输入的语句或者将其显示在屏幕上
RunOrEcho()
{
first_str=`echo ${2// /''}`
if [[ x${3} != x ]]; then
sec_str=`echo ${3// /''}`
fi
if [[ ${1} == 1 ]]; then
${first_str}
if [[ x${3} != x ]]; then
${sec_str}
fi
else
if [[ x${3} != x ]]; then
printf "%s\n%s\n" "${first_str}" "${sec_str}"
else
printf "%s\n" "${first_str}"
fi
fi
}
# GMT显示平面数据脚本输入文件为网格文件没有包含网格化语句因为网格化过程中的情况多样化建议在其他脚本中个别添加再调用此脚本
# 准备几个常用的排版参数组 记录的参数包括 图片的宽度 默认坐标标示设定 色标位置 色标的大小 标示字体大小
com_layouts=("1.5i,WesNZ,0.1i/-0.2i,1.3i/0.05i+h,10.5p" \
"1.5i,WeSnZ,1.65i/0.05i,1.3i/0.05i,10.5p" \
"1.2i,WesNZ,0.1i/-0.1i,1.0i/0.03i+h,7.5p" \
"1.2i,WeSnZ,1.3i/0.05i,1.0i/0.05i,7.5p")
# 初始化参数
data='null'
outpsfile='null'
gridData='null'
unit='m'
color='rainbow'
range='null'
overwriteRange='null'
overwrite=0
labels=("x (m)" "y (m)")
plotgrad=0
axistick=("a" "a") #横纵坐标轴间隔 a表示自动
bartick="a" #色标轴标轴间隔 a表示自动
frameset="null"
plot_colorbar=1
annotation="null"
keep_open=0
from_open=">"
image_origin=("0.5" "0.5")
bar_origin=("0" "0")
polyfile='null'
lay_num=0
run_command=1
close_poly='-L'
# 从命令行获取参数
while getopts "hi:o:L:r:u:c:l:G:t:v:f:a:p:X:x:gnKOSP" arg
do
case $arg in
h)
dispTitle "${0##*/}" "Image mapping of grid file(s) using the GMT software under the x-y plane. \
This script accepts one or several grid (.nc or .grd) files as inputs and outputs \
image files (.png and .eps files). The script could also output data specified shell scripts for \
further modification. For detailed explanations of relevant arguments, please look for the GMT's manuscripts. \
This software comes with ABSOLUTE NO WARRANTY. Use it at your own caution."
dispAuthorInfo "Dr. Yi Zhang (yizhang-geo@zju.edu.cn). Contact me if you have questions or suggestions."
dispUsage "${0##*/} -i<grid-data> [-o<img-file>] [-L<number>] [-r<xmin>/<xmax>/<ymin>/<ymax>] [-u<unit>] \
[-c<cpt-file>] [-t<number>,<number>] [-v<number>] [-l<x-label>;<y-label>] [-g] [-G<grad-file>] \
[-f<WesNZ>] [-p<poly-file>] [-n] [-a<x-coor>,<y-coor>,<annotation>] [-K] [-O] [-S] \
[-X<x-offset>,<y-offset>] [-x<x-offset>,<y-offset>]"
dispOptionShort "-i" "Input grid file. This is the only mandatory option that the script has. Just enter a filename \
and you get an image."
dispOptionShort "-o" "Output image files. Name extension is not needed. the input filename will be \
used if this option is not set. NOTE: set this option for plotting multiple images."
dispOptionShort "-L" "Select layout template that will be used for plotting. The script has four inbuilt templates \
you can choose from (0 (default) - 3). Please try out to see their difference. You can also add more templates by \
editing the script. Just make sure that you know what you are doing."
dispOptionShort "-r" "Range of the input data. The script will detect the full range of the data automatically and \
the use of this option will overwrite the argument."
dispOptionShort "-u" "Data unit. The default is meter (m)."
dispOptionShort "-c" "Filename of the color scale. The default is rainbow.cpt."
dispOptionShort "-t" "Intervals of axis's ticks. the script will set the intervals automatically \
if this option is not set."
dispOptionShort "-v" "Intervals of color bar's labels. the script will set the intervals \
automatically if this option is not set."
dispOptionShort "-l" "Axis' labels separated by semicolons. The defaults are x (m) and y (m)."
dispOptionShort "-g" "Plot a over layer of directional illumination. The default is false."
dispOptionShort "-G" "Use a input grid data for applying directional illumination. \
This should be used with the '-g' option at the same time."
dispOptionShort "-f" "Set frames for plotting axises. The default is WesNZ for layout template 0."
dispOptionShort "-n" "Don't plot color bar."
dispOptionShort "-a" "Add an annotation on the map. you need to enter a string contains both \
coordinates and text of the annotation, such as '10,10,(a)'. For more than one annotations, \
please output a data specified script using the -S option, then edit the script accordingly."
dispOptionShort "-p" "Take inputs from a table file and plot polygons. Each line of the file \
represents a 2D point's coordinates."
dispOptionShort "-P" "Keep the polygon plotted by the -p option open."
dispOptionShort "-K" "Keep the .ps file open. This option must be set for plotting multiple images to a single output file."
dispOptionShort "-O" "Continue from a previous .ps file. This option must be set for plotting multiple images to a single output file."
dispOptionShort "-X" "Move the starting point for mapping the image. Use this option for plotting sub-figures."
dispOptionShort "-x" "Move the starting point for mapping the color bar. Use this option for plotting sub-figures."
dispOptionShort "-S" "Instead of running the script for actually plotting, print an executable shell script on screen. \
This option is designed to help building data specified shell scripts for further modification."
exit 0;;
i)
data=$OPTARG;;
o)
outpsfile=$OPTARG;;
L)
lay_num=$OPTARG;;
u)
unit=$OPTARG;;
c)
color=$OPTARG;;
t)
axistick=(${OPTARG//,/ });;
X)
image_origin=(${OPTARG//,/ });;
x)
bar_origin=(${OPTARG//,/ });;
v)
bartick=$OPTARG;;
f)
frameset=$OPTARG;;
r)
overwrite=1
overwriteRange=$OPTARG;;
l)
OLD_IFS="${IFS}"
IFS=";"
labels=(${OPTARG})
IFS="${OLD_IFS}";;
g)
plotgrad=1;;
n)
plot_colorbar=0;;
G)
gridData=$OPTARG;;
a)
annotation=$OPTARG;;
p)
polyfile=$OPTARG;;
P)
close_poly='';;
K)
keep_open=1;;
O)
from_open="-O >>";;
S)
run_command=0;;
?)
printf "error: unknow argument\nuse -h option to see help information\n"
exit 1;;
esac
done
# 进行必要的参数检查
if [[ $data == "null" ]]; then
printf "error: no input file name\nuse -h option to see help information\n"
exit 1
else
# 初始化临时文件名
cptfile=user.cpt
if [[ ${outpsfile} == 'null' ]]; then
psfile=${data%.*}.ps
jpgfile=${data%.*}.png
else
psfile=${outpsfile}.ps
jpgfile=${outpsfile}.png
fi
# 获取排版
layouts=(${com_layouts[${lay_num}]//,/ })
# 获取坐标轴设置
if [[ ${frameset} == 'null' ]]; then
frameset=${layouts[1]}
fi
# 获取网格范围
if [[ $overwrite == 1 ]]; then
range=${overwriteRange}
else
range=$(gmt_get_gridregion ${data})
fi
#根据横纵坐标的范围计算图片的长和高 我们默认图片宽度为1.5i
range_coor=(${range//// })
pic_height=`echo "scale=4; ${layouts[0]%i} * (${range_coor[3]} - ${range_coor[2]})/(${range_coor[1]} - ${range_coor[0]})"|bc`
# 输出一个bash脚本的头部到屏幕
if [[ ${run_command} == 0 && ${from_open} == '>' ]]; then
echo "#!/bin/bash"
fi
# 设置绘图参数
if [[ ${from_open} == '>' ]]; then
RunOrEcho ${run_command} "gmt gmtset \
FONT_ANNOT_PRIMARY=${layouts[4]},Times-Roman,black \
MAP_FRAME_PEN=thinnest,black \
MAP_GRID_PEN_PRIMARY=thinnest,black \
MAP_TICK_PEN_PRIMARY=thinnest,black \
MAP_TICK_LENGTH_PRIMARY=1p/0.5p \
MAP_TITLE_OFFSET=7.5p \
MAP_GRID_CROSS_SIZE_PRIMARY=2p \
FONT_LABEL=${layouts[4]},Times-Roman,black \
MAP_LABEL_OFFSET=2.5p \
MAP_ANNOT_OFFSET_PRIMARY=2.5p"
fi
RunOrEcho ${run_command} "gmt gmtset MAP_FRAME_AXES=${frameset}"
# 处理-B选项的参数
if [[ ${axistick[0]} == 'a' ]]; then
axistick[0]=${axistick[0]}g
axistick[1]=${axistick[1]}g
else
axistick[0]=${axistick[0]}g${axistick[0]}
axistick[1]=${axistick[1]}g${axistick[1]}
fi
RunOrEcho ${run_command} "gmt grd2cpt ${data} -C${color} -R${range} -Z -D > ${cptfile}"
if [[ $plotgrad == 1 ]]; then
gradfile=${data%.*}Grad.nc
if [[ $gridData == 'null' ]]; then
RunOrEcho ${run_command} "gmt grdgradient ${data} -G${gradfile} -Nt -A0/45"
else
RunOrEcho ${run_command} "gmt grdgradient ${gridData} -G${gradfile} -Nt -A0/45"
fi
if [[ ${run_command} == 0 ]]; then
echo "gmt grdimage ${data} -R${range} -C${cptfile} -I${gradfile} \
-Bx${axistick[0]}+l\"${labels[0]}\" -By${axistick[1]}+l\"${labels[1]}\" \
-JX${layouts[0]}/${pic_height}i -X${image_origin[0]}i -Y${image_origin[1]}i -K -P ${from_open} $psfile"
else
gmt grdimage ${data} -R${range} -C${cptfile} -I${gradfile} \
-Bx${axistick[0]}+l"${labels[0]}" -By${axistick[1]}+l"${labels[1]}" \
-JX${layouts[0]}/${pic_height}i -X${image_origin[0]}i -Y${image_origin[1]}i -K -P ${from_open} $psfile
fi
else
if [[ ${run_command} == 0 ]]; then
echo "gmt grdimage ${data} -R${range} -C${cptfile} -Bx${axistick[0]}+l\"${labels[0]}\" \
-By${axistick[1]}+l\"${labels[1]}\" -JX${layouts[0]}/${pic_height}i \
-X${image_origin[0]}i -Y${image_origin[1]}i -K -P ${from_open} $psfile"
else
gmt grdimage ${data} -R${range} -C${cptfile} -Bx${axistick[0]}+l"${labels[0]}" \
-By${axistick[1]}+l"${labels[1]}" -JX${layouts[0]}/${pic_height}i \
-X${image_origin[0]}i -Y${image_origin[1]}i -K -P ${from_open} $psfile
fi
fi
# 画多边形
if [[ $polyfile != 'null' ]]; then
RunOrEcho ${run_command} "gmt psxy ${polyfile} -JX${layouts[0]}/${pic_height}i -W0.25p,black,- -R${range} ${close_poly} -K -O >> $psfile"
fi
#添加号码
if [[ ${annotation} != 'null' ]]; then
if [[ ${run_command} == 0 ]]; then
echo "gmt pstext -R${range} -JX${layouts[0]}/${pic_height}i -K -O <<- EOF >> $psfile"
echo ${annotation}
echo "EOF"
else
gmt pstext -R${range} -JX${layouts[0]}/${pic_height}i -K -O <<- EOF >> $psfile
${annotation}
EOF
fi
fi
#-C${cptfile}+Uk 使用km色标单位除1000
#如果unit等于km则在cptfile后面添加+Uk
if [[ ${plot_colorbar} == 1 ]]; then
if [[ ${unit} == 'km+Uk' ]]; then
RunOrEcho ${run_command} "gmt psscale -Dx${layouts[2]}+w${layouts[3]} -C${cptfile}+Uk -Bx${bartick} -By+l${unit} \
-X${bar_origin[0]}i -Y${bar_origin[1]}i -O -K >> $psfile"
else
RunOrEcho ${run_command} "gmt psscale -Dx${layouts[2]}+w${layouts[3]} -C${cptfile} -Bx${bartick} -By+l${unit} \
-X${bar_origin[0]}i -Y${bar_origin[1]}i -O -K >> $psfile"
fi
fi
if [[ ${keep_open} == 0 ]]; then
# 这一句什么都不做 唯一的作用就是关闭ps文件
if [[ ${run_command} == 0 ]]; then
echo "gmt pstext -R${range} -JX -O <<- EOF >> $psfile"
echo "EOF"
else
gmt pstext -R${range} -JX -O <<- EOF >> $psfile
EOF
fi
# 输出 eps 和 png 文件
RunOrEcho ${run_command} "gmt psconvert $psfile -A -TEG -E300"
# 删除临时文件 使用linux终端rm命令
RunOrEcho ${run_command} "rm $cptfile $psfile gmt.history gmt.conf"
if [[ $plotgrad == 1 ]]; then
RunOrEcho ${run_command} "rm $gradfile"
fi
# 在终端显示图像 此命令需要imgcat.sh脚本和iTerm终端
#imgcat $jpgfile
# 打开图片文件 此命令使用MacOS终端open命令 和Linux下的xdg-open命令
if [[ ${run_command} == 1 ]]; then
view_file ${jpgfile}
fi
fi
fi