卫星遥感图像瓦片处理发布——以高分二号卫星图像为例

本文目的是将高分二号卫星拍摄的图像映射到谷歌地图上，由于之前没了解过相关领域，解决起来颇为坎坷，因此将解决过程记录如下。本人非GIS相关专业，以下内容为个人理解，仅供参考。

1.数据获取

高分系列卫星，即高分专项工程。该工程全名为高分辨率对地观测系统重大专项。 该专项建立的初衷是建立一整套高时间分辨率、高空间分辨率、高光谱分辨率的自主可控卫星系列[1]，卫星状态如下：

高分二号轨道参数如下：

高分二号传感器参数如下：

如需获取指定范围的卫星数据，可在自然资源卫星遥感云服务平台（http://www.sasclouds.com/chinese/home/）进行购买，每一景影像的价格大致为660-880元[1]。

国家青藏高原科学数据中心提供了3张免费数据，下载时需要先申请，需要等待约3天左右。
国家青藏高原科学数据中心地址：https://www.tpdc.ac.cn/zh-hans/data/1b2ebe66-8389-4c9f-9756-1b29d83f851f

2.数据解析

以GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903为例，该组图包含以下文件：

其中MSS为多光谱图像，包含四个波段，分辨率为7300x6908像素，位深度为64位。
PAN为全色图像，单波段，分辨率为29200x27620像素，约为多光谱四倍，位深度为16位。

使用QGIS打开两幅图像的效果如下图所示，由于MSS为四波段图像，软件会将前三个波段视作RGB通道进行显示。

和其它遥感影响不同，其元数据信息写在xml文件中，以GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-MSS2.xml文件为例，数据内容和含义如下：

<?xml version="1.0" encoding="UTF-8" ?>
<ProductMetaData>
    <SatelliteID>GF2</SatelliteID> <!-- 卫星ID，GF2表示高分二号卫星 -->
    <SensorID>PMS2</SensorID> <!-- 传感器ID，PMS2表示传感器型号 -->
    <ReceiveTime>2020-12-01 03:17:12</ReceiveTime> <!-- 接收时间 -->
    <OrbitID>33977</OrbitID> <!-- 轨道ID，表示卫星在轨道上的位置编号 -->
    <ProduceType>STANDARD</ProduceType> <!-- 产品类型，STANDARD表示标准产品 -->
    <SceneID>8430602</SceneID> <!-- 场景ID，用于标识特定的影像场景 -->
    <ProductID>5278903</ProductID> <!-- 产品ID，用于唯一标识该产品 -->
    <ProductLevel>LEVEL1A</ProductLevel> <!-- 产品级别，LEVEL1A表示一级产品 -->
    <ProductQuality /> <!-- 产品质量信息（空白） -->
    <ProductQualityReport /> <!-- 产品质量报告（空白） -->
    <ProductFormat>GEOTIFF</ProductFormat> <!-- 产品格式，为GEOTIFF -->
    <ProduceTime>2020-12-03 13:35:50</ProduceTime> <!-- 生产时间 -->
    <Bands>1,2,3,4</Bands> <!-- 波段信息，包括1,2,3,4四个波段 -->
    <ScenePath>1014</ScenePath> <!-- 场景路径编号 -->
    <SceneRow>186</SceneRow> <!-- 场景行编号 -->
    <SatPath>1015</SatPath> <!-- 卫星路径编号 -->
    <SatRow>186</SatRow> <!-- 卫星行编号 -->
    <SceneCount>1</SceneCount> <!-- 场景数量 -->
    <SceneShift>1</SceneShift> <!-- 场景偏移 -->
    <StartTime>2020-12-01 11:17:10</StartTime> <!-- 开始时间 -->
    <EndTime>2020-12-01 11:17:13</EndTime> <!-- 结束时间 -->
    <CenterTime>2020-12-01 11:17:12</CenterTime> <!-- 中心时间 -->
    <ImageGSD>3.24</ImageGSD> <!-- 地面采样距离 -->
    <WidthInPixels>7300</WidthInPixels> <!-- 图像宽度，以像素为单位 -->
    <HeightInPixels>6908</HeightInPixels> <!-- 图像高度，以像素为单位 -->
    <WidthInMeters /> <!-- 图像宽度，以米为单位（空白） -->
    <HeightInMeters /> <!-- 图像高度，以米为单位（空白） -->
    <CloudPercent>0</CloudPercent> <!-- 云覆盖率 -->
    <QualityInfo></QualityInfo> <!-- 质量信息（空白） -->
    <PixelBits></PixelBits> <!-- 像素位深（空白） -->
    <ValidPixelBits></ValidPixelBits> <!-- 有效像素位深（空白） -->
    <RollViewingAngle>0</RollViewingAngle> <!-- 滚转观测角度 -->
    <PitchViewingAngle>0</PitchViewingAngle> <!-- 俯仰观测角度 -->
    <RollSatelliteAngle>7.9924</RollSatelliteAngle> <!-- 卫星滚转角度 -->
    <PitchSatelliteAngle>-0.00111633</PitchSatelliteAngle> <!-- 卫星俯仰角度 -->
    <YawSatelliteAngle>3.43866</YawSatelliteAngle> <!-- 卫星偏航角度 -->
    <SolarAzimuth>162.131</SolarAzimuth> <!-- 太阳方位角 -->
    <SolarZenith>46.7817</SolarZenith> <!-- 太阳天顶角 -->
    <SatelliteAzimuth>281.672</SatelliteAzimuth> <!-- 卫星方位角 -->
    <SatelliteZenith>80.4574</SatelliteZenith> <!-- 卫星天顶角 -->
    <GainMode>G2,G1,G1,G2</GainMode> <!-- 增益模式 -->
    <IntegrationTime>0.000473605</IntegrationTime> <!-- 积分时间 -->
    <IntegrationLevel>S5,S4,S5,S3</IntegrationLevel> <!-- 积分级别 -->
    <MapProjection /> <!-- 地图投影（空白） -->
    <EarthEllipsoid>WGS84</EarthEllipsoid> <!-- 地球椭球体模型 -->
    <ZoneNo /> <!-- 区号（空白） -->
    <ResamplingKernel /> <!-- 重采样核（空白） -->
    <HeightMode /> <!-- 高度模式（空白） -->
    <MtfCorrection>LAB</MtfCorrection> <!-- MTF校正方法 -->
    <RelativeCorrectionData></RelativeCorrectionData> <!-- 相对校正数据（空白） -->
    <TopLeftLatitude>23.0561</TopLeftLatitude> <!-- 左上角纬度 -->
    <TopLeftLongitude>113.933</TopLeftLongitude> <!-- 左上角经度 -->
    <TopRightLatitude>23.0098</TopRightLatitude> <!-- 右上角纬度 -->
    <TopRightLongitude>114.162</TopRightLongitude> <!-- 右上角经度 -->
    <BottomRightLatitude>22.8011</BottomRightLatitude> <!-- 右下角纬度 -->
    <BottomRightLongitude>114.112</BottomRightLongitude> <!-- 右下角经度 -->
    <BottomLeftLatitude>22.8474</BottomLeftLatitude> <!-- 左下角纬度 -->
    <BottomLeftLongitude>113.884</BottomLeftLongitude> <!-- 左下角经度 -->
    <TopLeftMapX /> <!-- 左上角地图X坐标（空白） -->
    <TopLeftMapY /> <!-- 左上角地图Y坐标（空白） -->
    <TopRightMapX /> <!-- 右上角地图X坐标（空白） -->
    <TopRightMapY /> <!-- 右上角地图Y坐标（空白） -->
    <BottomRightMapX /> <!-- 右下角地图X坐标（空白） -->
    <BottomRightMapY /> <!-- 右下角地图Y坐标（空白） -->
    <BottomLeftMapX /> <!-- 左下角地图X坐标（空白） -->
    <BottomLeftMapY /> <!-- 左下角地图Y坐标（空白） -->
