12.5米DEM下载、格式转换、裁剪与拼接

       以前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在一起就行了。