使用GDAL库按1000x1000米格网分幅处理卫图

概述

最近一位小伙伴私信我，说他在处理卫星影像时遇到一个需求：需要将大范围的栅格影像按固定网格（比如 1000x1000 米）分割成小块瓦片，确保接边无缝，还要保存为 PNG 格式方便查看。他希望有一个简单易用的工具，最好不用装一堆环境。于是，我写了一个 Python 脚本 split_raster_into_tiles.py，用 GDAL 实现了这个功能，并用 Pyinstaller 打包成了一个独立的 .exe 文件。这篇文章就来介绍一下这个工具的用法，希望能帮到有类似需求的朋友！

需求背景

小伙伴的需求很明确：
输入：一张大范围的卫星影像（比如 TIFF 格式）。
处理：按 1000x1000 米网格分割，边缘无缝对接。
输出：生成多个小瓦片，保存为 PNG 格式（方便查看），同时保留地理信息。
额外要求：最好是个独立程序，不用折腾 。

基于这个需求，我开发了 split_raster_into_tiles.py，并用 Pyinstaller 打包成了 split_raster_into_tiles.exe。下面就来详细介绍怎么用它。

功能简介

split_raster_into_tiles.exe 是一个命令行工具，主要功能包括：
读取栅格影像（如 TIFF）。
按指定范围（默认 1000 米）分割成小瓦片。

这里要说明一下，1000米是针对的是投影坐标系的影像，单位是米。如果你的影像是地理坐标系，这个参数的单位是度数。

确保瓦片之间无缝衔接，保留地理参考信息。
输出 PNG 格式图像（用于可视化）和 TIFF 格式文件（带地理信息）。
这个工具特别适合遥感影像预处理、地图瓦片生成等场景，尤其是对非技术用户来说，免去了配置环境的麻烦。

前置条件

操作系统：Windows（其他系统需重新打包）。
输入文件：支持 GDAL 可读取的栅格格式（推荐 TIFF）。
磁盘空间：确保输出目录有足够空间存放瓦片。
打包后无需安装 Python、GDAL 或其他依赖，所有库都已集成到 .exe 中。

  1. 准备工作
    把 split_raster_into_tiles.exe 放到一个工作目录，比如 C:\GIS_Tools。
    准备好输入影像，比如 input.tif。
    创建一个输出目录，比如 C:\output，用来保存生成的瓦片。

  2. 运行程序
    打开命令提示符（Windows 下按 Win + R，输入 cmd，回车）。
    切换到工具所在目录：
    cd C:\GIS_Tools运行程序：
    split_raster_into_tiles.exe3. 输入参数
    程序启动后会依次提示你输入：
    输入影像路径：键入影像文件的完整路径，比如 C:\data\input.tif。
    输入保存路径：键入输出目录，比如 C:\output。
    分幅范围：输入分割网格的大小（单位：米），比如 1000。
    输入完后按回车，程序就会开始处理。

  3. 输出结果
    处理完成后，输出目录中会生成一堆文件：
    tile_x_y.tif：带地理信息的 TIFF 瓦片。
    tile_x_y.png：对应的 PNG 图像。 控制台会显示完成信息，例如：

注意事项

路径问题：路径要写全，避免中文路径（可能会出编码问题）。
影像格式：输入文件要是 GDAL 支持的格式，TIFF 最稳。
分幅范围：默认 1000 米，可以改，比如 500 米，但太小可能会生成很多瓦片。
内存和空间：大影像处理时，确保电脑内存够用，磁盘空间充足。
错误排查：如果提示“无法打开输入影像”，检查文件路径或文件是否损坏。

总结

从一个小伙伴的需求出发，我开发并打包了 split_raster_into_tiles.exe，希望它能帮到更多处理栅格影像的朋友。工具虽小，但解决了实际问题：无缝分割、PNG 输出、免环境配置。如果你也有类似需求，不妨试试这个 .exe，或者直接告诉我你的想法，我再帮你调整代码。
技术是为需求服务的，有问题欢迎随时交流！

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 2025/3/5 10:59 
# @File : split_raster_into_tiles.py 
'''
使用GDAL库按1000x1000米格网分幅处理卫图的Python函数代码，确保影像接边无缝并保存为PNG格式
'''