</ProductMetaData>

从数据信息中，可以看出，两张图像均使用WGS84坐标系，即EPSG:4326坐标系。

3.动态瓦片显示

我的目的是用谷歌地图作为底图，在上面叠加显示更高分辨率的遥感图像。由于遥感图片分辨率过高，无法一次性完成图像加载。因此，需要使用动态瓦片显示的策略，即整幅图的显示范围是确定的，因此，只需要返回指定范围内的瓦片即可。

实现思路如下：

• 1.使用mercantile库计算出指定瓦片的地理范围。
• 2.计算地理变换矩阵的逆，以便将地理坐标转换为像素坐标。
• 3.读取图像数据：从TIFF文件中读取指定范围的像素数据。

利用flask搭建了web服务，通过folium构建图层管理和添加谷歌卫星图层。

import folium
import urllib.parse
import os
import numpy as np
import xml.etree.ElementTree as ET
import mercantile
from osgeo import gdal
from PIL import Image
from flask import Flask, render_template_string, request, jsonify, send_file
from io import BytesIO

app = Flask(__name__)



# 创建基础地图的函数
def create_map(bounds):
    center_lat = (bounds["top"] + bounds["bottom"]) / 2
    center_lon = (bounds["left"] + bounds["right"]) / 2
    print(f"Center: {center_lat}, {center_lon}")
    return folium.Map(location=[center_lat, center_lon], zoom_start=10)


# 更新地图的函数
def update_map_elements(resultpath, xml_path):
    # 检查文件是否存在
    if not os.path.exists(resultpath):
        print(f"Error: File {resultpath} does not exist")
        return None
    if not os.path.exists(xml_path):
        print(f"Error: File {xml_path} does not exist")
        return None
    # 获取TIFF文件的信息
    r = get_tiff_info(xml_path)
    # 获取边界信息
    bounds = r["bounds"]
    print(f"TIFF bounds: {bounds}")
    # 创建基础地图
    map = create_map(bounds)
    # 设置地图边界
    map.fit_bounds([[bounds["bottom"], bounds["left"]], [bounds["top"], bounds["right"]]])
    # 添加Google卫星图层
    folium.TileLayer(
        tiles="https://{s}.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
        attr="Google Satellite",
        name="Google Satellite",
        subdomains=['mt0', 'mt1', 'mt2', 'mt3'],
        overlay=False
    ).add_to(map)
    # 构建瓦片服务的URL
    tiff_tile_url_template = f'/tiff/tiles/{{z}}/{{x}}/{{y}}@1x?url=' + urllib.parse.quote(resultpath)
    # 添加TIFF图层，并包含瓦片参数
    folium.raster_layers.TileLayer(
        tiles=tiff_tile_url_template,
        attr="TIFF",
        name="Result Layer",
        overlay=True,  # 是否作为叠加层
        control=True,  # 是否在图层控制中显示
        show=True   # 是否默认显示
    ).add_to(map)

    # 添加图层控制
    folium.LayerControl().add_to(map)

    return map


# 从XML文件中获取TIFF文件信息的函数
def get_tiff_info(xml_path):
    try:
        tree = ET.parse(xml_path)
        root = tree.getroot()
        bounds = {
            "left": float(root.findtext("TopLeftLongitude")),
            "bottom": float(root.findtext("BottomRightLatitude")),
            "right": float(root.findtext("TopRightLongitude")),
            "top": float(root.findtext("TopLeftLatitude")),
        }

        width = int(root.findtext("WidthInPixels"))
        height = int(root.findtext("HeightInPixels"))

        # 计算每个像素的宽度和高度
        pixel_width = (bounds["right"] - bounds["left"]) / width
        pixel_height = (bounds["bottom"] - bounds["top"]) / height  # 通常为负数

        # 构建GeoTransform
        geotransform = [
            bounds["left"],  # top_left_x
            pixel_width,  # pixel_width
            0,  # rotation_x
            bounds["top"],  # top_left_y
            0,  # rotation_y
            pixel_height  # pixel_height
        ]
        return {"bounds": bounds, "geotransform": geotransform}
    except Exception as e:
        print(f"Error fetching TIFF info from XML: {str(e)}")
        return None


@app.route('/tiff/info')
def tiff_info():
    url = request.args.get('url')
    xml_path = url.replace('.tiff', '.xml')
    if not url:
        return jsonify({"error": "URL parameter is required"}), 400

    info = get_tiff_info(xml_path)
    if not info:
        return jsonify({"error": "Unable to fetch TIFF info"}), 500

    return jsonify(info)


def get_tile_image(url, inv_gt, left, bottom, right, top):
    dataset = gdal.Open(url)
    if not dataset:
        raise Exception(f"Unable to open TIFF file: {url}")

    # 将地理坐标转换为像素坐标
    ulx, uly = gdal.ApplyGeoTransform(inv_gt, left, top)
    lrx, lry = gdal.ApplyGeoTransform(inv_gt, right, bottom)

    ulx = max(0, int(ulx))
    uly = max(0, int(uly))
    lrx = min(dataset.RasterXSize, int(lrx))
    lry = min(dataset.RasterYSize, int(lry))

    width = lrx - ulx
    height = lry - uly

    if width <= 0 or height <= 0:
        raise Exception("Invalid window size")

    # 读取图像数据
    tile_data = dataset.ReadAsArray(ulx, uly, width, height)
    if tile_data is None or tile_data.size == 0:
        raise Exception("Failed to read tile data")

    # 处理无效值并缩放到 0-255 范围
    tile_image = tile_data.astype(np.float32)
    max_val = np.nanmax(tile_image)
    if max_val > 0:
        tile_image = (tile_image / max_val * 255).astype(np.uint8)
    else:
        tile_image = tile_image.astype(np.uint8)

    # 转换为 RGB 图像
    if len(tile_image.shape) == 2:
        tile_image = np.stack([tile_image] * 3, axis=-1)  # 灰度图转换为 RGB

    img = Image.fromarray(tile_image)
    buffer = BytesIO()
    img.save(buffer, format="PNG")
    buffer.seek(0)
    return buffer


