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

1. 引言

笔者在实现栅格数据的可视化的时间遇到了一个题目，计算栅格数据金字塔层级的地理变换信息错误导致可视化的时间存在微小的误差。地理变换信息指的就是栅格数据的地理坐标起点和分辨率，笔者在另外一篇文章中《GDAL读取的坐标起点在像素左上角还是像素中心？》论述了栅格数据会合坐标起点位置存在半个像素差的题目。但是栅格数据集的金字塔层级影像是如何处理这个题目的呢？
2. 详论

2.1 连续还是离散

从《GDAL读取的坐标起点在像素左上角还是像素中心？》这篇文章继续引申一个题目：栅格数据究竟是连续的还是离散的？从GIS的角度来看，栅格数据就是真实世界地理实体的表达，肯定应该是连续的。但是题目在于，现实世界并不存在如此完善的载体能够表达连续的实体对象，大多数都会离散化为网格。比如图像、屏幕这些数据载体，它们本质上就是栅格，你将它们放大来看，就能看到一个个放大的格子。GIS的栅格数据就是通过这些离散的数据载体来表达的，那你能说栅格数据一定就是连续的吗？所以笔者有一句结论：

GIS的栅格数据在宏观上是连续的，在微观上是离散的。

实在这个结论有的读者不一定能够接受，因为每个人都有自己的固有认知，强行接受别人的认知无异震碎自己的三观，比如笔者的前同事，就对我这个理论不以为然，总会纠结一个题目：如果栅格数据都离散成整数了，那么横坐标为0.6的像素在哪里？实在我觉得不妨如许明确，一个离散的像素网格确实对应了真实地理实体的一段长度，不过一个像素网格已经是最小的表达实体了，那么就只能使用像素网格中心的位置的值来表达代表整个像素网格的值。横坐标为0.6的像素确实在数据载体中不存在，但是其表达的数据中确实客观存在的，一定要获取这个值也有处理方法，那就是数值内插算法，比如《数字图像处理》中常常提到的最邻近、双线性以及三次卷积等内插算法。
2.2 起点位置题目

在笔者看来，栅格数据是连续的还是离散的看法实在关联着栅格数据中地理坐标起点的位置题目。tif数据内部存储的地理变换信息规定地理坐标起点在像素左上角，这是着重于栅格数据的连续性；tif外部数据tfw规定地理坐标起点在像素中心点，则着重于栅格数据的离散型。那么什么时间用像素左上角作为起点，什么时间用像素中心点作为起点呢？笔者的看法是，仅以GIS栅格数据来说：

举行空间计算的时间应该以像素左上角为起点；举行图像处理的时间应该以像素中心点为起点。

举例来说，笔者这里有一个栅格数据dom.tif，通过gdalinfo查询信息如下：

1. Driver: GTiff/GeoTIFF
2. Files: dom.tif
3.        dom.tif.aux.xml
4. Size is 19312, 22531
5. Origin = (1967768.351536701666191,570294.132588228210807)
6. Pixel Size = (0.250000000000000,-0.250000000000000)
7. Metadata:
8.   AREA_OR_POINT=Area
9. Image Structure Metadata:
10.   COMPRESSION=LZW
11.   INTERLEAVE=PIXEL
12. Corner Coordinates:
13. Upper Left  ( 1967768.352,  570294.133)
14. Lower Left  ( 1967768.352,  564661.383)
15. Upper Right ( 1972596.352,  570294.133)
16. Lower Right ( 1972596.352,  564661.383)
17. Center      ( 1970182.352,  567477.758)
18. Band 1 Block=19312x32 Type=Byte, ColorInterp=Red
19.   Min=0.000 Max=255.000
20.   Minimum=0.000, Maximum=255.000, Mean=110.556, StdDev=57.195
21.   Overviews: 9656x11266, 4828x5633, 2414x2817, 1207x1409, 604x705, 302x353, 151x177
22.   Metadata:
23.     STATISTICS_MAXIMUM=255
24.     STATISTICS_MEAN=110.55552426626
25.     STATISTICS_MINIMUM=0
26.     STATISTICS_STDDEV=57.19535628196
27.     STATISTICS_VALID_PERCENT=100
28. Band 2 ...
29. 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读取栅格数据的时间可以重采样读取，如下所示：

1. //申请buf
2. size_t imgBufNum = (size_t) bufWidth * bufHeight * bandNum * depth;
3. size_t imgBufOffset = (size_t) bufWidth * (bufHeight-1) * bandNum * depth;
4. GByte *imgBuf = new GByte[imgBufNum];       
5. //读取
6. img->RasterIO(GF_Read, 0, 0, imgWidth, imgHeight, imgBuf + imgBufOffset, bufWidth, bufHeight,
7.         GDT_Byte, bandNum, nullptr, bandNum*depth, -bufWidth*bandNum*depth, depth);
8. //写入
9. dst->RasterIO(GF_Write, 0, 0, bufWidth, bufHeight, imgBuf + imgBufOffset, bufWidth, bufHeight,
10.         GDT_Byte, bandNum, nullptr, bandNum*depth, -bufWidth*bandNum*depth, depth);
11. //释放
12. delete[] imgBuf;
13. imgBuf = nullptr;

复制代码

这里读取的影像块宽高imgWidth、imgHeight与buffer块宽高bufWidth、bufHeight是可以不同等的，而且GDAL会通过重采样主动处理这个过程。影像块的地理变换信息我们很清晰，但是buffer块地理变换信息则需要自己举行计算。原理也是一样的，要优先包管buffer块的四至范围与影像块的四至范围同等，先算分辨率，在偏移半个像素的距离算出具体的起点位置。

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。