1. 引言

筆者在實現柵格數據的可視化的時候遇到了一個問題，計算柵格數據金字塔層級的地理變換信息錯誤導致可視化的時候存在微小的誤差。地理變換信息指的就是柵格數據的地理座標起點和分辨率，筆者在另外一篇文章中《GDAL讀取的座標起點在像素左上角還是像素中心？》論述了柵格數據集中座標起點位置存在半個像素差的問題。但是柵格數據集的金字塔層級影像是如何處理這個問題的呢？

2. 詳論

2.1 連續還是離散

從《GDAL讀取的座標起點在像素左上角還是像素中心？》這篇文章繼續引申一個問題：柵格數據究竟是連續的還是離散的？從GIS的角度來看，柵格數據就是真實世界地理實體的表達，肯定應該是連續的。但是問題在於，現實世界並不存在如此完美的載體能夠表達連續的實體對象，大多數都會離散化為網格。比如圖像、屏幕這些數據載體，它們本質上就是柵格，你將它們放大來看，就能看到一個個放大的格子。GIS的柵格數據就是通過這些離散的數據載體來表達的，那你能説柵格數據一定就是連續的嗎？所以筆者有一句結論：

GIS的柵格數據在宏觀上是連續的，在微觀上是離散的。

其實這個結論有的讀者不一定能夠接受，因為每個人都有自己的固有認知，強行接受別人的認知無異震碎自己的三觀，比如筆者的前同事，就對我這個理論不以為然，總會糾結一個問題：如果柵格數據都離散成整數了，那麼橫座標為0.6的像素在哪裏？其實我覺得不妨這樣理解，一個離散的像素網格確實對應了真實地理實體的一段長度，不過一個像素網格已經是最小的表達實體了，那麼就只能使用像素網格中心的位置的值來表達代表整個像素網格的值。橫座標為0.6的像素確實在數據載體中不存在，但是其表達的數據中確實客觀存在的，一定要獲取這個值也有處理方法，那就是數值內插算法，比如《數字圖像處理》中經常提到的最鄰近、雙線性以及三次卷積等內插算法。

2.2 起點位置問題

在筆者看來，柵格數據是連續的還是離散的看法其實關聯着柵格數據中地理座標起點的位置問題。tif數據內部存儲的地理變換信息規定地理座標起點在像素左上角，這是着重於柵格數據的連續性；tif外部數據tfw規定地理座標起點在像素中心點，則着重於柵格數據的離散型。那麼什麼時候用像素左上角作為起點，什麼時候用像素中心點作為起點呢？筆者的看法是，僅以GIS柵格數據來説：

進行空間計算的時候應該以像素左上角為起點；進行圖像處理的時候應該以像素中心點為起點。

舉例來説，筆者這裏有一個柵格數據dom.tif，通過gdalinfo查詢信息如下：

Driver: GTiff/GeoTIFF
Files: dom.tif
       dom.tif.aux.xml
Size is 19312, 22531
Origin = (1967768.351536701666191,570294.132588228210807)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 1967768.352,  570294.133)
Lower Left  ( 1967768.352,  564661.383)
Upper Right ( 1972596.352,  570294.133)
Lower Right ( 1972596.352,  564661.383)
Center      ( 1970182.352,  567477.758)
Band 1 Block=19312x32 Type=Byte, ColorInterp=Red
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=110.556, StdDev=57.195
  Overviews: 9656x11266, 4828x5633, 2414x2817, 1207x1409, 604x705, 302x353, 151x177
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=110.55552426626
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=57.19535628196
    STATISTICS_VALID_PERCENT=100
Band 2 ...
Band 3 ...

這裏的四至範圍（Upper Left、Lower Left、Upper Right、Lower Right）的計算就是以像素左上角為起點（Origin）來進行計算的，有興趣的讀者可以進行驗算一下。然而，如果你需要可視化這個柵格數據，通常需要將柵格值傳遞到GUI畫布的Image對象值中，將柵格像素對齊GUI像素，每個像素的值應該是其像素中心的地理座標的像素值，這時最好以像素中心作為起點進行計算。

2.3 金字塔層級影像

最後回到金字塔層級的地理變換信息計算的問題。在進行可視化的時候，使用像素中心作為起點進行空間座標的計算，重採樣出合適的像素值。但是首先，任何金字塔層級影像的四至範圍是不能有變化的。以上例來説，原始圖像的寬高是19312×22531，第一層金字塔影像是9656×11266，如果以像素左上角為起點，計算其地理變換信息：

起點位置不會變化，Origin = (1967768.351536701666191,570294.132588228210807)。
分辨率則需要重新計算，在X方向，（19312 * 0.25）/ 9656 = 0.5；在Y方向，（22531 * 0.25）/ 11266 = 0.49997780933783064。也就是Pixel Size = (0.500000000000000,-0.49997780933783064)。

這時會發生一個略顯詭異的現象，就是原始柵格影像上X方向Y方向分辨率一致都是0.25，在第一級金字塔層級影像上卻已經有微小的差異了。其實產生這個現象的原因並不奇怪，因為金字塔層級的寬高通常以2的倍數遞減，但是原始柵格的寬高卻不是偶數。在這種情況下，要優先保證地理座標的四至範圍的一致性。

接下來計算以像素中心為起點，地理變換信息：

分辨率不會發生變化，同樣是Pixel Size = (0.500000000000000,-0.49997780933783064)。
起點位置則需要偏移半個像素，移動到像素中心，Origin = (1967768.601536701666191,570293.88259932354189168)

一定要注意，在以像素中心為起點的情況下，每個金字塔層級影像的起點座標是會變化的，因為每個金字塔層級的分辨率不同。

3. 其他

通過上述方法，可以計算任意金字塔層級影像的地理變換信息。不止是金字塔層級影像，筆者在《GDAL關於讀寫圖像的簡明總結》這篇文章中介紹過GDAL讀取柵格數據的時候可以重採樣讀取，如下所示：

//申請buf
size_t imgBufNum = (size_t) bufWidth * bufHeight * bandNum * depth;
size_t imgBufOffset = (size_t) bufWidth * (bufHeight-1) * bandNum * depth;
GByte *imgBuf = new GByte[imgBufNum];	
//讀取
img->RasterIO(GF_Read, 0, 0, imgWidth, imgHeight, imgBuf + imgBufOffset, bufWidth, bufHeight,
	GDT_Byte, bandNum, nullptr, bandNum*depth, -bufWidth*bandNum*depth, depth);
//寫入
dst->RasterIO(GF_Write, 0, 0, bufWidth, bufHeight, imgBuf + imgBufOffset, bufWidth, bufHeight,
	GDT_Byte, bandNum, nullptr, bandNum*depth, -bufWidth*bandNum*depth, depth);
//釋放
delete[] imgBuf;
imgBuf = nullptr;

這裏讀取的影像塊寬高imgWidth、imgHeight與buffer塊寬高bufWidth、bufHeight是可以不一致的，並且GDAL會通過重採樣自動處理這個過程。影像塊的地理變換信息我們很清楚，但是buffer塊地理變換信息則需要自己進行計算。原理也是一樣的，要優先保證buffer塊的四至範圍與影像塊的四至範圍一致，先算分辨率，在偏移半個像素的距離算出具體的起點位置。