@app.route('/tiff/tiles/<int:z>/<int:x>/<int:y>@1x')
def tiff_tiles(z, x, y):
    url = request.args.get('url')
    xml_path = url.replace('.tiff', '.xml')
    if not url:
        return jsonify({"error": "URL parameter is required"}), 400

    # 获取边界信息
    info = get_tiff_info(xml_path)
    if not info:
        return jsonify({"error": "Unable to fetch TIFF info"}), 500

    bounds = info["bounds"]
    geotransform = info["geotransform"]

    # 计算逆地理变换
    inv_gt = gdal.InvGeoTransform(geotransform)
    if not inv_gt:
        return jsonify({"error": "Failed to compute inverse GeoTransform"}), 500

    # 计算瓦片的地理边界
    tile_bounds = mercantile.bounds(mercantile.Tile(x, y, z))

    try:
        # 确保窗口不超出图像范围
        left = max(tile_bounds.west, bounds["left"])
        bottom = max(tile_bounds.south, bounds["bottom"])
        right = min(tile_bounds.east, bounds["right"])
        top = min(tile_bounds.north, bounds["top"])

        buffer = get_tile_image(url, inv_gt, left, bottom, right, top)
        return send_file(buffer, mimetype="image/png")
    except Exception as e:
        print(f"Error fetching TIFF tile: {str(e)}")
        return jsonify({"error": str(e)}), 500


@app.route('/')
def index():
    # 从URL参数中获取TIFF文件路径
    tiff_file_path = request.args.get('tiff_file', r'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2.tiff')
    xml_file_path = tiff_file_path.replace('.tiff', '.xml')

    # 更新地图元素并生成HTML
    map = update_map_elements(tiff_file_path, xml_file_path)
    if map is None:
        return "Error generating map. Please check the console for details."
    map_html = map._repr_html_()

    # 使用模板展示地图
    return render_template_string('''
        <!DOCTYPE html>
        <html>
        <head>
            <title>Map</title>
            <meta charset="utf-8">
            <meta name="viewport" content="width=device-width, initial-scale=1.0">
            {{ map_html|safe }}
        </head>
        <body>
        </body>
        </html>
    ''', map_html=map_html)


if __name__ == '__main__':
    app.run(host='0.0.0.0', port=5000)

效果显示如下：

虽然显示效果还可以，不过速度过于缓慢。主要原因是每次请求区域都需要进行计算，并且没有预先设定瓦片的大小，导致显示的瓦片没有经过下采样，而是用原分辨率直接显示，内存消耗也很大。

4.COG动态显示

为了解决瓦片显示过慢的问题，一个解决思路就是是否可以将瓦片提前处理好，这样加载的时候就不用重新计算了，一种思路方式是将数据预处理为COG形式。

4.1 COG格式简介

COG（Cloud Optimized GeoTIFF）是一种优化的GeoTIFF文件格式，专门用于在云存储和Web应用中高效访问和处理地理空间数据。

COG格式具有以下特点：

• 分块和压缩
  COG文件内部数据被划分为多个独立的块（tiles），每个块都可以独立压缩。这使得可以只读取和处理感兴趣的特定区域，而无需加载整个文件。

• 内置索引（Overviews）
  COG文件通常包含多个分辨率级别的预览图（overviews），允许快速访问不同缩放级别的数据，提高了读取效率和显示速度。

• HTTP Range 请求支持：
  COG文件支持HTTP Range请求，这意味着可以通过HTTP协议请求文件的特定部分。这在云存储和Web服务中尤为重要，可以显著减少数据传输量和访问时间。

• 元数据标准化：
  COG文件包含标准化的元数据，方便数据共享和互操作性。

优势：

• 高效读取：因为支持分块和多分辨率预览，COG文件可以快速读取特定区域的数据，适合大规模地理空间数据的在线处理和显示。
• 云存储优化：支持HTTP Range请求，适合在云存储中高效存储和访问，减少了数据传输和存储成本。
• 兼容性：COG文件是标准的GeoTIFF文件，因此与现有的GeoTIFF工具和库兼容，同时优化了在线使用的性能。

4.2 COG格式转换

下面通过Gdal将原始数据转成COG格式。

根据Gdal官方文档[2]，gdal 3.1版本更新了COG栅格driver，对于一般包含元数据信息的tif文件，可通过gdalwarp和gdal_translate进行转换。

gdalwarp src1.tif src2.tif out.tif -of COG
gdal_translate world.tif world_webmerc_cog.tif -of COG -co TILING_SCHEME=GoogleMapsCompatible -co COMPRESS=JPEG

不过由于高分二号卫星的元信息数据主要在xml中，因此需要将转换所需信息从xml中提取，转换代码如下：

import time
from osgeo import gdal, osr
import xml.etree.ElementTree as ET


def translateToCOG(in_ds, out_path):
    if in_ds is None:
        print("Error: Input dataset is None")
        return

    im_bands = in_ds.RasterCount
    for i in range(im_bands):
        band = in_ds.GetRasterBand(i + 1)
        if band is None:
            print(f"Error: Band {i + 1} is None")
            continue

        nodataVal = band.GetNoDataValue()
        maxBandValue = band.GetMaximum()

        if maxBandValue is None:
            band.ComputeStatistics(0)
        if nodataVal is None:
            band.SetNoDataValue(0.0)

    try:
        in_ds.BuildOverviews("NEAREST", [2, 4, 8, 16])
    except Exception as e:
        print(f"Error building overviews: {e}")
        return

    driver = gdal.GetDriverByName('GTiff')
    if driver is None:
        print("Error: Could not get GTiff driver")
        return

    try:
        driver.CreateCopy(out_path, in_ds, options=[
            "COPY_SRC_OVERVIEWS=YES",
            "TILED=YES",
            "COMPRESS=DEFLATE",
            "INTERLEAVE=BAND"
        ])
    except Exception as e:
        print(f"Error creating COG: {e}")


def warpDataset(in_ds, proj, resampling=1):
    RESAMPLING_MODEL = ['', gdal.GRA_NearestNeighbour,
                        gdal.GRA_Bilinear, gdal.GRA_Cubic]

    resampleAlg = RESAMPLING_MODEL[resampling]

    try:
        return gdal.AutoCreateWarpedVRT(in_ds, None, proj, resampleAlg)
    except Exception as e:
        print(f"Error warping dataset: {e}")
        return None


