COG（Cloud optimized GeoTIFF——云优化GeoTiff）简介与实践

参考
简介部分翻译自cogeo
关于GeoTIFF文件格式，参考了osgeo和awaresystems
关于HTTP请求范围，参考MDN-HTTP范围请求

云优化GeoTIFF(COG)是一个常规的GeoTIFF文件，旨在托管在 HTTP 文件服务器上，其内部组织可以在云上实现更高效的工作流。它通过利用客户端发出HTTP来仅请求他们需要的文件部分来做到这一点。

常见的问题：

1. 不是很明白这个东西有啥用？
   简单来说，可以代替geoserver/arcgis server发布栅格数据，比wms/wmts节省服务器空间和内存
2. 支持多大的tif？
   我试过最大4G多的tif，一点问题没有，缩放平移十分流畅
3. 缺点和局限？
   目前前端API里只知道OpenLayers支持COG直接渲染，其他api不清楚。
   OL使用webgl渲染cog，对于用canvas渲染的图层适用的比如后处理功能，有些不适用。
   需要对数据做转化处理。

技术概述

COG依赖两种技术

• 在GeoTIFF图片中保存除原始像元数据外的瓦片（tile）和概览（overview）
• 使用HTTP范围请求获取单一文件的所需部分

GeoTIFF构建

COG所需的GeoTIFF，主要技术是在文件内部构建瓦片和概览，辅以图像压缩。

关于瓦片

熟悉地图开发的小伙伴应该对瓦片都很熟悉了，不清楚的可以自己搜一下关键字地图瓦片，简单介绍就是将完整的影像像拼图一样按照横竖的格网拆分成瓦片，只加载需要范围内的一个个瓦片而不是完整的巨大TIFF。
常用的影像瓦片技术，比如ArcMap/ArcGIS Server，GeoServer的切片，都需要预先生成好切片存储到磁盘上，而COG是以瓦片格式存到TIFF文件中。

关于概览

概览创建同一图像的低分辨率重采样版本。这意味着它的细节要少得多（原始图像可能有 100 或 1000 像素对应到概览的 1 个像素），但也小得多。通常一个 GeoTIFF 会有很多概览，以匹配不同的缩放级别。这些增加了整个文件的大小，但能够更快地提供服务。
用过ArcMap的可能更熟悉影像金字塔这个概念，和同事讨论了下，影像金字塔应该和概览意思差不多。

完整的文件结构

• TIFF / BigTIFF 签名
• 全分辨率图像的 IFD ( ​Image File Directory文件图像目录 )
• 不适合内嵌在 IFD 目录中的 TIFF 标签的值，例如TileOffsets，TileByteCounts和GeoTIFF键
• 可选：第一个概览的 IFD（通常以 2 倍子采样），后跟不适合内联的标签值
• 可选：第二个概览的 IFD（图像文件目录）（通常以 4 倍子采样），后跟不适合内联的标签值
• ···
• 可选：最后概览的 IFD（通常以 2的N次幂进行二次采样），后跟不适合内联的标签值
• 可选：最后一个概览级别的瓦片内容
• ···
• 可选：第一个概览级别的瓦片内容
• 全分辨率图像的瓦片内容

HTTP范围请求

HTTP 协议范围请求允许服务器只发送 HTTP 消息的一部分到客户端。范围请求在传送大的媒体文件，或者与文件下载的断点续传功能搭配使用时非常有用。

检测服务器端是否支持范围请求

假如在响应中存在 Accept-Ranges 头部（并且它的值不为 “none”），那么表示该服务器支持范围请求。

这里请求的是我本地Nginx里的COG图像。
在上面的响应中， Accept-Ranges: bytes 表示界定范围的单位是 bytes 。这里 Content-Length 也是有效信息，因为它提供了要检索的图片的完整大小。

从服务器端请求特定的范围

假如服务器支持范围请求的话，你可以使用 Range 头来生成该类请求。该首部指示服务器应该返回文件的哪一或哪几部分。

这里就请求的该TIFF的0-65536字节。

两者结合

GeoTIFF将瓦片和概览按规则进行组织，存放放在云端（比如我本地测试用的Nginx），以便HTTP范围请求可以只请求相关的文件部分。
当客户端想要渲染整个文件的概览图像时，它不必下载每个像素，它可以只请求更小的、已经创建的概览。HTTP范围请求支持服务器上的 GeoTIFF 文件的结构使客户端能够轻松找到所需的整个文件的一部分。
当整个文件的一小部分需要处理或可视化时，瓦片就会发挥作用。这可能是概览的一部分，也可能是全分辨率。但是瓦片将一个区域的所有相关字节组织在文件的同一部分，因此范围请求可以获取它需要的内容。
如果 GeoTIFF 没有通过概览和瓦片进行“云优化”，那么对数据进行远程操作仍然有效。只是当实际只需要很小一部分数据时，他们可能会下载整个文件或大部分文件。

COG生成

可以使用GDAL生成COG。这里简单介绍下GDAL，它是一个基于C++的空间数据功能库，我这里有基于JAVA和Python的安装介绍，见Java安装GDAL依赖，python安装GDAL。
基于Java或Python开发可以将COG生成的功能比较方便地集成到现有工程里，除此之外也可以使用GDAL现成的脚本，这些脚本可以先按照Java安装GDAL依赖里下载的压缩包，解压后在release-1928-x64-dev\release-1928-x64\bin\gdal\apps下找到。

