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塊的四至范圍與影像塊的四至范圍一致，先算分辨率，在偏移半個像素的距離算出具體的起點位置。