import os
import numpy as np
from osgeo import gdal
from PIL import Image, ImageDraw

def split_raster_into_tiles(input_raster_path, output_dir):
    # 打开原始影像
    ds = gdal.Open(input_raster_path)
    if ds is None:
        raise ValueError("无法打开输入影像")

    # 获取影像基本信息
    geotransform = ds.GetGeoTransform()
    proj = ds.GetProjection()
    x_size = ds.RasterXSize
    y_size = ds.RasterYSize
    num_bands = ds.RasterCount
    px_width = abs(geotransform[1])  # 像素宽度（米）
    px_height = abs(geotransform[5])  # 像素高度（米）

    # 计算每个分幅的像素尺寸
    tile_pixel_width = int(round(1000 / px_width))
    tile_pixel_height = int(round(1000 / px_height))

    # 计算分幅数量
    num_x_tiles = int(np.ceil(x_size / tile_pixel_width))
    num_y_tiles = int(np.ceil(y_size / tile_pixel_height))

    # 创建输出目录
    os.makedirs(output_dir, exist_ok=True)

    # 遍历所有分幅
    for x in range(num_x_tiles):
        for y in range(num_y_tiles):
            # 计算当前分幅的像素范围
            x_offset = x * tile_pixel_width
            y_offset = y * tile_pixel_height

            # 计算实际读取尺寸
            read_x = min(tile_pixel_width, x_size - x_offset)
            read_y = min(tile_pixel_height, y_size - y_offset)

            # 初始化数据数组
            tile_data = np.zeros((tile_pixel_height, tile_pixel_width, num_bands), dtype=np.uint8)

            # 读取每个波段数据
            for b in range(num_bands):
                band = ds.GetRasterBand(b + 1)
                data = band.ReadAsArray(x_offset, y_offset, read_x, read_y)

                # 处理边缘分幅
                if data is None:
                    continue

                # 填充数据到分幅数组
                if data.ndim == 2:
                    tile_data[:read_y, :read_x, b] = data
                else:
                    tile_data[:read_y, :read_x, b] = data[..., 0]

            # 计算当前分幅的地理坐标
            ulx = geotransform[0] + x_offset * geotransform[1]
            uly = geotransform[3] + y_offset * geotransform[5]

            # 创建输出文件
            output_path = os.path.join(output_dir, f'tile_{x}_{y}.tif')
            driver = gdal.GetDriverByName('GTiff')
            target_ds = driver.Create(output_path,
                                      tile_pixel_width,
                                      tile_pixel_height,
                                      num_bands,
                                      gdal.GDT_Byte)

            # 设置地理参考
            new_geotransform = (
                ulx,  # 左上角x坐标
                geotransform[1],  # 像素宽度
                geotransform[2],  # 旋转参数
                uly,  # 左上角y坐标
                geotransform[4],  # 旋转参数
                geotransform[5]  # 像素高度
            )

            target_ds.SetGeoTransform(new_geotransform)
            target_ds.SetProjection(proj)

            # 写入数据
            for b in range(num_bands):
                target_ds.GetRasterBand(b + 1).WriteArray(tile_data[..., b])


            image = Image.fromarray(tile_data, mode="RGB")
            target_path = os.path.join(output_dir, f'tile_{x}_{y}.png')
            image.save(target_path, 'PNG')  # 保存

            target_ds.FlushCache()
            target_ds = None

    ds = None
    print(f"分幅完成，共生成 {num_x_tiles}x{num_y_tiles}={num_x_tiles * num_y_tiles} 个分幅")



if __name__ == '__main__':
    file = r'F:\out1_qgis.tif'
    output_dir = r'F:\output'
    split_raster_into_tiles(file, output_dir)

获取工具

公众号回复分幅使用方法

结果对比

输入影像如下：

image-20250306101719172

输出结果如下：

image-20250306101634629