def getTifDataset(fileDir, srid=None):
    dataset = gdal.Open(fileDir, gdal.GA_ReadOnly)
    if dataset is None:
        print(fileDir + "文件无法打开")
        return None

    fileSrs = osr.SpatialReference()
    fileSrs.ImportFromWkt(dataset.GetProjection())

    if srid is None:
        return dataset
    else:
        outSrs = osr.SpatialReference()
        outSrs.ImportFromEPSG(srid)

        if fileSrs.IsSame(outSrs):
            return dataset
        else:
            return warpDataset(dataset, outSrs.ExportToWkt())


def set_geo_transform_from_xml(dataset, xml_path):
    if dataset is None:
        print("Error: Dataset is None, cannot set geo transform")
        return

    tree = ET.parse(xml_path)
    root = tree.getroot()

    top_left_lon = float(root.findtext("TopLeftLongitude"))
    top_left_lat = float(root.findtext("TopLeftLatitude"))
    top_right_lon = float(root.findtext("TopRightLongitude"))
    bottom_right_lat = float(root.findtext("BottomRightLatitude"))

    width_in_pixels = int(root.findtext("WidthInPixels"))
    height_in_pixels = int(root.findtext("HeightInPixels"))

    lon_res = (top_right_lon - top_left_lon) / width_in_pixels
    lat_res = (top_left_lat - bottom_right_lat) / height_in_pixels

    geo_transform = [
        top_left_lon, lon_res, 0,
        top_left_lat, 0, -lat_res
    ]

    try:
        dataset.SetGeoTransform(geo_transform)
    except Exception as e:
        print(f"Error setting geo transform: {e}")


def set_projection_from_xml(dataset, xml_path):
    if dataset is None:
        print("Error: Dataset is None, cannot set projection")
        return

    tree = ET.parse(xml_path)

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    try:
        dataset.SetProjection(srs.ExportToWkt())
    except Exception as e:
        print(f"Error setting projection: {e}")


if __name__ == '__main__':
    start = time.process_time()
    inPath = "GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2.tiff"
    xmlPath = "GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2.xml"
    outputPath = "GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_cog.tiff"

    originDataset = getTifDataset(inPath)
    if originDataset is None:
        print(f"Error: Unable to open dataset {inPath}")
    else:
        set_geo_transform_from_xml(originDataset, xmlPath)
        set_projection_from_xml(originDataset, xmlPath)
        translateToCOG(originDataset, outputPath)

    originDataset = None
    end = time.process_time()
    print('Running time: %s Seconds' % (end - start))

转换完之后，可以用下面的脚本对cog格式进行验证。
cog格式的验证思路来自gdal源码：https://github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/osgeo_utils/samples/validate_cloud_optimized_geotiff.py

主要检查文件中IFD偏移量及尺寸是否符合COG格式规范，具体内容可参考注释信息：

from osgeo import gdal
import struct

def get_block_offset(band):
    blockxsize, blockysize = band.GetBlockSize()
    for y in range(int((band.YSize + blockysize - 1) / blockysize)):
        for x in range(int((band.XSize + blockxsize - 1) / blockxsize)):
            block_offset = band.GetMetadataItem(
                "BLOCK_OFFSET_%d_%d" % (x, y), "TIFF"
            )
            if block_offset:
                return int(block_offset)
    return 0

