因工作需求,需要绘制华北地区所有符合一定条件的地电阻率台站,过程如下。

  1. D:\MapSIS\A13系统_台站参数\ 目录下找到全国地电阻率台站参数文件 A国家_国家地震前兆台网中心_地电阻率.ORA

  2. 筛选华北地区地电阻率台站,并去除重复信息。

  3. 编写脚本,绘制华北地区地电阻率台站位置示意图。

  4. 进一步修改脚本,美化图件。

筛选台站

因为参数文件给出了所有台站的所有测项信息,但是我们只需要台站的经纬度和台站名,所以需要编写 bat 脚本提取华北地区地电阻率台站并且去除重复信息。

rem 以 | 为分隔符需要转义,去除重复行
gawk -F ^| "{if (NR>2) print $16,$17,$6}" A国家_国家地震前兆台网中心_地电阻率.ORA | gawk "!a[$0]++" > 全国.txt
rem 提取华北台站
gawk  "{if ($1>=109 && $1<=121 && $2>=33 && $2 <=43) print $1,$2,$3}" 全国.txt  > 华北.txt

绘制示意图

华北.txt 文件进行手工调整,删除不需要的台站,然后在 Windows 下编写 gmt 脚本绘制台站位置示意图(gmt 版本为 5.4.2)

gmt set MAP_TICK_LENGTH_PRIMARY 2p
gmt set FONT_LABEL 10p
set PS=华北-colorTopo.ps
set R=109/121/33/43
gmt grdcut earth_relief_01m.grd -R%R% -GcutTopo.grd
gmt grd2cpt cutTopo.grd -Cglobe -S-6800/8200/10 -D > colorTopo.cpt
gmt psbasemap -R%R% -JM5i -B2 -BWesN -K > %PS%
gmt grdimage cutTopo.grd  -R -J -CcolorTopo.cpt  -O -K >>%PS%
gmt psxy 华北-调格式.txt -R -J -St0.2c -Gblack -K -O >> %PS%
gmt pstext -R -J -F+f9p,39+jBC -Dj0/0.2c 华北-调格式.txt -K -O >> %PS%
gmt psscale -Bx3000+u"m" -By+l"Topo" -CcolorTopo.cpt -Dx1.3c/0+w4i/0.4c+o0/-1c+h -O --FONT_ANNOT_PRIMARY=10p >> %PS%
gmt psconvert %PS% -P -A -Tg -Z -C-sFONTPATH=C:/windows/fonts
del gmt*

效果见 华北地区台站分布1:

华北地区台站分布1

这种把华北平原弄成一片汪洋的配色方式被人吐槽了。。。将 850m 以下的配色手动修改了一下

华北.cpt
0 193/228/193 300 193/228/193
300 193/228/193 500 215/241/255
500 215/241/255 770 51/109/7
770 51/108/7 850 178/226/143

效果见 华北地区台站分布2:

华北地区台站分布2

起码不违反常识了。。华北地区的海拔高度大约为 -60~2848m ,以后有时间了再好好修改一下华北地区地形配色吧。