栅格數據有效/無效掩膜的(de)計算(suàn)與表達

因爲多(duō)種原因，栅格數據中經常存在一些無效像素。在航天測繪遙感中，遙感器獲取的(de)單景影(yǐng)像在經過投影(yǐng)坐(zuò)标變換後往往會存在黑(hēi)邊（見圖 1(a)），黑(hēi)邊内的(de)像素爲無效像素。在航空測繪遙感中，位于測區(qū)外的(de)像素也(yě)可(kě)能被拍(pāi)攝/掃描到，但是在成果交付時(shí)往往隻需要輸出測區(qū)内的(de)像素，測區(qū)外的(de)像素即爲無效像素（見圖 2(a)）。栅格數據常采用(yòng)二值的(de)掩膜來(lái)表達像素有效或無效。

基于像素值的(de)掩膜計算(suàn)

單個(gè)波段的(de)掩膜的(de)計算(suàn)非常簡單。設對(duì)于某個(gè)波段，值爲0代表無效值，值非0代表有效值，則利用(yòng)Envi的(de)band math功能計算(suàn)掩膜的(de)方法爲：

byte((b1 ne 0)*255)

設值範圍在[M,N]内的(de)爲有效值，其它值爲無效值，則計算(suàn)掩膜的(de)方法爲：

byte((b1 ge M and b1 le N)*255)

基于各個(gè)波段的(de)掩膜，可(kě)以計算(suàn)出整個(gè)數據集的(de)掩膜。對(duì)于3波段真彩色影(yǐng)像，計算(suàn)方法爲：

byte((b1 eq 255) and (b2 eq 255) and (b3 eq 255))*byte(255)

在生成的(de)掩膜中，255代表有效像素，0代表無效像素。相關示例可(kě)見圖 1。

圖 1 基于像素值的(de)掩膜計算(suàn)示例

基于矢量邊界的(de)掩膜計算(suàn)

無論是基于測區(qū)邊界，圖幅邊界還(hái)是ROI生成掩膜，均可(kě)以歸結爲基于矢量邊界的(de)掩膜計算(suàn)。如果矢量邊界對(duì)應的(de)面狀幾何元素存在多(duō)個(gè)不重疊的(de)部分(fēn)，例如測區(qū)由多(duō)個(gè)不聯通(tōng)、不重疊的(de)子測區(qū)構成[1]，則矢量邊界可(kě)以表達爲OGC簡單特征規範定義的(de)多(duō)多(duō)邊形（MultiPolygon），否則可(kě)以表達爲多(duō)邊形（Polygon），多(duō)邊形允許帶有空洞[2]。

基于矢量邊界計算(suàn)掩膜，可(kě)以采用(yòng)矢量到栅格轉換的(de)方法（需要确保兩種栅格數據的(de)分(fēn)辨率一緻），也(yě)可(kě)以采用(yòng)幾何分(fēn)析的(de)方法。下(xià)面介紹幾何分(fēn)析的(de)方法。

文獻[3]提出一種基于四叉樹的(de)方法，其核心思想是對(duì)栅格數據進行四叉樹分(fēn)割，判斷分(fēn)割出的(de)每個(gè)部分(fēn)對(duì)應的(de)矩形和(hé)矢量邊界對(duì)應的(de)面狀幾何元素之間是否存在重疊。如果存在重疊，則對(duì)矩形再次進行四叉樹分(fēn)割，直至分(fēn)割出的(de)每個(gè)部分(fēn)所含像素數目低于指定阈值。由于該算(suàn)法是遞歸執行的(de)，在栅格數據行列數較多(duō)時(shí)，對(duì)棧空間的(de)需求較大(dà)，容易發生棧溢出。爲提高(gāo)計算(suàn)穩定性，可(kě)以逐塊（Block）計算(suàn)。對(duì)于支持按塊存儲的(de)格式，例如Tiff，按塊存儲栅格數據可(kě)以有效降低計算(suàn)對(duì)棧空間的(de)需求。設每塊大(dà)小爲256*256，則四叉樹最大(dà)深度爲8層。對(duì)于Asc格式，雖然其不支持按塊存儲，但該方法仍然可(kě)用(yòng)，隻不過一行爲一個(gè)塊，四叉樹分(fēn)割變爲二叉樹分(fēn)割。

判斷單個(gè)像素是否在矢量邊界内，可(kě)采用(yòng)如下(xià)兩種判斷方法：

(1)基于像素中心的(de)方法。如果像素中心在矢量邊界對(duì)應的(de)面狀幾何元素内，則認爲像素爲有效像素，否則爲無效像素；