def validate_cog(file_path):
    ds = gdal.Open(file_path)
    if ds is None:
        return False, "Cannot open file"
    # 检查 GDAL version Cog功能在GDAL 2.2以上引入
    if int(gdal.VersionInfo("VERSION_NUM")) < 2020000:
        return False, "GDAL 2.2 or above required"
    # 检查文件是否为GeoTIFF
    if ds.GetDriver().ShortName != "GTiff":
        return False, "The file is not a GeoTIFF"
    errors = []
    warnings = []
    filename = ds.GetDescription()
    main_band = ds.GetRasterBand(1)
    ovr_count = main_band.GetOverviewCount()
    filelist = ds.GetFileList()
    if filelist is not None and filename + ".ovr" in filelist:
        errors.append("Overviews found in external .ovr file. They should be internal")
    # 检查图像块是否大于512x512
    if main_band.XSize > 512 or main_band.YSize > 512:
        block_size = main_band.GetBlockSize()
        if block_size[0] == main_band.XSize and block_size[0] > 1024:
            errors.append("The file is greater than 512xH or Wx512, but is not tiled")
        if ovr_count == 0:
            warnings.append(
                "The file is greater than 512xH or Wx512, it is recommended "
                "to include internal overviews"
            )
    # 检查主IFD偏移量
    # 主IFD偏移量用来告诉文件解析器从文件的哪个位置开始读取图像的元数据
    ifd_offset = int(main_band.GetMetadataItem("IFD_OFFSET", "TIFF"))
    ifd_offsets = [ifd_offset]

    # 检查文件的结构化元数据
    block_order_row_major = False  # 是否按行主序排列数据块
    block_leader_size_as_uint4 = False  # 数据块前是否有一个4字节的大小前缀
    block_trailer_last_4_bytes_repeated = False  #  数据块后4字节是否重复
    mask_interleaved_with_imagery = False  # 是否图像和掩膜数据交错存储

    # 检查主IFD偏移量是否在预期的8或16字节位置(经典TIFF或BigTIFF)，如果不在，则需进一步检查文件结构
    if ifd_offset not in (8, 16):
        f = gdal.VSIFOpenL(filename, "rb")
        if not f:
            return False, "Cannot open file"
        # 读取文件的前4字节来判断文件是TIFF还是BigTIFF
        signature = struct.unpack("B" * 4, gdal.VSIFReadL(4, 1, f))
        bigtiff = signature in ((0x49, 0x49, 0x2B, 0x00), (0x4D, 0x4D, 0x00, 0x2B))
        expected_ifd_pos = 16 if bigtiff else 8
        gdal.VSIFSeekL(f, expected_ifd_pos, 0)
        pattern = "GDAL_STRUCTURAL_METADATA_SIZE=%06d bytes\n" % 0
        got = gdal.VSIFReadL(len(pattern), 1, f).decode("LATIN1")
        if len(got) == len(pattern) and got.startswith("GDAL_STRUCTURAL_METADATA_SIZE="):
            size = int(got[len("GDAL_STRUCTURAL_METADATA_SIZE="):][0:6])
            extra_md = gdal.VSIFReadL(size, 1, f).decode("LATIN1")
            block_order_row_major = "BLOCK_ORDER=ROW_MAJOR" in extra_md
            block_leader_size_as_uint4 = "BLOCK_LEADER=SIZE_AS_UINT4" in extra_md
            block_trailer_last_4_bytes_repeated = "BLOCK_TRAILER=LAST_4_BYTES_REPEATED" in extra_md
            mask_interleaved_with_imagery = "MASK_INTERLEAVED_WITH_IMAGERY=YES" in extra_md
            if "KNOWN_INCOMPATIBLE_EDITION=YES" in extra_md:
                errors.append("KNOWN_INCOMPATIBLE_EDITION=YES is declared in the file")
            expected_ifd_pos += len(pattern) + size
            expected_ifd_pos += (expected_ifd_pos % 2)
        gdal.VSIFCloseL(f)

        if expected_ifd_pos != ifd_offsets[0]:
            errors.append(
                "The offset of the main IFD should be %d. It is %d instead" % (expected_ifd_pos, ifd_offsets[0])
            )

    # 检查概览
    for i in range(ovr_count):
        ovr_band = ds.GetRasterBand(1).GetOverview(i)
        # 检查第一个概览的尺寸。如果第一个概览的尺寸大于主波段的尺寸则添加错误信息。因为概览的分辨率应该总是低于或等于主波段的分辨率
        if i == 0:
            if ovr_band.XSize > main_band.XSize or ovr_band.YSize > main_band.YSize:
                errors.append("First overview has larger dimension than main band")
        # 检查后续概览的尺寸，检查它们的尺寸是否小于前一个概览的尺寸，概览的尺寸应逐渐减小。
        else:
            prev_ovr_band = ds.GetRasterBand(1).GetOverview(i - 1)
            if ovr_band.XSize > prev_ovr_band.XSize or ovr_band.YSize > prev_ovr_band.YSize:
                errors.append(
                    "Overview of index %d has larger dimension than overview of index %d" % (i, i - 1)
                )
        # 检查概览的块大小，如果块大小等于概览的宽度且大于1024，表示这个概览不是按块存储的，而是按行存储的，这不符合COG标准
        block_size = ovr_band.GetBlockSize()
        if block_size[0] == ovr_band.XSize and block_size[0] > 1024:
            errors.append("Overview of index %d is not tiled" % i)

        ifd_offset = int(ovr_band.GetMetadataItem("IFD_OFFSET", "TIFF"))
        ifd_offsets.append(ifd_offset)

        # 检查当前概览的IFD偏移量是否大于前一个概览或主图像的IFD偏移量，因为IFD偏移量应该按顺序递增
        if ifd_offsets[-1] < ifd_offsets[-2]:
            if i == 0:
                errors.append(
                    "The offset of the IFD for overview of index %d is %d, whereas it should be greater than the one of the main image, which is at byte %d" % (i, ifd_offsets[-1], ifd_offsets[-2])
                )
            else:
                errors.append(
                    "The offset of the IFD for overview of index %d is %d, whereas it should be greater than the one of index %d, which is at byte %d" % (i, ifd_offsets[-1], i - 1, ifd_offsets[-2])
                )
    # 检查数据块的偏移量
    block_offset = get_block_offset(main_band)
    data_offsets = [block_offset]
    for i in range(ovr_count):
        ovr_band = ds.GetRasterBand(1).GetOverview(i)
        block_offset = get_block_offset(ovr_band)
        data_offsets.append(block_offset)
    # 检查最小概览的第一个数据块偏移量是否在其IFD之后
    if data_offsets[-1] != 0 and data_offsets[-1] < ifd_offsets[-1]:
        if ovr_count > 0:
            errors.append("The offset of the first block of the smallest overview should be after its IFD")
        else:
            errors.append("The offset of the first block of the image should be after its IFD")
    # 检查概览的数据块偏移量顺序
    for i in range(len(data_offsets) - 2, 0, -1):
        if data_offsets[i] != 0 and data_offsets[i] < data_offsets[i + 1]:
            errors.append(
                "The offset of the first block of overview of index %d should be after the one of the overview of index %d" % (i - 1, i)
            )
    # 检查主图像的数据块偏移量顺序
    if len(data_offsets) >= 2 and data_offsets[0] != 0 and data_offsets[0] < data_offsets[1]:
        errors.append(
            "The offset of the first block of the main resolution image should be after the one of the overview of index %d" % (ovr_count - 1)
        )
    if errors:
        return False, "Errors found: " + "; ".join(errors)
    return True, "File is a valid Cloud Optimized GeoTIFF"

# 输入文件路径
file_path = "GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_cog.tiff"

is_cog, message = validate_cog(file_path)
if is_cog:
    print(f"{file_path} is a Cloud Optimized GeoTIFF.")
else:
    print(f"{file_path} is NOT a Cloud Optimized GeoTIFF. Reason: {message}")

下面这段代码可以辅助理解COG格式具体是如何运作的，选取全图范围(也可以自定义一块范围)，COG将该范围预先处理成不同分辨率的图像，以满足不同缩放比例下的渲染。

from osgeo import gdal
import matplotlib.pyplot as plt

# 输入文件路径
file_path = "GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_cog.tiff"

# 打开文件
dataset = gdal.Open(file_path, gdal.GA_ReadOnly)
if not dataset:
    raise Exception(f"Unable to open file: {file_path}")

# 获取主波段
main_band = dataset.GetRasterBand(1)

# 获取概览图像数量
overview_count = main_band.GetOverviewCount()

# 定义感兴趣区域 (例如一个瓦片)
tile_x_offset = 1000
tile_y_offset = 1000
tile_width = 512
tile_height = 512

if overview_count == 0:
    print("No overviews found in the file.")
else:
    fig, axs = plt.subplots(1, overview_count, figsize=(15, 5))
    fig.suptitle('Overview Images at Different Resolutions')

    # 读取并显示主波段中的感兴趣区域
    main_data = main_band.ReadAsArray(tile_x_offset, tile_y_offset, tile_width, tile_height)
    axs[0].imshow(main_data, cmap='gray')
    axs[0].set_title(f"Main\n{main_band.YSize} x {main_band.XSize}")
    axs[0].axis('off')

    for i in range(overview_count):
        overview_band = main_band.GetOverview(i)
        # overview_data = overview_band.ReadAsArray()

        # 计算对应概览中的区域
        ov_x_offset = tile_x_offset // (2 ** (i + 1))
        ov_y_offset = tile_y_offset // (2 ** (i + 1))
        ov_width = tile_width // (2 ** (i + 1))
        ov_height = tile_height // (2 ** (i + 1))

        overview_data = overview_band.ReadAsArray(ov_x_offset, ov_y_offset, ov_width, ov_height)

        # 获取分辨率
        x_res = overview_band.XSize
        y_res = overview_band.YSize

        # 显示概览图像
        ax = axs[i]
        ax.imshow(overview_data, cmap='gray')
        ax.set_title(f"Level {i + 1}\n{y_res} x {x_res}")
        ax.axis('off')

    plt.tight_layout()
    plt.show()

print("Done displaying overviews.")

4.3 坐标系格式转换(EPSG:4326->EPSG:3857)