使用脚本将普通tif转成COG可以这样操作：

1. CD到脚本目录
2. 执行以下脚本，生成概览

gdaladdo -r average  你的tif路径 2 4 8 16 32

3. 执行以下脚本，生成COG

gdal_translate 输入tif的路径 输出tif的路径 -co TILED=YES -co COMPRESS=DEFLATE -co COPY_SRC_OVERVIEWS=YES

主要参数的含义：

• TILED=YES 生成瓦片
• COMPRESS=DEFLATE 使用DEFLATE方法压缩
• COPY_SRC_OVERVIEWS=YES 从原数据中拷贝概览

使用Python程序生成的代码见这篇文章

成果查看

生成的COG有多种方法检验和查看，比如使用validate_cloud_optimized_geotiff.py检验

python validate_cloud_optimized_geotiff.py test.tif

至于查看，个人感觉比较方便的是使用QGIS（3.2以上版本）。
QGIS是一个免费开源的GIS桌面端软件，和常用的ArcMap相比最大的优点肯定是免费（不需要破解 ），而且体量较轻，打开比较快，而且基本功能也比较全。
当你生成了一个COG并且启用网络容器，如Nginx，就可以用QGIS打开COG的网络地址来在线查看。



这样就可以看到了。右键图层窗口里的图层对象，查看Properties，可以查看该COG的属性信息，如范围，波段信息，空间参考等

QGIS还可以加载在线地图服务，可以看下影像位置准不准，和其他资源对比下之类的

这里做了波段拉伸，不然影像一片黑

在线渲染（基于Openlayers）

Openlayers更新了WebGLTile图层和GeoTIFFSource数据源，用geotiffjs解析tif源数据，再用WebGL渲染，实现了前端渲染GeoTIFF图像并进行波段组合，拉伸等功能。
这里的例子用的是VUE2+TypeScript
OL的引入不赘述了，下面是用到的组件

import GeoTIFFSource from "ol/source/GeoTIFF";
import Map from "ol/Map";
import WebGLTile from "ol/layer/WebGLTile";
import View from "ol/View";

获取元信息

COG文件里存储的波段信息和nodata值可以通过geotiffjs获取。
注意，这里获取是为了演示功能，如果仅为了显示不需要提前获取并使用自己指定的像素真实值范围，OL的GeoTIFFSource默认开启了normalize属性，将波段像素范围自动标准化为0-1之间的比例。

/**
 * 获取tif元信息
 * @param url
 * @returns  noData,波段信息
 */
export async function getGeoTifMetaData(url: string) {
  // 返回的是GeoTIFF对象 https://geotiffjs.github.io/geotiff.js/module-geotiff-GeoTIFF.html
  const tiff = await fromUrl(url);
  // 返回GeoTIFFImage https://geotiffjs.github.io/geotiff.js/module-geotiffimage-GeoTIFFImage.html
  const image = await tiff.getImage();
  const bands = [];
  // 获取波段数
  const samples = image.getSamplesPerPixel();
  for (let i = 0; i < samples; i++) {
    // 获取该波段信息
    const element = image.getGDALMetadata(i);
    bands.push(element);
  }
  // 获取nodata值
  const noData = image.getGDALNoData();
  return { noData, bands };
}

获取的元信息如下

编写渲染表达式

OL的WebGL图层使用表达式配置图层的渲染样式，具体的语法详见ExpressionValue

const url = "http://localhost:8082/static/cogtest/cog_origin.tif";

const metaData = await this.getGeoTifMetaData(url);
// 取出元信息里的波段最大最小值
const bands = metaData.bands.map((item) => {
  return {
    min: parseFloat(item.STATISTICS_MINIMUM),
    max: parseFloat(item.STATISTICS_MAXIMUM),
  };
});

const style = {
  color: [
    "array",
    // 这里123波段对应RGB，分别做线性拉伸，最大值*0.1为了效果比较好
    // 若初始化GeoTIFFSource时normalize为true则不需要使用像素值拉伸，而是0-1之间的比例
    ["interpolate", ["linear"], ["band", 1], bands[0].min, 0, bands[0].max * 0.1, 1],
    ["interpolate", ["linear"], ["band", 2], bands[1].min, 0, bands[1].max * 0.1, 1],
    ["interpolate", ["linear"], ["band", 3], bands[2].min, 0, bands[2].max * 0.1, 1],
    1,
  ],
};

加载图层

const source = new GeoTIFFSource({
  // 默认为true，这里为了演示使用真实值拉伸设为false
  normalize: false,
  sources: [
    {
      url,
      // 这里可以手动重设min，max，nodata等属性以覆盖tif原始信息
      nodata: 0,
    },
  ],
});
// GeoTIFF的view是异步的，里面存了COG的范围，空间参考等信息
source.getView().then((view) => {
  const layer = new WebGLTile({
    className: "WebGLTile",
    style,
    source,
    extent: view.extent,
  });

  const map = new Map({
    target: "map",
    layers: [layer],
    view: new View({
      projection: "EPSG:3857",
    }),
  });

  map.getView().fit(view.extent!, { padding: [50, 50, 50, 50] });
});

显示效果

其他功能示例可以参考OL官网的示例