1、2023 年 2 月第 1 期城市勘测Urban Geotechnical Investigation SurveyingFeb2023No1引文格式:佘佐明，程晓光，王艳军 基于 GDAL 开源库的海量 DOM 分幅裁剪 J 城市勘测，2023(1):7478文章编号:16728262(2023)017405中图分类号:TP751文献标识码:A基于 GDAL 开源库的海量 DOM 分幅裁剪佘佐明1*，程晓光2，王艳军3*收稿日期:20221026作者简介:佘佐明(1965)，男，高级工程师，注册测绘师，主要从事大地测量、工程测量、摄影测量与遥感专业方面的研究与管理工作。Email:65448
2、0322 qqcom通信作者:程晓光(1985)，男，博士，高级工程师，主要从事航空摄影测量、遥感、测绘内业软件研发。Email:chengxg outlookcom基金项目:湖南省自然科学基金项目(2020JJ3020、2020JJ5164)(1.贵阳市测绘院，贵州 贵阳550000;2.飞燕航空遥感技术有限公司，湖北 武汉430061;3.湖南科技大学地理空间信息技术国家地方联合工程实验室，湖南 湘潭411201)摘要:海量数字正射影像(DOM)的分幅裁剪是测绘与 GIS 领域的基础问题。现有商业软件的 DOM 分幅裁剪功能存在自动化、批量化程度和效率低下的问题。针对海量 DOM 分幅裁剪
3、过程中涉及的源 DOM 预处理、图幅生成、数据部署、影像裁剪和结果检查等各个环节，提出了基于开源软件 GDAL 的解决方案。通过综合使用 GDAL 命令和开源软件编程，在最大程度上提高了海量 DOM 分幅裁剪的效率和自动化程度。本文方法在贵阳市航空航天摄影项目中得到了充分验证。关键词:海量 DOM;分幅裁剪;GDAL;开源软件1引言数字正射影像(Digital Orthophoto Map，DOM)是基础测绘领域主要的数据产品类型之一。在测区范围大、空间分辨率高时，DOM 数据量十分巨大，因而也通常称为海量 DOM。如果海量 DOM 以不分幅的形式存储，则在加载、显示、管理、调度、分析上十分困
4、难，还容易造成存储空间浪费。因此，一般需要对 DOM 进行分幅裁剪，方便管理与使用。GB/T 139892012 国家基本比例尺地形图分幅和编号1 规定了地形图分幅和编号的方法。除国家标准外，也可以自定义分幅方法。目前各类商业 GIS、遥感及摄影测量软件均支持对 DOM 基于空间坐标、像素坐标、矢量边界、空间关系等进行裁剪26。但是，DOM 裁剪往往存在源 DOM数量多，图幅数量大，影像格式、数据类型、波段数目、空间分辨率、坐标系、无数据值等属性不统一，坐标不对齐等问题。商业软件一般对单张 DOM、单个图幅的裁剪支持得很好，但是对海量 DOM、大量图幅的分幅裁剪的支持不够完善，操作烦琐，效率低
5、下。当前，测绘行业往往采用国外商业软件进行自动化、批量化的 DOM 分幅裁剪，例如，文献 7 9利用ArcGIS 的模型构建器制作批量裁剪工具;文献 10 利用 Java 和 ArcGIS 的二次开发组件库 ArcObject 开发裁剪工具;文献 11 利用 Arcpy 脚本修改 ArcGIS 裁剪函数模板，将模板应用到栅格数据裁剪。但国外商业软件价格昂贵，不利于测绘软件的自主可控。随着开源软件的兴起，借助开源软件进行海量DOM 的自动化分幅裁剪的难度大大降低。本文提出了基于开源软件 GDAL12(Geospatial Data AbstractionLibrary)对 DOM 进行分幅裁剪的
6、方法，通过在各步骤中综合使用 GDAL 命令和开源软件编程，实现了海量DOM 裁剪的自动化、批量化和高效化。贵阳市航空航天摄影项目验证了本文方法的有效性。2GDAL 简介GDAL 是 GIS 和遥感领域对栅格和矢量空间数据的读写支持得最好的开源软件之一12，提供了功能强大的处理命令可用于 DOM 裁剪。但大部分 GDAL 命令一般针对单个文件进行处理，在利用多个图幅裁剪时需要为各个图幅单独编写命令，在图幅多时十分不便。此外，GDAL 命令对于中文的支持不好。因此，海量 DOM的批量分幅裁剪需要引入 GDAL 的二次开发编程。除 C/C+的二次开发接口外，GDAL 还提供了 Py-thon，Ja
7、va，C#等的调用接口12。在 C+接口中，对栅格/矢量数据格式的支持是通过类 GDALDriver 实现的，栅格/矢量数据集抽象为类 GDALDataset，栅格数据集的单个波段抽象为类 GDALasterBand，空间坐标系抽象为类 OGSpatialeference，GDALOpen 和 GDAL-Close 函数用于打开和关闭数据集，CPLMalloc 和 CPL-第 1 期佘佐明，程晓光，王艳军.基于 GDAL 开源库的海量 DOM 分幅裁剪Free 常用于申请和释放内存。3方法3.1源 DOM 预处理由于 GDAL 几乎支持所有的主流栅格数据格式，且 GDALDataset 屏蔽了
8、格式差异，使用相同的方法即可读取绝大部分格式的源 DOM12，因此一般不需要考虑文件格式统一的问题。源 DOM 只需要做到数据类型、波段、坐标系、空间分辨率和无数据值的统一及像素坐标对齐。(1)数据类型统一DOM 常用 1 字节无符号整型(下文简称为 UInt8_t)存储单个波段的 DN 值。但是，DOM 的数据类型不限于此，例如，卫星遥感影像常采用 2 字节无符号整型。在 gdal_translate 命令中，可通过 ot 选项设置输出的数据类型，通过 scale 和 exponent 选项采用幂函数将位于指定范围的 DN 值拉伸到指定的输出范围12。gdal_trans-late 对应的
9、C/C+接口函数为 GDALTranslate()。(2)波段统一常见的真彩色 DOM 一般包含 3 个波段。对只有1 个波段、采用颜色查找表记录 DN 值的彩色 DOM，可在 gdal _ translate 命令中设置 expand 选 项 的 值 为rgb12。该命令的 b 选项可用于设置输出波段、调整波段顺序。(3)坐标系统一源 DOM 可能具有不同的空间坐标系。编程实现影像坐标转换较为复杂，可使用 gdalwarp 命令将源DOM 的坐标系统一到指定坐标系。gdalwarp 命令的 r选项用于指定重采样方法，s_srs 和 t_srs 选项分别用于设置源 DOM 和目标 DOM 坐标
10、系12。指定坐标系的方式包括 EPSG 编码、WKT(Well Known Text)文本、POJ4 声明或含 WKT 格式的坐标系定义的 prj 文件。作者建议，对 EPSG 已定义的坐标系，优先使用 EPSG编码;否则可在与目标坐标系接近且 EPSG 已定义的坐标系的 WKT 的基础上进行修改。(4)空间分辨率统一源 DOM 的空间分辨率统一非常重要。gdal_trans-late 和 gdalwarp 都有 tr 选项可用于设置输出分辨率12。作者建议将 tr 选项设置为确定的值，非必要不设置 ts、outsize、a_ullr 选项。(5)无数据值统一为表示某像素是否有效，测绘领域常用
11、的方法是将无效像素的 DN 值设为无数据(NoData)值。gdalwarp 命令的 dstnodata 选项和 gdal_translate 的 a_nodata 选项可用于设置输出影像的无数据值，既可以所有波段共用相同的无数据值，也可以逐波段设置12。但该方法只是设置影像的无数据值属性，并不能修改像素 DN 值。针对 DOM 单个波段具有多个无数据值以及影像存在伪无数据区域的问题，以数据类型为 UInt8_t 的DOM 为例，本文给出一种解决方法，可统一无数据值，并排除伪无数据区域，步骤如下:将所有可能的无数据值放入无重复集合 Snodata;创建一个单通道、行列数和源 DOM 相同的二值
12、图像 Mbin，值初始化为 0;遍历 DOM 所有像素，查找各波段 DN 值均在Snodata内的像素，Mbin中对应像素的值设为 255;对 Mbin进行联通分析，获得联通区域;遍历各个联通区域，采用图 1 所示方法进行处理。区域邻接像素计算方法如下:采用 33 的结构元对 Mbin进行膨胀，将膨胀结果减去 Mbin，差值图像中值为 255 的像素即为区域邻接像素。Abig可简单取整幅DOM 面积 5%左右的值，DNhigh可取大于等于 200 的值，DNlow可取小于等于 20 的值，phigh和 plow可取大于等于 50%的值。图 1无数据值修复流程图(6)坐标对齐GDAL 定义的像素
13、坐标对齐是确保图像角点 XY坐标是分辨率的整数倍，在 gdalwarp、gdalbuildvrt 和gdal_merge py 中，tap 选项可用于设置像素坐标对齐12。但如果 DOM 无旋转，且分辨率已统一，则可以认为像素坐标对齐是所有 DOM 的角点坐标差为分辨率的整数倍。此处给出一种对源 DOM 进行坐标对齐57城市勘测2023 年 2 月的方法。设基准坐标为 xbase，ybase()，源 DOM 左上角坐标为 xul，yul()，则对齐后的左上角坐标 xul，yul()的计算方法可见式，其中 round 为四舍五入函数，x和 y分别为源 DOM 在列和行方向上的分辨率。xul=xr
14、ound(xulxbase)/x+xbaseyul=yround(yulybase)/y+ybase(1)3.2图幅生成文献 13给出了一种基于测区边界进行地图分幅的方法，可以自动、快速地生成位于测区内部和边界上的图幅。3.3数据部署鉴于数据量大，DOM 可能被存储在多台电脑的外存储器上。虽然可通过网络读写 DOM，但速度受带宽制约。显然，读写本地数据、在多台电脑上并行裁剪的方式效率高，在单台电脑上应尽量使用只覆盖本地源DOM 的图幅进行裁剪。此处给出一种方法:(1)使用 gdaltindex12 或其他工具生成含有本地源 DOM 覆盖范围的矢量文件;(2)将前一步生成文件与图幅接合表叠加显示
15、;(3)将完全落在(需要考虑到图幅外扩)本地源DOM 覆盖范围内的图幅单独保存为一个图幅接合表，用于本地 DOM 裁剪;(4)将源 DOM 分布在多台电脑上的图幅保存为图幅接合表，通过网络裁剪。3.4分幅裁剪DOM 分幅裁剪的流程如图 2 所示。如果求得所有源 DOM 包容盒的并集并生成对应的矩形图幅，则该流程得到的就是拼接后的图幅。图 2影像裁剪流程图(1)重叠集定义DOM 分幅裁剪的核心问题是计算图幅对应的目标 DOM 的像素来自于哪些源 DOM 的哪些行和列，本文称为重叠集。重叠集抽象出的结构体 OverlapSet 的示例定义代码如下:struct OverlapSetstring s
16、trSrcPath;/源 DOM 路径string strDstPath;/目标 DOM 路径int nSrcow0;/在源 DOM 的起始行int nSrcCol0;/在源 DOM 的起始列int nDstow0;/在目标 DOM 的起始行int nDstCol0;/在目标 DOM 的起始列int nowNum;/重叠的行数int nColNum;/重叠的列数;(2)重叠集计算设某图幅的二维包容盒为 Eb，某源 DOM 的二维包容盒为 Es，计算二者的交集 Eo。如果 Eo为空，则该源DOM 不含构成该图幅的像素;否则，需要计算重叠集。单个像素跨越图幅边界是裁剪中的常见问题。如果图幅不外扩，一种处理方法是将它们作为所有跨越的图幅的一部分，即生成 DOM 的覆盖范围是包含图幅的最小矩形，另一种处理方法是采用四舍五入的原则确定它们属于哪一个跨越的图幅。前者会导致相邻图幅间存在一行或一列像素的重叠，称为有重叠法，后者不会，称为无重叠法。无重 叠 法 计 算 目 标 DOM 的 左 上 角 坐 标xd，min，yd，max()的方法可见式(2)，有重叠法的计算方法可见式(3)，其中 xt，m
