【GIS教程】使用Arcpy+GDAL批量投影并批量生成COG文件-附完整代码

项目场景：

之前写过一个帖子将tif转为COG【GIS教程】使用GDAL-Python将tif转为COG并在ArcGIS Js前端加载-附完整代码，现在有一个需求：需要批量处理多张栅格数据，因此我将其改成了一个批处理脚本。

由于最后的COG是需要在web端进行展示，因此在转COG之前还需要对栅格数据进行一个批量投影。我们可以编写Arcpy脚本进行批量投影。

在web端展示的GIS数据我们一般使用的坐标系为Web墨卡托（EPSG:3857） ，是 Web 地图的一种常用切片格式。

实现代码：

一、 使用Arcpy进行批量投影

（1）环境配置

添加python环境，在IDEA或Pycharm中，【File】-【Project Structure】-【SDKS】中，添加python环境。选择ArcGISPro安装目录下的arcpy环境 (python-3)，路径如下：D:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe

ps：如果是ArcMap的arcpy的环境路径为 (python-2)：
D:\Program Files\ArcGIS\Python2.7\ArcGIS10.8\python.exe

（2）批处理脚本

使用函数arcpy.ProjectRaster_management：

import os
import arcpy

# 栅格数据批量投影
def batch_raster_prj(raster_path, save_path, reference_file):
	"""
	raster_path ：需要批量投影的栅格文件的文件夹路径
	save_path ： 输出投影结果的文件夹路径
	reference_file ： 一个具有目标投影的tif文件路径，作为投影模板文件: "./data/reference.tif"
	
	"""
    try:
        if not os.path.exists(save_path):
            os.mkdir(save_path)
        for file_list in os.listdir(raster_path):
            if file_list.endswith('.tif') or file_list.endswith('.tiff'):
                raster = os.path.join(raster_path, file_list)# 待投影文件
                print(raster)
                outPath = os.path.join(save_path, file_list.split('.')[0]+'_prj.tif')# 投影结果文件
                print(outPath)
                coordinate_system = arcpy.Describe(reference_file).spatialReference
                arcpy.ProjectRaster_management(raster, outPath, coordinate_system, "BILINEAR", "", "", "#", "#")
    except arcpy.ExecuteError:
        print(arcpy.GetMessages())

#执行函数
batch_raster_prj(raster_path, save_path, reference_file)

当如果没有坐标模板文件 ，以下提供几种替代方法：
①使用具有水平坐标系的 WKT2 字符串创建 SpatialReference 对象。

spatialReference = """
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Mercator_Auxiliary_Sphere"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
PARAMETER["Standard_Parallel_1",0.0],
PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
"""

coordinate_system = arcpy.SpatialReference(text=spatialReference)

②使用坐标系代码（或 WKID）创建 SpatialReference 对象。

coordinate_system = arcpy.SpatialReference(3857)

（3）Tips：报错解决（未遇到就跳过）

在进行以上脚本批处理时候我遇到了一个报错：

报错原因为我的文件名中有空格，但是我在ArcgisPro的工具箱中执行栅格投影工具是不会报错这个问题的 。

解决办法为：批量修改文件名

代码如下：

import os

#批量替换文件名，移除文件名中空格
def remove_spaces_in_filenames(directory):
    """
    param directory: The path to the directory to process.
    """
    for root, dirs, files in os.walk(directory):
        for filename in files:
            if ' ' in filename:
                new_filename = filename.replace(' ', '')
                old_file = os.path.join(root, filename)
                new_file = os.path.join(root, new_filename)
                os.rename(old_file, new_file)
                print(f'Renamed "{old_file}" to "{new_file}"')

#执行函数
remove_spaces_in_filenames(raster_path)

二、 使用GDAL进行批量COG转换

（1）环境配置

Anaconda新建python环境

conda create --name gdalEnv python=3.10

激活虚拟环境：

activate gdalEnv

GDAL环境安装，直接pip会有很多报错，推荐使用whl离线安装。

conda install GDAL-3.6.4-cp310-cp310-win_amd64.whl

（2）批处理代码

直接上批处理代码：

from osgeo import gdal, gdalconst
import os

'''获取相应统计数据'''
def compute_stat(srcDS):
    outBand = srcDS.GetRasterBand(1)
    stat = outBand.ComputeStatistics(True)

    outBand = None
    stat = None
    del outBand
    del stat


# 定义输入和输出文件夹路径
input_folder = r'C:\yourpath\Desktop\web墨卡托投影'  # 替换为您的输入文件夹路径
output_folder = r'C:\yourpath\Desktop\COG'  # 替换为您的输出文件夹路径

# 确保输出文件夹存在
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

# 遍历输入文件夹中的所有文件
for filename in os.listdir(input_folder):
    if filename.endswith('.tif') or filename.endswith('.tiff'):
        input_path = os.path.join(input_folder, filename)
        output_path = os.path.join(output_folder, "c"+filename)

        # 读取输入文件
        cal_area = gdal.Open(input_path)

        # 将文件转换为COG格式
        desImg = gdal.Translate(output_path, cal_area, format="COG")

        # 计算COG的栅格值
        n_desImg = gdal.Open(output_path)
        compute_stat(n_desImg)

        print(f"Converted {filename} to COG format and computed statistics.")

        # 关闭数据集
        n_desImg = None
        cal_area = None

print("All files have been converted to COG format.")