一种利用gdal命令行提取栅格像元值生成等值面的方法

概述

地学工作者在工作中常常会用到根据栅格文件提取等值线作图的工作，如利用地形DEM文件提取栅格等工作。目前主流的GIS软件都有相关的分析提取等值线功能。但是有时出于渲染地图效果的目的，需要将等值线转换为矢量面文件，这时有一部分桌面GIS软件自带的功能实现起来比较麻烦，尤其是在图幅、要素很多，需要多次提取的时候。

本文介绍了gdal命令行程序来完成提取等值面并添加属性的工作方法。

实现方法

gdal是开源地理空间基金会（OsGeo）发布的一个处理栅格和矢量空间数据的编译库。具有供多种编程语言调用的API，同时也可以使用命令行程序完成一些不复杂的操作。

本方法需要使用gdal 2.4版本及以上，可以在官网下载源代码编译下载地址Download — GDAL documentation。也可以下载编译好的程序。作者是C++小白，所以选择了编译好的程序，可以在这个地址下载​GISInternals。本文使用到的程序在\bin\gdal\apps目录下。可以选择设置环境变量，或者直接cd到该目录执行。

根据gdal官网介绍，gdal_contour是用来提取等值线的工具，有很多命令行参数。其中参数 -a 为添加间距值的属性； -p为生成等值面，但是-p生成的等值面时会忽略-a参数，这种情况下需要使用-amin来指定面要素上的最小值对应的属性名称，用-amax指定面要素上最大值对应的属性名称。

gdal_contour -b 1 -amin minELEV -amax maxELEV -p -i 20.0 -f "ESRI Shapefile"  dem30.tif contour.shp

以上命令的意思为：

读取当前目录下的dem30.tif的第1个波段值为间距20的等值面，文件指定格式为ESRI Shapefile，每个等值面要素的下限字段名称为minELEV、上限字段名称为maxELEV，在当前目录下生成contour.shp

命令执行后，屏幕会输出如下字样

0...10...20...30...40...50...60...70...80...90...100 - done.

至此，根据gdal命令行提取一个等值面的操作就完成了。可以用QGIS查看效果。个人认为在大量处理等值线图时，命令行的效果要好于GIS软件自带工具箱操作。

图形效果

 属性效果

批量处理

若要多次、批量操作，可以用操作系统的命令行脚本调用gdal_contour实现。Windows系统下可以使用批处理脚本实现。

for /f  %%i in ('dir *.tif /b') do (gdal_contour -b 1 -amin minELEV -amax maxELEV -p -i 20.0 -f "ESRI Shapefile"  %%i %%i.contour.shp )

以上是一个十分简易的批量提取多个栅格文件的处理脚本，不是很灵活，可以复制到记事本保存为bat后缀，然后复制到栅格文件所在的目录下执行。

题外话

制图过程中有可能需要等值线与等值面结合使用，所以这里也写一下生成等值线的命令。

gdal_contour -b 1 -a ELEV -i 20.0 -f "ESRI Shapefile"  dem30.tif Test.shp