以前GMT画图,用的都是90米或30米的数字高程模型(DEM),最近在搞InSar监测形变,偶然发现居然有12.5米的DEM,而且可以免费获得。网上一艘,好多介绍。但是,我一般不用windows系统,几乎都在Linux下干活,而网上说的下载下来后需要用ArcGis软件处理拼接等,我这强迫症就犯了,必须在Linux下搞定所有处理过程。
一,12.5米地形数据下载
这个按网上的方法下载,网址:ASF Data Search
怎么下载就不多说了,网上搜12.5m DEM,一大堆。这里简单说下,就是选ALOS PALSAR影像,然后在filter里筛选下:File type选 Hi-Res Terrain Corrected,这个应该就是ALOS的高分辨率地形修正文件,然后Beam Model选FBS和FBD。
这里注意下,如果你想要的区域很大,会有很多文件,确保相互重叠覆盖(看图选就行),每个小区域选一幅数据就可以了,数据非常大,下载速度看自己网络了。
二,格式转换
下载的文件解压后,*.dem.tif这个文件就是DEM。为了便于后续处理,要将tif转换为grd格式
用gdal,没有的装下 sudo yum install gdal,
格式转换: gdal_translate -of GMT file.dem.tif file.dem.grd
三,WGS84地心直角坐标转经纬度
原始地形数据是WGS84直角坐标,即(X米,Y米),为了便于使用,我一般转为经纬度。
这里我全用GMT处理。
核心命令是grd2xyz和mapproject。
gmt grd2xyz $file -di0 | gmt mapproject -Ju$zone/1:1 -I -C -F --PROJ_SCALE_FACTOR=0.9996 > $name.xyz
先用grd2xyz转成ascii文件,这里注意,要用-di0选项,否则没有数据的地方会输出0,这样DEM最小值会出现0,我们需要输出NaN。
mapproject命令用-Ju进行坐标系转换,这里注意,-Ju$zone,zone为你选区域的投影带,带号在解压文件中 .aux.xml文件里,可以自己查看,想脚本处理的可以:
zone=`grep " <SRS>" ./$file.aux.xml | awk '{printf("%s\n",$6)}'`
zone=${zone:0:2}
读出来。然后在用xyz2grd转会二进制grd文件,这一步的shell脚本如下,自己根据情况稍微改下就行(这里用多线程处理)。
if [ $# -ne 1 ];
then
echo ""
echo "Usage: wgs2ll.sh thread_num"
echo ""
exit
fi
tmp_fifofile="/tmp/$$.fifo"
mkfifo $tmp_fifofile # 新建一个FIFO类型的文件
exec 6<>$tmp_fifofile # 将FD6指向FIFO类型
rm $tmp_fifofile #删也可以,
thread_num=$1 # 定义最大线程数
for ((i=0;i<${thread_num};i++));do
echo
done >&6
t1=`date`
for file in *.dem.grd
do
read -u6
{
name=${file%%.*}
echo $name
zone=`grep " <SRS>" ./$file.aux.xml | awk '{printf("%s\n",$6)}'`
zone=${zone:0:2}
xnum=`gmt grdinfo $file -C | awk '{printf("%d",$10);}'`
ynum=`gmt grdinfo $file -C | awk '{printf("%d",$11);}'`
echo "WGS 84 / UTM zone: $zone;" xnum=$xnum ynum=$ynum
gmt grd2xyz $file -di0 | gmt mapproject -Ju$zone/1:1 -I -C -F --PROJ_SCALE_FACTOR=0.9996 > $name.xyz
string=`gmt minmax $name.xyz -C`
xmax=`echo $string | awk '{printf("%s",$2);}'`
xmin=`echo $string | awk '{printf("%s",$1);}'`
ymax=`echo $string | awk '{printf("%s",$4);}'`
ymin=`echo $string | awk '{printf("%s",$3);}'`
dx=`echo $xmax $xmin $xnum | awk '{printf("%.9lf",($1-$2)/$3);}'`
dy=`echo $ymax $ymin $ynum | awk '{printf("%.9lf",($1-$2)/$3);}'`
echo lonmax=$xmax lonmin=$xmin latmax=$ymax latmin=$ymin dx=$dx dy=$dy
gmt xyz2grd $name.xyz -G$name.grd -R/$xmin/$xmax/$ymin/$ymax -I$dx/$dy
rm $name.xyz
echo >&6
} &
done
wait
exec 6>&- # 关闭FD6
t2=`date`
echo $t1 $t2
三,裁剪和拼接
这一步最麻烦,要根据自己情况处理。总体思路就是,先横向(同一纬度方向)处理,将同一纬度的几个数据裁剪拼接。裁剪的时候,下端纬度边界为这一排数据中纬度下界最大值。
这一步就形成几个“横条“,然后把所有“横条“的东、西两边经度裁剪齐,最后把横条拼接起来就好了。
这一步比较麻烦,因为根据不同情况,需要自己写脚本排序等确定经纬度边界。
主要用grdsample,统一采样率,不然grdpaste不能用。
grdcut,用来将数据裁剪为有共同的经纬度边界。
grdpaste用来经有共同经度或纬度边界的数据拼接起来。
生成多个横条*.grd,最后把横条的经度裁剪齐,然后grdpaste在一起就行了。