在地理信息系统（GIS）中，常用的两个坐标系是WGS84（EPSG:4326）和Web Mercator（EPSG:3857）。WGS84是全球通用的地理坐标系，使用经纬度（度）来表示位置；而Web Mercator是投影坐标系，通常用于在线地图服务，使用米为单位表示位置。

因此，在发布之前，需要先将EPSG:4326转换成EPSG:3857坐标系。

转换公式如下：

# 地球赤道半径（米）
EARTH_EQUATORIAL_RADIUS = 6378137.0
def from4326_to3857(lat, lon):
    xtile = math.radians(lon) * EARTH_EQUATORIAL_RADIUS
    ytile = math.log(math.tan(math.radians(45 + lat / 2.0))) * EARTH_EQUATORIAL_RADIUS
    return xtile, ytile

• xtile：将经度转换为弧度后乘以地球赤道半径，得到Web Mercator坐标系中的x坐标。
• ytile：将纬度转换为弧度后通过数学公式得到y坐标。

下面采用gdal的Warp进行转换

import math
from osgeo import gdal


def reproject_tiff(input_file, output_file, target_srs='EPSG:3857'):
    # 打开输入TIFF文件
    dataset = gdal.Open(input_file)

    if not dataset:
        print(f"无法打开文件: {input_file}")
        return

    # 获取输入文件的投影信息
    input_srs = dataset.GetProjection()

    # 设置转换选项
    warp_options = gdal.WarpOptions(dstSRS=target_srs, format='GTiff')

    # 执行转换
    gdal.Warp(output_file, dataset, options=warp_options)

    # 再次打开输出文件并更新其元数据
    output_dataset = gdal.Open(output_file, gdal.GA_Update)
    output_dataset.SetProjection(target_srs)

    print(f"文件已成功转换为{target_srs}坐标系，并保存为{output_file}")


if __name__ == '__main__':
    input_file = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_cog.tiff'
    output_file = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_cog_3857.tiff'

    # 转换为3857坐标系
    reproject_tiff(input_file, output_file)

4.4 发布显示

通过Flask框架发布显示，代码如下：

import folium
import urllib.parse
import os
import numpy as np
import mercantile
from osgeo import gdal, osr
from PIL import Image
from flask import Flask, render_template_string, request, jsonify, send_file
from io import BytesIO

app = Flask(__name__)

# 创建基础地图的函数
def create_map(bounds):
    center_lat = (bounds["top"] + bounds["bottom"]) / 2
    center_lon = (bounds["left"] + bounds["right"]) / 2
    return folium.Map(location=[center_lat, center_lon], zoom_start=10)

# 更新地图的函数
def update_map_elements(resultpath):
    if not os.path.exists(resultpath):
        print(f"Error: File {resultpath} does not exist")
        return None

    r = get_tiff_info(resultpath)
    if not r:
        return None

    bounds = r["bounds"]
    print(f"TIFF bounds: {bounds}")

    map = create_map(bounds)
    map.fit_bounds([[bounds["bottom"], bounds["left"]], [bounds["top"], bounds["right"]]])

    folium.TileLayer(
        tiles="https://{s}.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
        attr="Google Satellite",
        name="Google Satellite",
        subdomains=['mt0', 'mt1', 'mt2', 'mt3'],
        overlay=False
    ).add_to(map)

    tiff_tile_url_template = f'/tiff/tiles/{{z}}/{{x}}/{{y}}@1x?url=' + urllib.parse.quote(resultpath)

    folium.raster_layers.TileLayer(
        tiles=tiff_tile_url_template,
        attr="TIFF",
        name="Result Layer",
        overlay=True,
        control=True,
        show=True
    ).add_to(map)

    folium.LayerControl().add_to(map)
    return map

# 获取TIFF文件信息的函数
def get_tiff_info(tiff_path):
    dataset = gdal.Open(tiff_path)
    if not dataset:
        raise Exception(f"Unable to open TIFF file: {tiff_path}")

    gt = dataset.GetGeoTransform()
    proj = dataset.GetProjection()
    width = dataset.RasterXSize
    height = dataset.RasterYSize

    min_x = gt[0]
    max_y = gt[3]
    max_x = min_x + width * gt[1]
    min_y = max_y + height * gt[5]

    bounds = {
        "left": min_x,
        "bottom": min_y,
        "right": max_x,
        "top": max_y
    }

    return {"bounds": bounds, "geotransform": gt, "projection": proj}

def reproject_to_epsg3857(src_dataset):
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(3857)  # EPSG:3857 for Web Mercator

    reprojected_ds = gdal.AutoCreateWarpedVRT(src_dataset, None, dst_srs.ExportToWkt(), gdal.GRA_NearestNeighbour)
    return reprojected_ds

def get_tile_image(url, x, y, z):
    # 使用GDAL的Translate功能直接从COG文件中读取瓦片数据
    cog_file = gdal.Open(url)
    if not cog_file:
        raise Exception(f"Unable to open TIFF file: {url}")

    # 重新投影为EPSG:3857
    cog_file = reproject_to_epsg3857(cog_file)

    # 计算瓦片的地理边界
    tile_bounds = mercantile.bounds(mercantile.Tile(x, y, z))

    # 确保边界在图像范围内
    gt = cog_file.GetGeoTransform()
    min_x = gt[0]
    max_y = gt[3]
    max_x = min_x + cog_file.RasterXSize * gt[1]
    min_y = max_y + cog_file.RasterYSize * gt[5]

    left = max(tile_bounds.west, min_x)
    bottom = max(tile_bounds.south, min_y)
    right = min(tile_bounds.east, max_x)
    top = min(tile_bounds.north, max_y)

    # 确保窗口大小为正值
    if right <= left or bottom >= top:
        print("Window size is invalid, returning empty tile.")
        return BytesIO()

    translate_options = gdal.TranslateOptions(
        format='MEM',
        projWin=[left, top, right, bottom],
        width=256,
        height=256
    )

    # 将瓦片数据转换为内存文件
    mem_ds = gdal.Translate('', cog_file, options=translate_options)
    if mem_ds is None:
        raise Exception("Failed to translate COG to memory dataset")

    data = mem_ds.ReadAsArray()
    if data is None:
        raise Exception("Failed to read array from memory dataset")

    if data.ndim == 2:
        data = np.stack([data] * 3, axis=-1)
    elif data.ndim == 3:
        if data.shape[0] == 1:
            data = np.stack([data[0]] * 3, axis=-1)
        elif data.shape[0] == 3 or data.shape[0] == 4:
            data = data.transpose(1, 2, 0)
        else:
            raise Exception("Unexpected data shape")
    else:
        raise Exception("Unexpected data shape")

    max_val = np.nanmax(data)
    if max_val > 0:
        data = (data / max_val * 255).astype(np.uint8)
    else:
        data = data.astype(np.uint8)

    buffer = BytesIO()
    img = Image.fromarray(data)
    img.save(buffer, format="PNG")
    buffer.seek(0)

    return buffer

