柵格數據有效/無(wú)效掩膜的計算與表達

因為多種原因，柵格數據中經(jīng)常存在一些無(wú)效像素。在航天測繪遙感中，遙感器獲取的單景影像在經(jīng)過(guò)投影坐標變換后往往會(huì )存在黑邊（見(jiàn)圖 1(a)），黑邊內的像素為無(wú)效像素。在航空測繪遙感中，位于測區外的像素也可能被拍攝/掃描到，但是在成果交付時(shí)往往只需要輸出測區內的像素，測區外的像素即為無(wú)效像素（見(jiàn)圖 2(a)）。柵格數據常采用二值的掩膜來(lái)表達像素有效或無(wú)效。

基于像素值的掩膜計算

單個(gè)波段的掩膜的計算非常簡(jiǎn)單。設對于某個(gè)波段，值為0代表無(wú)效值，值非0代表有效值，則利用Envi的band math功能計算掩膜的方法為：

byte((b1 ne 0)*255)

設值范圍在[M,N]內的為有效值，其它值為無(wú)效值，則計算掩膜的方法為：

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

基于各個(gè)波段的掩膜，可以計算出整個(gè)數據集的掩膜。對于3波段真彩色影像，計算方法為：

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

在生成的掩膜中，255代表有效像素，0代表無(wú)效像素。相關(guān)示例可見(jiàn)圖 1。

圖 1 基于像素值的掩膜計算示例

基于矢量邊界的掩膜計算

無(wú)論是基于測區邊界，圖幅邊界還是ROI生成掩膜，均可以歸結為基于矢量邊界的掩膜計算。如果矢量邊界對應的面狀幾何元素存在多個(gè)不重疊的部分，例如測區由多個(gè)不聯(lián)通、不重疊的子測區構成[1]，則矢量邊界可以表達為OGC簡(jiǎn)單特征規范定義的多多邊形（MultiPolygon），否則可以表達為多邊形（Polygon），多邊形允許帶有空洞[2]。

基于矢量邊界計算掩膜，可以采用矢量到柵格轉換的方法（需要確保兩種柵格數據的分辨率一致），也可以采用幾何分析的方法。下面介紹幾何分析的方法。

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

判斷單個(gè)像素是否在矢量邊界內，可采用如下兩種判斷方法：

(1)基于像素中心的方法。如果像素中心在矢量邊界對應的面狀幾何元素內，則認為像素為有效像素，否則為無(wú)效像素；

(2)基于像素邊界的方法。計算像素邊界對應的矩形和矢量邊界對應的面狀幾何元素之間的重疊面積。如果重疊面積占像素面積的比例大于指定閾值p，則認為像素為有效像素，否則為無(wú)效像素。其中，如果p取0，則意味著(zhù)如果像素有任何一部分在矢量邊界內，像素即為有效像素；如果p取0.5，且在單個(gè)像素內的矢量邊界為直線(xiàn)段，則此時(shí)該判斷方法和基于像素中心的方法等價(jià)。

圖 2是利用矢量邊界裁剪DEM數據的效果圖。

圖 2 基于矢量邊界的裁剪DEM的效果，紅線(xiàn)為矢量邊界

掩膜的表達形式

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

(1)掩膜文件。GDAL采用后綴為msk，但實(shí)質(zhì)為T(mén)IFF格式的柵格文件（簡(jiǎn)稱(chēng)掩膜文件）存儲掩膜數據。掩膜文件的行列數、塊大小、路徑與對應的柵格數據相同，但文件名比柵格數據多“.msk”。掩膜文件存儲的數據類(lèi)型為1字節無(wú)符號整型（下文簡(jiǎn)稱(chēng)uint_8t），0代表像素無(wú)效，非0（一般是255）代表像素有效。掩膜文件的波段數要么為1，適用于整個(gè)柵格數據，要么和柵格數據的波段數相同，第K個(gè)波段適用于柵格數據的第K個(gè)波段。

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

圖 3 常見(jiàn)GIS遙感軟件對無(wú)數據值的支持

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

(4)全部像素均有效。在無(wú)特殊標記/說(shuō)明的情況下，全部像素均為有效像素。

在GDAL中，按(1)、(2)、(3)、(4)的優(yōu)先順序識別掩膜。

圖 4 采用不同方式的掩膜的效果

參考文獻

[1]程曉光,魏婧,嚴明,等. 應用地表覆蓋數據估算LiDAR內業(yè)工作量的方法研究[J]. 地理信息世界, 2021, 28(2): 89-94.

[2]程曉光,陸泉源,鄭誠慧,等. 使用開(kāi)源軟件的地圖分幅軟件設計與實(shí)現[J]. 測繪地理信息, 2021, 46(4): 122-127.

[3]佘佐明,程曉光,王艷軍. 基于GDAL開(kāi)源庫的海量DOM分幅裁剪[J]. 城市勘測, 2023, No.195(1): 74-78.