使用 Python 将矢量数据重投影到 Krasovsky_1940_Albers 坐标系

ytkz2025-02-272025-02-27

在地理信息系统（GIS）中，坐标系转换是常见的操作，尤其在处理历史数据或特定区域测绘数据时。本文将介绍如何使用 Python 和 GDAL/OGR 库，将矢量文件（如 Shapefile）从任意坐标系重投影到基于 Krasovsky 1940 椭球的 Albers 等面积圆锥投影坐标系（简称 Krasovsky_1940_Albers）。该方法适用于需要处理中国或其他使用 Krasovsky 1940 椭球数据的场景，具有较强的实用性。

背景

Krasovsky_1940_Albers 是一种基于 Krasovsky 1940 椭球的投影坐标系，结合 Albers 等面积圆锥投影的特点，广泛用于保持面积精确的测绘任务。Krasovsky 1940 椭球（长半轴 6,378,245 米，扁率 1/298.3）曾是苏联及中国早期测绘的基础，而 Albers 投影通过指定标准纬线和中央经线，适用于中纬度地区的长条形区域。

由于 Krasovsky_1940_Albers 并非 EPSG 数据库中的标准坐标系，其具体参数（如标准纬线、中央经线）需根据数据来源定义。本文提供一个通用的实现方案，并以中国常用参数为例。

Krasovsky_1940_Albers投影介绍:

实现步骤

1. 环境准备

依赖库：安装 GDAL（包含 OGR 和 OSR 模块），可通过以下命令安装：
输入数据：一个矢量文件（如 Shapefile），需明确其源坐标系
输出目标：重投影后的 Shapefile，坐标系为 Krasovsky_1940_Albers。

2. 代码实现

以下是一个完整的 Python 脚本，用于将矢量文件重投影到 Krasovsky_1940_Albers：

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 2025/2/27 20:28 
# @File : transfer_other_project.py 
"""
将矢量文件重投影到 Krasovsky_1940_Albers 坐标系
依赖库：osgeo (gdal, ogr, osr)
"""

import os
from osgeo import ogr, osr

def reproject_to_krasovsky_albers(input_shp_path, output_shp_path, source_epsg=4326):
    """
    将输入矢量文件重投影到 Krasovsky_1940_Albers 坐标系
    Args:
        input_shp_path (str): 输入Shapefile路径
        output_shp_path (str): 输出Shapefile路径
        source_epsg (int): 输入数据的EPSG代码，默认为4326（WGS84）
    """
    # 定义 Krasovsky_1940_Albers 投影参数（需根据实际情况调整）
    albers_proj4 = "+proj=aea +ellps=krass +lat_1=25 +lat_2=47 +lon_0=105 +x_0=0 +y_0=0  +units=m +no_defs"

    # 打开输入Shapefile
    driver = ogr.GetDriverByName("ESRI Shapefile")
    input_ds = ogr.Open(input_shp_path)
    if input_ds is None:
        raise ValueError(f"无法打开输入文件: {input_shp_path}")
    input_layer = input_ds.GetLayer()

    # 定义源坐标系
    source_srs = osr.SpatialReference()
    source_srs.ImportFromEPSG(source_epsg)  # 默认WGS84，可以根据输入数据的实际EPSG修改
    # 设置输入空间参考的坐标轴映射策略为传统GIS顺序（经度，纬度）
    source_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)


    # 定义目标坐标系（Krasovsky_1940_Albers）
    target_srs = osr.SpatialReference()
    target_srs.ImportFromProj4(albers_proj4)
    # # 设置目标空间参考的坐标轴映射策略为传统GIS顺序
    target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    # 创建坐标转换
    coord_transform = osr.CoordinateTransformation(source_srs, target_srs)

    # 创建输出Shapefile
    if os.path.exists(output_shp_path):
        driver.DeleteDataSource(output_shp_path)
    output_ds = driver.CreateDataSource(output_shp_path)
    output_layer = output_ds.CreateLayer("reprojected", srs=target_srs,
                                        geom_type=input_layer.GetLayerDefn().GetGeomType())

    # 复制字段定义
    input_defn = input_layer.GetLayerDefn()
    for i in range(input_defn.GetFieldCount()):
        output_layer.CreateField(input_defn.GetFieldDefn(i))

    # 遍历并重投影每个要素
    for feature in input_layer:
        geom = feature.GetGeometryRef()
        geom.Transform(coord_transform)  # 重投影到 Krasovsky_1940_Albers

        # 创建新要素并写入
        new_feature = ogr.Feature(output_layer.GetLayerDefn())
        new_feature.SetGeometry(geom)
        for i in range(input_defn.GetFieldCount()):
            new_feature.SetField(i, feature.GetField(i))
        output_layer.CreateFeature(new_feature)
        new_feature = None

    # 清理
    input_ds = None
    output_ds = None
    print(f"重投影完成，输出文件: {output_shp_path}")

if __name__ == "__main__":

    input_shp = r"C:\roi.shp" 
    output_shp = r"C:\roi_7024.shp" 

    source_epsg = 4547     # 输入数据的EPSG代码（这里假设是WGS84）

    reproject_to_krasovsky_albers(input_shp, output_shp, source_epsg)

    # ogr2ogr -t_srs "+proj=aea +ellps=krass +lat_1=25 +lat_2=47 +lon_0=105 +x_0=0 +y_0=0 +units=m +no_defs" output.shp

3. 运行与验证

运行脚本：将 input_shp 和 output_shp 替换为实际路径，设置正确的 source_epsg，然后执行脚本。
验证结果：使用 QGIS 或 ArcGIS 打开输出文件，检查坐标单位（米）和投影参数是否符合预期。

4.实用性提升

支持多个文件重投影，只需稍作修改：

input_dir = r"C:\path\to\shp_folder"
for shp in os.listdir(input_dir):
    if shp.endswith(".shp"):
        input_shp = os.path.join(input_dir, shp)
        output_shp = input_shp.replace(".shp", "_albers.shp")
        reproject_to_krasovsky_albers(input_shp, output_shp, source_epsg=4547)

注意事项

源坐标系确认：脚本假设输入为 EPSG:4547
投影参数验证：错误的 lat_1, lat_2, 或 lon_0 会导致投影失真，务必与项目要求一致。

性能优化：对于大型数据集，可考虑使用 ogr2ogr 命令行工具：

ogr2ogr -t_srs "+proj=aea +ellps=krass +lat_1=25 +lat_2=47 +lon_0=105 +x_0=0 +y_0=0 +units=m +no_defs" output.shp input.shp

结论

通过本文提供的脚本，你可以轻松将矢量数据重投影到 Krasovsky_1940_Albers 坐标系。代码灵活支持参数调整和编码处理，适用于多种实际场景。无论是单一文件转换还是批量处理，该方法都能显著提升工作效率。建议根据你的具体数据调整投影参数，并在转换后使用 GIS 软件验证结果。