@app.route('/tiff/tiles/<int:z>/<int:x>/<int:y>@1x')
def tiff_tiles(z, x, y):
    url = request.args.get('url')
    if not url:
        return jsonify({"error": "URL parameter is required"}), 400

    try:
        buffer = get_tile_image(url, x, y, z)
        return send_file(buffer, mimetype="image/png")
    except Exception as e:
        print(f"Error fetching TIFF tile: {str(e)}")
        return jsonify({"error": str(e)}), 500

@app.route('/')
def index():
    tiff_file_path = request.args.get('tiff_file', r'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_cog_3857.tiff')
    map = update_map_elements(tiff_file_path)
    if map is None:
        return "Error generating map. Please check the console for details."
    map_html = map._repr_html_()
    return render_template_string('''
        <!DOCTYPE html>
        <html>
        <head>
            <title>Map</title>
            <meta charset="utf-8">
            <meta name="viewport" content="width=device-width, initial-scale=1.0">
            {{ map_html|safe }}
        </head>
        <body>
        </body>
        </html>
    ''', map_html=map_html)

if __name__ == '__main__':
    app.run(host='0.0.0.0', port=5000)

加载速度很快，不过仍存在图像尺寸异常的问题。

5.XYZ动态显示

本想顺着cog的思路往下做，可惜一直未能解决问题。归根原因是我对cog的格式还不够了解，因此想换一种更直观的方式，比如，是否可以将瓦片以文件的形式存储出来，这样看起来就比较清晰。于是，查阅到了使用XYZ进行切片的方式：

5.1 xyz简介

XYZ瓦片机制（XYZ tiling scheme）是一种用于管理和组织大规模地理空间数据的方法。它主要用于将地球表面的三维地理空间数据，如卫星影像、地形数据和建筑模型等，分割成小块状的瓦片（tiles），以便于存储、传输和处理。

这种机制通常涉及将地球表面划分为不同级别的瓦片，每个级别的瓦片分辨率逐渐增加。例如，XYZ瓦片可以根据所需分辨率动态加载和显示，从而支持多层级的地理信息展示和交互。XYZ代表了在地理空间数据中使用的坐标系统：X表示水平位置，Y表示垂直位置，Z表示缩放级别。

经纬度到瓦片坐标的转换公式：

tile_x = floor((longitude + 180) / 360 * (2^zoom))
tile_y = floor((1 - log(tan(latitude * pi / 180) + 1 / cos(latitude * pi / 180)) / pi) / 2 * (2^zoom))

• longitude：经度，取值范围为 -180 到 180 度。
• latitude：纬度，取值范围为 -90 到 90 度。
• zoom：缩放级别，通常是一个整数值，表示地图的放大倍数。

瓦片坐标到经纬度的逆转换公式：

n = 2^zoom
lon_deg = tile_x / n * 360.0 - 180.0
lat_rad = atan(sinh(pi * (1 - 2 * tile_y / n)))
lat_deg = lat_rad * 180.0 / pi

• tile_x：瓦片的X坐标。
• tile_y：瓦片的Y坐标。
• zoom：缩放级别，与转换时相同的缩放级别。

5.2 切分瓦片

使用gdal2tiles可以很方便的切分瓦片，比如，我想要将全色图像切分成1-8的缩放层级的瓦片，一行代码就可以完成：

import gdal2tiles
import time

if __name__ == '__main__':
    s_t = time.time()
    gdal2tiles.generate_tiles('GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2.tiff', 'demo/', np_processes=16, zoom='1-17')
    e_t = time.time()
    print(f'Time: {e_t - s_t}s')

完成之后，它会生成以下这些文件：

打开leaflet.html，就可以查看效果：

用同样的方式试了下多光谱图像，颜色很奇怪，可能是无法直接处理4波段的信息：

5.3 QGIS切分

为了解决多光谱的颜色问题，使用QGIS自带的工具箱进行XYZ瓦片切分，可以在这里找到：

由于高分二号图像的坐标系写在xml文件中，在处理时，QGIS无法获取到，因此先将xml信息假如到tiff图像中：

from osgeo import gdal, osr
import xml.etree.ElementTree as ET

def add_geoinfo_from_xml_to_tiff(xml_file, input_tiff, output_tiff):
    # 解析XML文件
    tree = ET.parse(xml_file)
    root = tree.getroot()

    # 提取经纬度信息
    latitudes = []
    longitudes = []
    for tag in ['TopLeftLatitude', 'TopRightLatitude', 'BottomRightLatitude', 'BottomLeftLatitude']:
        latitudes.append(float(root.find(tag).text))
    for tag in ['TopLeftLongitude', 'TopRightLongitude', 'BottomRightLongitude', 'BottomLeftLongitude']:
        longitudes.append(float(root.find(tag).text))

    # 打开TIFF文件
    dataset = gdal.Open(input_tiff, gdal.GA_Update)
    if not dataset:
        print("无法打开TIFF文件.")
        return

    # 设置坐标系为WGS84
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    dataset.SetProjection(srs.ExportToWkt())

    # 设置地理变换参数
    geotransform = [
        longitudes[0],
        (longitudes[1] - longitudes[0]) / dataset.RasterXSize,
        0,
        latitudes[0],
        0,
        (latitudes[3] - latitudes[0]) / dataset.RasterYSize
    ]
    dataset.SetGeoTransform(geotransform)

    # 创建一个拷贝，保存到输出文件
    driver = gdal.GetDriverByName('GTiff')
    driver.CreateCopy(output_tiff, dataset)

    # 关闭TIFF文件
    dataset = None

if __name__ == '__main__':
    xml_file = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2.xml'
    input_tiff = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2.tiff'
    output_tiff = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_withxml.tiff'

    add_geoinfo_from_xml_to_tiff(xml_file, input_tiff, output_tiff)

现在QGIS就可以直接读取图层坐标信息：

Minimum zoom表示最小缩放比例，设置为1。
Maximum zoom表示最大缩放比例，为了追求精细可以设置为17。
默认情况y轴是从下往上增大，为了适配谷歌地图，需要勾选Use inverted tile Y axis，即Y轴翻转。
Ouput directory设置保存文件夹路径，即可点击运行，让它进行切片。

为了进行直观可视化，对之前的html进行改进，代码如下：

