根据经纬度范围自动进行UTM投影【代码实现】

GISer35

提供一个平时工作中经常使用的代码：根据经纬度范围，自动将数据找到对应的UTM投影，并进行投影变换。

注意：

1. 输入图层需采用地理坐标系（如EPSG:4326）
2. 对应修改输出目录和图层名称。

代码如下

# -*- coding: utf-8 -*-
from qgis.core import (
    QgsProject,
    QgsVectorLayer,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsVectorFileWriter,
    QgsRectangle,
)
import os
import re

# === 1. 获取图层 ===
layer = QgsProject.instance().mapLayersByName("河流t")[0]
print(f"✅ 使用图层：{layer.name()}（CRS: {layer.crs().authid()}）")

# === 2. 计算 UTM 带（WGS84）===
ext = layer.extent()
center_x = (ext.xMinimum() + ext.xMaximum()) / 2
center_y = (ext.yMinimum() + ext.yMaximum()) / 2
utm_zone = int((center_x + 180) // 6) + 1
epsg_code = 32600 + utm_zone if center_y >= 0 else 32700 + utm_zone
utm_crs = QgsCoordinateReferenceSystem(f"EPSG:{epsg_code}")
print(f"🎯 UTM 带：第 {utm_zone} 带 → EPSG:{epsg_code}")

# === 3. 构建坐标变换 ===
transform_context = QgsProject.instance().transformContext()
transform = QgsCoordinateTransform(layer.crs(), utm_crs, transform_context)

# === 4. 准备导出选项 ===
options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "GPKG"
options.fileEncoding = "UTF-8"
options.layerName = f"{layer.name()}_UTM{utm_zone}"

# === 5. 执行导出 ===
output_path = r"D:\output\河流t_UTM48.gpkg"
os.makedirs(os.path.dirname(output_path), exist_ok=True)

error_code, new_filename, new_layername, error_msg = QgsVectorFileWriter.writeAsVectorFormatV3(
    layer=layer,
    fileName=output_path,
    transformContext=transform_context,
    options=options,
)

# === 6. 输出结果 ===
if error_code == QgsVectorFileWriter.NoError:
    print(f"✅ 成功导出：{new_filename}")
    # 自动加载结果图层
    result_layer = QgsVectorLayer(new_filename, f"[UTM{utm_zone}] {layer.name()}", "ogr")
    if result_layer.isValid():
        QgsProject.instance().addMapLayer(result_layer)
        print("🖼️ 已加载结果图层")
else:
    print(f"❌ 导出失败（{error_code}）：{error_msg}")

计算UTM带号原理

对于任意一点的经度（λ，单位：度），其所属 UTM 带号按如下公式计算：
带号 = floor((λ + 180) / 6) + 1

λ 是东经为正、西经为负的十进制度数（例如：北京 ≈ 116.4°，拉萨 ≈ 91.1°，纽约 ≈ −74.0°）；
floor() 表示向下取整（即舍去小数部分）；
+180 是为了将经度范围从 [−180, +180) 平移到 [0, 360)，便于整除；
/6 是因为每带宽 6°；
+1 是因带号从 1 开始编号（不是从 0）

更多算法问题，欢迎留言或联系我们。转载须注明出处。

GISer35

补充批量对所有图层执行UTM投影的代码：

# -*- coding: utf-8 -*-
from qgis.core import (
    QgsProject,
    QgsVectorLayer,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsVectorFileWriter,
)
import os
import re

# === 1. 配置输出目录 ===
output_dir = r"D:\output"
os.makedirs(output_dir, exist_ok=True)

# === 2. 获取所有有效矢量图层 ===
project = QgsProject.instance()
vector_layers = [
    lyr for lyr in project.mapLayers().values()
    if isinstance(lyr, QgsVectorLayer) and lyr.isValid() and lyr.crs().authid() == "EPSG:4326"
]

if not vector_layers:
    print("⚠️ 未找到 EPSG:4326 的有效矢量图层")
else:
    print(f"✅ 发现 {len(vector_layers)} 个待处理图层：")
    for lyr in vector_layers:
        print(f"   • {lyr.name()}")

# === 3. 目标 CRS：UTM Zone 48N (WGS84) ===
target_crs = QgsCoordinateReferenceSystem("EPSG:32648")
if not target_crs.isValid():
    raise Exception("❌ EPSG:32648 坐标系无效，请检查 QGIS CRS 数据库")

transform_context = project.transformContext()

# === 4. 批量处理每个图层 ===
success_count = 0
for layer in vector_layers:
    try:
        print(f"\n🔄 处理：{layer.name()}...")

        # ✅ 关键修复：必须实例化 SaveVectorOptions！
        options = QgsVectorFileWriter.SaveVectorOptions()  # ← 第 16 行正确写法！

        options.driverName = "GPKG"
        options.fileEncoding = "UTF-8"
        options.layerName = f"{layer.name()}_UTM48"

        # 构建输出路径（清理非法字符）
        safe_name = re.sub(r'[^\w\u4e00-\u9fa5\-_]', '_', layer.name())
        output_path = os.path.join(output_dir, f"{safe_name}_UTM48.gpkg")

        # 执行导出（V3 接口：layer, fileName, transformContext, options）
        result = QgsVectorFileWriter.writeAsVectorFormatV3(
            layer=layer,
            fileName=output_path,
            transformContext=transform_context,
            options=options,
        )

        error_code, new_filename, new_layername, error_msg = result

        if error_code == QgsVectorFileWriter.NoError:
            print(f"✅ 成功：{new_filename}")
            # 自动加载结果图层
            result_layer = QgsVectorLayer(new_filename, f"[UTM48] {layer.name()}", "ogr")
            if result_layer.isValid():
                project.addMapLayer(result_layer)
            success_count += 1
        else:
            print(f"❌ 失败：{error_msg}")

    except Exception as e:
        print(f"❌ 异常：{str(e)}")

# === 5. 总结 ===
print(f"\n🎉 批量完成：{success_count}/{len(vector_layers)} 个图层成功")
print(f"📁 输出位置：{output_dir}")