(2)基于像素邊界的(de)方法。計算(suàn)像素邊界對(duì)應的(de)矩形和(hé)矢量邊界對(duì)應的(de)面狀幾何元素之間的(de)重疊面積。如果重疊面積占像素面積的(de)比例大(dà)于指定阈值p，則認爲像素爲有效像素，否則爲無效像素。其中，如果p取0，則意味著(zhe)如果像素有任何一部分(fēn)在矢量邊界内，像素即爲有效像素；如果p取0.5，且在單個(gè)像素内的(de)矢量邊界爲直線段，則此時(shí)該判斷方法和(hé)基于像素中心的(de)方法等價。

圖 2是利用(yòng)矢量邊界裁剪DEM數據的(de)效果圖。

圖 2 基于矢量邊界的(de)裁剪DEM的(de)效果，紅線爲矢量邊界

掩膜的(de)表達形式

掩膜的(de)主要表達形式包括：

(1)掩膜文件。GDAL采用(yòng)後綴爲msk，但實質爲TIFF格式的(de)栅格文件（簡稱掩膜文件）存儲掩膜數據。掩膜文件的(de)行列數、塊大(dà)小、路徑與對(duì)應的(de)栅格數據相同，但文件名比栅格數據多(duō)“.msk”。掩膜文件存儲的(de)數據類型爲1字節無符号整型（下(xià)文簡稱uint_8t），0代表像素無效，非0（一般是255）代表像素有效。掩膜文件的(de)波段數要麽爲1，适用(yòng)于整個(gè)栅格數據，要麽和(hé)栅格數據的(de)波段數相同，第K個(gè)波段适用(yòng)于栅格數據的(de)第K個(gè)波段。

(2)無數據值（nodata value）。在GeoTiff、Asc、ERDAS IMAGINE等格式中，支持無數據值屬性。有的(de)格式支持對(duì)單個(gè)波段設置無數據值（如GeoTiff），有的(de)格式隻支持對(duì)數據集整體設置無數據值（如Asc）。常見的(de)GIS和(hé)遙感軟件均支持該屬性（見圖 3），在讀取波段數據時(shí)，如果某像素在某波段的(de)值爲指定的(de)無數據值，則認爲該像素在該波段無效。在使用(yòng)該方式時(shí)，需要将無效像素的(de)值設置爲無數據值。在DEM/DSM中，通(tōng)常選擇一個(gè)明(míng)顯不合理(lǐ)的(de)高(gāo)程值作爲無數據值，例如-9999 m。在數據類型爲uint_8t的(de)3波段真彩色DOM中，通(tōng)常選擇0或255作爲無數據值。但是該做(zuò)法的(de)風險在于，DOM中的(de)有效像素的(de)波段值也(yě)可(kě)能爲0或255。例如，陰影(yǐng)、水(shuǐ)體等低光(guāng)部分(fēn)在某些波段可(kě)能值爲0（見圖 4(a)），白車車頂、雲等高(gāo)光(guāng)部分(fēn)在某些波段可(kě)能值爲255。所以，爲确保有效像素不會出現無數據值，有可(kě)能需要修改像素值[3]。

圖 3 常見GIS遙感軟件對(duì)無數據值的(de)支持

(3)alpha波段。在數據類型爲uint_8t的(de)影(yǐng)像中，alpha波段常被用(yòng)來(lái)存儲不透明(míng)度（0代表完全透明(míng)，255代表完全不透明(míng)）。當使用(yòng)alpha波段存儲掩膜時(shí)，GDAL規定無效像素表示爲0值，有效像素表示爲非0值（效果見圖 4(b)）。該方法在非遙感GIS領域内應用(yòng)較多(duō)。但缺點在于，一來(lái)隻能用(yòng)于數據類型爲uint_8t的(de)真彩色影(yǐng)像或灰度影(yǐng)像；二來(lái)會變相增加1個(gè)波段，占用(yòng)較多(duō)的(de)存儲空間。

(4)全部像素均有效。在無特殊标記/說明(míng)的(de)情況下(xià)，全部像素均爲有效像素。

在GDAL中，按(1)、(2)、(3)、(4)的(de)優先順序識别掩膜。

圖 4 采用(yòng)不同方式的(de)掩膜的(de)效果

參考文獻

[1]程曉光(guāng),魏婧,嚴明(míng),等. 應用(yòng)地表覆蓋數據估算(suàn)LiDAR内業工作量的(de)方法研究[J]. 地理(lǐ)信息世界, 2021, 28(2): 89-94.

[2]程曉光(guāng),陸泉源,鄭誠慧,等. 使用(yòng)開源軟件的(de)地圖分(fēn)幅軟件設計與實現[J]. 測繪地理(lǐ)信息, 2021, 46(4): 122-127.

[3]佘佐明(míng),程曉光(guāng),王豔軍. 基于GDAL開源庫的(de)海量DOM分(fēn)幅裁剪[J]. 城(chéng)市勘測, 2023, No.195(1): 74-78.