基于GDAL开源库的海量DOM分幅裁剪_佘佐明.pdf

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:654480322 qqcom通信作者:程晓光(1985)，男，博士，高级工程师，主要从事航空摄影测量、遥感、测绘内业软件研发。Email:chengxg outlookcom基金项目:湖南省自然科学基金项目(2020JJ3020、2020JJ5164)(1.贵阳市测绘院，贵州 贵阳550000;2.飞燕航空遥感技术有限公司，湖北 武汉430061;3.湖南科技大学地理空间信息技术国家地方联合工程实验室，湖南 湘潭411201)摘要:海量数字正射影像(DOM)的分幅裁剪是测绘与 GIS 领域的基础问题。现有商业软件的 DOM 分幅裁剪功能存在自动化、批量化程度和效率低下的问题。针对海量 DOM 分幅裁剪过程中涉及的源 DOM 预处理、图幅生成、数据部署、影像裁剪和结果检查等各个环节，提出了基于开源软件 GDAL 的解决方案。通过综合使用 GDAL 命令和开源软件编程，在最大程度上提高了海量 DOM 分幅裁剪的效率和自动化程度。本文方法在贵阳市航空航天摄影项目中得到了充分验证。关键词:海量 DOM;分幅裁剪;GDAL;开源软件1引言数字正射影像(Digital Orthophoto Map，DOM)是基础测绘领域主要的数据产品类型之一。在测区范围大、空间分辨率高时，DOM 数据量十分巨大，因而也通常称为海量 DOM。如果海量 DOM 以不分幅的形式存储，则在加载、显示、管理、调度、分析上十分困难，还容易造成存储空间浪费。因此，一般需要对 DOM 进行分幅裁剪，方便管理与使用。GB/T 139892012 国家基本比例尺地形图分幅和编号1 规定了地形图分幅和编号的方法。除国家标准外，也可以自定义分幅方法。目前各类商业 GIS、遥感及摄影测量软件均支持对 DOM 基于空间坐标、像素坐标、矢量边界、空间关系等进行裁剪26。但是，DOM 裁剪往往存在源 DOM数量多，图幅数量大，影像格式、数据类型、波段数目、空间分辨率、坐标系、无数据值等属性不统一，坐标不对齐等问题。商业软件一般对单张 DOM、单个图幅的裁剪支持得很好，但是对海量 DOM、大量图幅的分幅裁剪的支持不够完善，操作烦琐，效率低下。当前，测绘行业往往采用国外商业软件进行自动化、批量化的 DOM 分幅裁剪，例如，文献 7 9利用ArcGIS 的模型构建器制作批量裁剪工具;文献 10 利用 Java 和 ArcGIS 的二次开发组件库 ArcObject 开发裁剪工具;文献 11 利用 Arcpy 脚本修改 ArcGIS 裁剪函数模板，将模板应用到栅格数据裁剪。但国外商业软件价格昂贵，不利于测绘软件的自主可控。随着开源软件的兴起，借助开源软件进行海量DOM 的自动化分幅裁剪的难度大大降低。本文提出了基于开源软件 GDAL12(Geospatial Data AbstractionLibrary)对 DOM 进行分幅裁剪的方法，通过在各步骤中综合使用 GDAL 命令和开源软件编程，实现了海量DOM 裁剪的自动化、批量化和高效化。贵阳市航空航天摄影项目验证了本文方法的有效性。2GDAL 简介GDAL 是 GIS 和遥感领域对栅格和矢量空间数据的读写支持得最好的开源软件之一12，提供了功能强大的处理命令可用于 DOM 裁剪。但大部分 GDAL 命令一般针对单个文件进行处理，在利用多个图幅裁剪时需要为各个图幅单独编写命令，在图幅多时十分不便。此外，GDAL 命令对于中文的支持不好。因此，海量 DOM的批量分幅裁剪需要引入 GDAL 的二次开发编程。除 C/C+的二次开发接口外，GDAL 还提供了 Py-thon，Java，C#等的调用接口12。在 C+接口中，对栅格/矢量数据格式的支持是通过类 GDALDriver 实现的，栅格/矢量数据集抽象为类 GDALDataset，栅格数据集的单个波段抽象为类 GDALasterBand，空间坐标系抽象为类 OGSpatialeference，GDALOpen 和 GDAL-Close 函数用于打开和关闭数据集，CPLMalloc 和 CPL-第 1 期佘佐明，程晓光，王艳军.基于 GDAL 开源库的海量 DOM 分幅裁剪Free 常用于申请和释放内存。3方法3.1源 DOM 预处理由于 GDAL 几乎支持所有的主流栅格数据格式，且 GDALDataset 屏蔽了格式差异，使用相同的方法即可读取绝大部分格式的源 DOM12，因此一般不需要考虑文件格式统一的问题。源 DOM 只需要做到数据类型、波段、坐标系、空间分辨率和无数据值的统一及像素坐标对齐。(1)数据类型统一DOM 常用 1 字节无符号整型(下文简称为 UInt8_t)存储单个波段的 DN 值。但是，DOM 的数据类型不限于此，例如，卫星遥感影像常采用 2 字节无符号整型。在 gdal_translate 命令中，可通过 ot 选项设置输出的数据类型，通过 scale 和 exponent 选项采用幂函数将位于指定范围的 DN 值拉伸到指定的输出范围12。gdal_trans-late 对应的 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 坐标系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)无数据值统一为表示某像素是否有效，测绘领域常用的方法是将无效像素的 DN 值设为无数据(NoData)值。gdalwarp 命令的 dstnodata 选项和 gdal_translate 的 a_nodata 选项可用于设置输出影像的无数据值，既可以所有波段共用相同的无数据值，也可以逐波段设置12。但该方法只是设置影像的无数据值属性，并不能修改像素 DN 值。针对 DOM 单个波段具有多个无数据值以及影像存在伪无数据区域的问题，以数据类型为 UInt8_t 的DOM 为例，本文给出一种解决方法，可统一无数据值，并排除伪无数据区域，步骤如下:将所有可能的无数据值放入无重复集合 Snodata;创建一个单通道、行列数和源 DOM 相同的二值图像 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 定义的像素坐标对齐是确保图像角点 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=xround(xulxbase)/x+xbaseyul=yround(yulybase)/y+ybase(1)3.2图幅生成文献 13给出了一种基于测区边界进行地图分幅的方法，可以自动、快速地生成位于测区内部和边界上的图幅。3.3数据部署鉴于数据量大，DOM 可能被存储在多台电脑的外存储器上。虽然可通过网络读写 DOM，但速度受带宽制约。显然，读写本地数据、在多台电脑上并行裁剪的方式效率高，在单台电脑上应尽量使用只覆盖本地源DOM 的图幅进行裁剪。此处给出一种方法:(1)使用 gdaltindex12 或其他工具生成含有本地源 DOM 覆盖范围的矢量文件;(2)将前一步生成文件与图幅接合表叠加显示;(3)将完全落在(需要考虑到图幅外扩)本地源DOM 覆盖范围内的图幅单独保存为一个图幅接合表，用于本地 DOM 裁剪;(4)将源 DOM 分布在多台电脑上的图幅保存为图幅接合表，通过网络裁剪。3.4分幅裁剪DOM 分幅裁剪的流程如图 2 所示。如果求得所有源 DOM 包容盒的并集并生成对应的矩形图幅，则该流程得到的就是拼接后的图幅。图 2影像裁剪流程图(1)重叠集定义DOM 分幅裁剪的核心问题是计算图幅对应的目标 DOM 的像素来自于哪些源 DOM 的哪些行和列，本文称为重叠集。重叠集抽象出的结构体 OverlapSet 的示例定义代码如下:struct OverlapSetstring strSrcPath;/源 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，min、yt，max分别是图幅的最小 X、最大Y 坐标，xs，min、ys，max分别是源 DOM 的最小 X、最大 Y 坐标，floor 为向下取整函数。xd，min=xroundxt，minxs，min()/x()+xs，minyd，max=yroundyt，maxys，max()/y()+ys，max(2)xd，min=xfloorxt，minxs，min()/x()+xs，minyd，max=yflooryt，maxys，max()/y()+ys，max(3)Eo需要对齐到 DOM 像素。设 Eo的最小 X、最大X、最小 Y、最大 Y 坐标分别为 xo，min、xo，max、yo，min、yo，max，无重叠法采用式(4)对齐，有重叠法采用式(5)对齐，xo，min、xo，max、yo，min、yo，max为对齐后的坐标。xo，min=xroundxo，minxs，min()/x()+xs，minxo，max=xroundxo，maxxs，min()/x()+xs，minyo，min=yroundyo，minys，max()/y()+ys，maxyo，max=yroundyo，maxys，max()/y()+ys，max(4)xo，min=xfloorxo，minxs，min()/x()+xs，minxo，max=xceilxo，maxxs，min()/x()+xs，minyo，min=yceilyo，minys，max()/y()+ys，maxyo，max=yflooryo，maxys，max()/y()+ys，max(5)67第 1 期佘佐明，程晓光，王艳军.基于 GDAL 开源库的海量 DOM 分幅裁剪重叠集属性 nSrcow0(s，0)、nSrcCol0(Cs，0)、nowNum(ol)、nColNum(Col)的计算方法见式(6)，其中 int 为浮点数取整函数，min 为计算最小值函数，s和 Cs分别是源 DOM 的行数和列数。s，0=intyo，maxys，max()/y()Cs，0=intxo，minxs，min()/x()ol=min intyo，minyo，max()/y()，ss，0()Col=min intxo，maxxo，min()/x()，CsCs，0()(6)重叠集属性 nDstow0(d，0)和 nDstCol0(Cd，0)的计算方法见式(7)。d，0=intyo，maxyd，max()/y()Cd，0=intxo，minxd，min()/x()(7)(3)基于矢量边界的裁剪有时 DOM 裁剪要求只保留测区或图幅边界内的像素，将位于边界外的像素设为无数据值。DOM 一般像素很多，逐像素判断是否在矢量边界内的效率很低。针对闭合、不自相交的矢量边界，此处给出一种高效的基于四叉树的裁剪方法，步骤包括:对影像进行迭代四叉树分割，除叶节点外，其他节点均含 4 棵子树，叶节点包含的像素数目大于等于 1;生成闭合矢量边界对应的多边形(Multipoly-gon)对象 Pbound;以四叉树根节点为当前节点，进入步骤;设当前节点为 Tnode，计算 Tnode对应的二维包容盒 Enode，生成 Enode对应的多边形(Polygon)对象 Pnode，计算 Pnode的面积 Anode、Pbound与 Pnode的交集 Pinter、Pinter的面积 Ainter。如果 Ainter=Anode，说明 Tnode的所有像素均在矢量边界内;如果 Ainter=0，说明 Tnode的所有像素均在矢量边界外;如果 Ainter 0，Anode()，说明 Tnode部分像素在矢量边界内，部分像素在矢量边界外。对最后一种情况，如果 Tnode为叶节点，则遍历它包含的所有像素，判断像素中心是否在 Pbound内;如果 Tnode含 4 个子节点，则将 Tnode的每个子节点 Tchild，kk 1，2，3，4()作为当前节点进入步骤迭代计算。对不在矢量边界内的像素，将其各波段 DN 值设为指定的无数据值。(4)DOM 读写GDALDriver 支持使用 Create 或 CreateCopy 函数创建栅格数据文件12。鉴于源 DOM 和目标 DOM 之间一般不存在精确的一对一关系，所以只能采用 Create 函数创建目标 DOM。GDALDriver 允许在 Create 函数中设置创建文件选项，如对 GeoTiff，设置 TFW=YES 以生成对应的 TFW 文件、TILED=YES 以按瓦片存储等。在实际写入波段 DN 值数据前，可调用 GDALas-terBand 的 Fill 函数，用该波段的无数据值进行填充。GDALDataset 和 GDALasterBand 均提供了具有普适性的 asterIO 函数用于数据读写。但是，要想在全分辨率(不涉及重采样)的数据读写上取得最好的性能，GDAL推荐使用 GDALasterBand 的 eadBlock 和 WriteBlock函数读写数据块12。此外，使用这两个函数也可尽量避免因 DOM 数据量大导致的内存分配失败。超大 DOM 的像素或 DN 值数目可能超过 4 字节整型可表达的最大值，应使用 8 字节整型。3.5结果检查在裁 剪 初 步 完 成 后，需 要 对 结 果 进 行 检 查。GDALDataset 的 GetSpatialef、GetasterCount、GetGeo-Transform 函数可用于坐标系、波段数目、分辨率的检查，GDALasterBand 的 GetasterDataType、GetNoData-Value 函数可用于数据类型、无数据值的检查12。此处给出一种判断裁剪后的 DOM 之间是否存在缝隙和重叠的方法:(1)计算所有裁剪后的 DOM 的二维包容盒;(2)对各个二维包容盒略微外扩(适用于缝隙检查)或内缩(适用于重叠检查);(3)对所有二维包容盒求空间并集(适用于缝隙检查)或两两求空间交集(适用于重叠检查)，将结果保存为矢量文件;(4)打开矢量文件，目视检查是否存在缝隙或重叠。该方法引入包容盒内缩和外扩的目的是避免包容盒在边界处重叠导致的交并集计算误差。内缩或外扩的距离取 x的 1%10%即可。4方法验证本文利用 C+和 GDAL 研发了 DOM 分幅裁剪软件，并应用在 2020 年度贵阳市航空航天摄影项目中。该项目总航测面积为 8 525 km2，源 DOM 为 02 m分辨率、UInt8_t 数据类型、按行政区划存储的 Erdas Imagine格式的真彩色影像。对贵阳 6 个市辖区和双龙新区的DOM(文件总大小约为 435 GB)采用 1 2 000比例尺的正方形图幅进行了裁剪。在 CPU 为 Intel Core i97940X、主频 31 GHz的单台服务器上使用测区矢量边界和 9 个线程并行裁剪为 3 323张图幅。各线程运行时间之和约为 12 h，生成的 GeoTiff 格式的分幅影像、金字塔和 TFW 文件总大小约为 301 GB，单个线程使用内存不超过 2 GB。分幅裁剪成果样例如图 3 所示。77城市勘测2023 年 2 月图 3DOM 分幅裁剪成果样例5结论本文提出了基于 GDAL 的海量 DOM 分幅裁剪方法，通过综合使用 GDAL 命令和开源软件编程，在最大程度上提高了 DOM 裁剪的效率和自动化、批量化水平。同时，本文提出的方法还适用于 DEM、DSM 等其他栅格数据。参考文献 1GBT 139902012 国家基本比例尺地形图分幅和编号 S 2任宇征，韩磊，董敏 ArcGIS 和 FME 在权籍调查工作中若干应用对比分析 J 测绘与空间地理信息，2021，44(10):113115 3白洪伟，徐洋洋 基于 ENVI 和 ArcGIS 的合肥市土地利用/覆被变化分析J 沈阳大学学报自然科学版，2016，28(5):365372 4宋伟东，赵芳杰，赵开会 基于遥感影像的交通专题图制作 J 测绘与空间地理信息，2018，41(10):1214+17 5朱典潞 FME 在 GIS 数据处理中的应用 J 城市勘测，2021(4):6063+67 6董正国，丁云 ArcGIS 影像服务在国情监测中的应用 J 测绘与空间地理信息，2021，44(6):129131+134 7周春峰 浅谈基于 ArcGIS 的影像批量裁剪J 北京测绘，2015(6):136137+146 8昝建春，冷鸿天，吴颜奎 基于 GIS 模型构建器的影像批量裁剪方法研究 J 林业调查规划，2022，47(2):1619+64 9黄国森，阮芳 基于镶嵌数据集的多源影像标准分幅裁剪 J 测绘与空间地理信息，2021，44(3):124126 10 张广耀，王蕾 基于 ArcObject 的遥感影像裁剪探讨 J 中华建设，2018(11):144145 11 高磊 遥感影像动态裁剪程序化实现方法J 北京测绘，2020，34(2):172175 12Contributors GDAL/OG GDAL/OG Geospatial DataAbstraction Software LibraryDB/OL:Open Source Geo-spatial Foundation，2022 https:/gdalorg/105281/zeno-do5884351 13 程晓光，陆泉源，郑诚慧等 使用开源软件的地图分幅软件设计与实现J 测绘地理信息，2021，46(4):122127The Subdivision Clipping of Massive DOM Based on GDALOpen Source LibrarySHE Zuoming1，CHENG Xiaoguang2，WANG Yanjun3(1.Guiyang Institute of Surveying and Mapping，Guizhou 550000，China;2.Feiyan Aerial emote Sensing and Tech-nology Co，Ltd，Wuhan 430061，China;3.Nationallocal Joint Engineering Laboratory of Geospatial InformationTechnology，Hunan University of Science and Technology，Xiangtan 411201，China)Abstract:The subdivision clipping of massive digital orthophoto map(DOM)is a fundamental problem in the fieldof mapping and GIS The DOM clipping function of existing commercial software has the problems of low automation andlow efficiency For the steps involved in massive DOM clipping，including DOM preprocessing，subdivision，data deploy-ment，image clipping and result checking，a solution based on open source library GDAL is proposed in this paperThrough the comprehensive use of GDAL command and open source software programming，the efficiency and autonomouslevel of massive DOM subdivision clipping can be improved to the maximum extent The feasibility and effectiveness ofthe proposed solution are well verified by Guiyang aerial photogrammetry projectsKey words:massive DOM;subdivision clipping;GDAL;open source software87