栅格数据金字塔层级的地理变换信息

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.tifdom.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=LZWINTERLEAVE=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=RedMin=0.000 Max=255.000Minimum=0.000, Maximum=255.000, Mean=110.556, StdDev=57.195Overviews: 9656x11266, 4828x5633, 2414x2817, 1207x1409, 604x705, 302x353, 151x177Metadata:STATISTICS_MAXIMUM=255STATISTICS_MEAN=110.55552426626STATISTICS_MINIMUM=0STATISTICS_STDDEV=57.19535628196STATISTICS_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块的四至范围与影像块的四至范围一致，先算分辨率，在偏移半个像素的距离算出具体的起点位置。