<!DOCTYPE html>
<html lang="en">
<head>
    <meta charset="utf-8">
    <meta name='viewport' content='width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=no' />
    <title>Layer Switcher</title>
    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.5/leaflet.css" />
    <script src="http://cdn.leafletjs.com/leaflet-0.7.5/leaflet.js"></script>
    <style>
        body { margin:0; padding:0; }
        body, table, tr, td, th, div, h1, h2, input { font-family: "Calibri", "Trebuchet MS", "Ubuntu", Serif; font-size: 11pt; }
        #map { position:absolute; top:0; bottom:0; width:100%; }
        .ctl {
            padding: 2px 10px 2px 10px;
            background: white;
            background: rgba(255,255,255,0.9);
            box-shadow: 0 0 15px rgba(0,0,0,0.2);
            border-radius: 5px;
            text-align: right;
        }
        .title {
            font-size: 18pt;
            font-weight: bold;
        }
        .src {
            font-size: 10pt;
        }
    </style>
</head>
<body>

<div id="map"></div>

<script>
// Base layers
var osm = L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {attribution: '&copy; <a href="http://osm.org/copyright">OpenStreetMap</a> contributors'});
var googleSat = L.tileLayer('https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}', {
    maxZoom: 20,
    attribution: 'Google Satellite'
});

// Overlay layers (TMS)
var lyr1 = L.tileLayer('./layer1/{z}/{x}/{y}.png', {tms: true, opacity: 0.7, attribution: ""});
var lyr2 = L.tileLayer('./layer2/{z}/{x}/{y}.png', {tms: true, opacity: 0.7, attribution: ""});

// Map
var map = L.map('map', {
    center: [35.99680347453017, 116.82861328125],
    zoom: 20,
    minZoom: 1,
    maxZoom: 20,
    layers: [osm, lyr1]
});

var basemaps = {"OpenStreetMap": osm, "Google Satellite": googleSat};
var overlaymaps = {"Layer 1": lyr1, "Layer 2": lyr2};

// Title
var title = L.control();
title.onAdd = function(map) {
    this._div = L.DomUtil.create('div', 'ctl title');
    this.update();
    return this._div;
};
title.update = function(props) {
    this._div.innerHTML = "Layer Switcher Example";
};
title.addTo(map);

// Note
var src = 'Generated by GDAL2Tiles, Copyright &copy; 2008 Klokan Petr Pridal, GDAL & OSGeo GSoC';
var title = L.control({position: 'bottomleft'});
title.onAdd = function(map) {
    this._div = L.DomUtil.create('div', 'ctl src');
    this.update();
    return this._div;
};
title.update = function(props) {
    this._div.innerHTML = src;
};
title.addTo(map);

// Add base layers
L.control.layers(basemaps, overlaymaps, {collapsed: false}).addTo(map);

// Fit to overlay bounds
map.fitBounds([[35.875698032496665, 116.9549560546875], [36.117908916563664, 116.7022705078125]]);
</script>

</body>
</html>

html文件放置在templates文件夹下。

设置flask启动脚本：

from flask import Flask, render_template

app = Flask(__name__, static_folder='static', static_url_path='/')

@app.route('/')
def index():
    return render_template('index.html')


if __name__ == '__main__':
    app.run(debug=True)

生成的切片资源放置在stastic静态资源文件夹下，结构如下图所示：

为了将多光谱图像和全色图像对齐，将全色图像进行下采样：

from osgeo import gdal

def downsample_tiff(input_tiff, output_tiff, width, height):
    # 打开输入文件
    dataset = gdal.Open(input_tiff)
    if not dataset:
        print(f"无法打开文件: {input_tiff}")
        return

    # 设置转换选项
    options = gdal.TranslateOptions(
        width=width,
        height=height,
        resampleAlg='bilinear'  # 使用双线性重采样
    )

    # 执行转换
    gdal.Translate(output_tiff, dataset, options=options)
    print(f"生成下采样文件: {output_tiff}")

if __name__ == '__main__':
    input_tiff = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2.tiff'
    output_tiff = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_downsampled.tiff'
    width = 7300
    height = 6908

    downsample_tiff(input_tiff, output_tiff, width, height)

但由于坐标参数的差异，这两幅图像的并不能严格对齐，导入QGIS中，可明显发现无法完全重合。

因此，想的一个办法就是将多光谱的地理信息完全拷贝给全色图像，代码如下：

from osgeo import gdal

def copy_geoinfo(source_tiff, target_tiff, output_tiff):
    # 打开源TIFF文件和目标TIFF文件
    src_ds = gdal.Open(source_tiff, gdal.GA_ReadOnly)
    tgt_ds = gdal.Open(target_tiff, gdal.GA_Update)

    if not src_ds or not tgt_ds:
        print("无法打开TIFF文件")
        return

    # 复制坐标系统信息
    tgt_ds.SetProjection(src_ds.GetProjection())
    tgt_ds.SetGeoTransform(src_ds.GetGeoTransform())

    # 复制元数据
    tgt_ds.SetMetadata(src_ds.GetMetadata())

    # 关闭TIFF文件
    tgt_ds = None
    src_ds = None

    # 保存输出文件
    gdal.Translate(output_tiff, target_tiff, format='GTiff')

if __name__ == '__main__':
    source_tiff = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-MSS2.tiff'  # 包含正确地理信息的TIFF文件
    target_tiff = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_downsampled.tiff'  # 需要更新地理信息的TIFF文件
    output_tiff = 'GF2_PMS2_E114.0_N22.9_20201201_L1A0005278903-PAN2_downsampled_geoinfo.tiff'  # 输出文件
    copy_geoinfo(source_tiff, target_tiff, output_tiff)

最终效果：

仔细看还是会有一些位置偏差，可能是由于传感器的误差，或者是未对图像做精确校正。

总结

1. 使用QGIS主要为了解决多光谱图像的颜色问题，由于QGIS是开源软件，因此如果需要进一步深究QGIS是如何进行处理的，可参考QGIS仓库[4]的
   QGIS/src/analysis/processing/qgsalgorithmxyztiles.cpp 相关代码：
   链接：https://github.com/qgis/QGIS/blob/711c75d87377d0c94786559a67000ae26b9407ca/src/analysis/processing/qgsalgorithmxyztiles.cpp#L256

2. 本文的一些操作纯属个人理解，如需更专业的GIS操作，在处理之前，可以对图像进行辐射定标、大气校正、正射校正等操作，具体可参考[5]。另外，全色图像舍弃了颜色信息换来了更高的分辨率，而多光谱图像则补充了颜色信息，如何进行全色图像和多光谱图像进行融合，可能是更有意义的一个处理方向。

参考资料

[1] 高分系列卫星介绍与下载教程：https://www.shangyexinzhi.com/article/4472540.html
[2] Gdal官方文档：https://gdal.org/drivers/raster/cog.html#raster-cog
[3] 使用python GDAL生成COG(Cloud Optimized GeoTIFF):https://blog.csdn.net/Neuromancerr/article/details/123185553
[4] QGIS仓库：https://github.com/qgis/QGIS
[5] ENVI操作：GF2影像全色与多光谱融合：https://blog.csdn.net/weixin_45905741/article/details/138153273