转换shapefile投影：将Shapefile文件中的字段中的投影坐标系转换为经纬度坐标系

         这次的需求是将投影为Gauss_Kruger，投影坐标系为CGCS2000_3_Degree_GK_CM_108E的shapefile（shp）文件中某些字段的值进行投影坐标系的转换，将投影坐标转换为经纬度坐标（WGS84）

1.算法思路

        （1）使用相应的地理信息系统（GIS）库读取Shapefile文件的几何信息和属性表。

        （2）利用Pyproj库进行投影坐标系的转换，将投影坐标转换为经纬度坐标。

        （3）更新Shapefile文件中对应字段的值为转换后的经纬度坐标。

        （4）保存更新后的Shapefile文件

2.算法步骤

        （1）读取Shapefile文件并获取投影坐标系信息。

        （2）定义投影坐标系和经纬度坐标系之间的转换函数

        （3）对Shapefile中的每个要素进行循环处理：

                读取该要素的投影坐标

                转换为经纬度坐标

                更新要素的属性表中的经纬度字段值

        （4）保存更新后的Shapefile文件

3.Python代码实现

import shapefile
from pyproj import Proj, Transformer

# 读取Shapefile文件
shp_path = "shp/cut.shp"
sf = shapefile.Reader(shp_path)

# 读取投影信息
prj_file = open("shp/cut.prj", "r")
prj_content = prj_file.read()
prj_file.close()

# 构建投影对象
projection = Proj(prj_content)

# 定义投影坐标系和经纬度坐标系之间的转换函数
def project_to_latlon(x, y):
    transformer = Transformer.from_proj(projection, "epsg:4326")
    lon, lat = transformer.transform(x, y)
    return lon, lat

# 获取要素字段名
fields = [field[0] for field in sf.fields[1:]]  # 排除第一个字段 FID

# 确定要素字段中经纬度字段的索引(我需要转换的字段名为center_x和center_y)
center_x_index = fields.index("center_x")
center_y_index = fields.index("center_y")

# 更新Shapefile中的字段值为经纬度坐标
records = sf.records()
for record in records:
    center_x = record[center_x_index]
    center_y = record[center_y_index]
    lon, lat = project_to_latlon(center_x, center_y)
    # 更新经纬度字段值
    record[center_x_index] = lon
    record[center_y_index] = lat

# 保存更新后的Shapefile文件
with shapefile.Writer("transform_cut.shp") as w:
    # 添加字段
    w.fields = sf.fields[1:]
    # 添加几何信息和更新后的记录
    for i, shape in enumerate(sf.shapes()):
        w.record(*records[i])
        w.shape(shape)