为什么要干这件事
甲方来一句"成果统一交 CGCS2000 3 度带",你打开文件夹一看:87 个 Excel、40 多个 shp,还有几个同事从 CAD 导出来的。手动一个个另存为,一天就交代进去了,而且第 60 个的时候你一定会点错一次带号——这种错还查不出来,因为图形看着都在。
这篇给你一段能直接跑的脚本,把整个文件夹一次转完,顺手统一小数位和字段编码。
适用
- 会一点 Python(能改变量、会运行 .py 就够了),不会也能照抄
- 手上有成堆的 xlsx / csv / shp 要换坐标系
- 目标坐标系明确(问不清楚就别开始转,见下文「常见坑」第 4 条)
准备
pip install pyproj pandas openpyxl geopandas
只转 Excel 表格的话,geopandas 可以不装——它体积大,装起来慢。脚本里做了懒加载,遇到 shp 才会导入。
本文所有代码都在 Python 3.13 + pyproj 3.7 + geopandas 1.1 上实测跑通,输出结果就是文中贴的数字。
步骤
1. 先搞清楚「从哪来、到哪去」
这是整件事最容易出错、也最省时间的一步。别信文件名,别信同事口述,用代码把坐标系的官方名字打出来看:
from pyproj import CRSprint(CRS.from_epsg(4490).name)print(CRS.from_epsg(4547).name)print(CRS.from_epsg(4526).name)
输出:
China Geodetic Coordinate System 2000CGCS2000 / 3-degree Gauss-Kruger CM 114ECGCS2000 / 3-degree Gauss-Kruger zone 38
看清楚了——4547 和 4526 是同一个投影,区别只在 Y 值带不带带号前缀。选错了,坐标会差 3800 万米。
国内常用的几个记住就够用:
记不住也没关系,对话里那个 EPSG 速查器选一选就出码,选完直接把那行 Transformer.from_crs(...) 复制走。
2. 建转换器:一次定义,全程复用
from pyproj import Transformertf = Transformer.from_crs(”EPSG:4490”, ”EPSG:4547”, always_xy=True)x, y = tf.transform(114.5, 30.5)print(round(x, 3), round(y, 3))# 547999.761 3375648.033
两个细节:
always_xy=True 一定要写。 不写的话 pyproj 按 EPSG 官方轴序收参数,4490 的官方轴序是「纬度在前」,你传 (114.5, 30.5) 它当成纬度 114.5——直接算飞。- Transformer 建在循环外面。 它初始化时要查参数表,放循环里每个文件建一次,几百个文件能慢出感觉来。
3. 让脚本自己去翻文件夹
from pathlib import PathIN = Path(r”D:\成果\原始”)OUT = Path(r”D:\成果\转换后”)OUT.mkdir(parents=True, exist_ok=True)for f in IN.rglob(”*.xlsx”): print(f.relative_to(IN))
glob 只翻当前一层,rglob 连子文件夹一起翻。路径前面的 r 别丢,否则 \成果 里的 \成 会被当转义字符。
输出目录一定跟输入目录分开。原地覆盖翻车了没有后悔药。
4. 表格走 Transformer,矢量走 to_crs
表格(整列一次性传,比逐行 apply 快一两个量级):
x, y = tf.transform(df[”经度”].values, df[”纬度”].values)df[”X”] = pd.Series(x).round(3)df[”Y”] = pd.Series(y).round(3)
矢量(geopandas 会顺手把 .prj 写对):
import geopandas as gpdgpd.read_file(f).to_crs(”EPSG:4547”).to_file(out, encoding=”utf-8”)
实测一条:(114.5, 30.5) → (547999.761, 3375648.033)。
5. 抽点核验,再交付
chk = gpd.read_file(out)print(chk.crs.name)print(chk.head(3).geometry.tolist())
挑首、中、尾三个点,跟已知控制点或者天地图比一比:
- 差 几十厘米 → 多半是大地基准选错了(西安 80 当成了 CGCS2000)
6. 完整脚本
存成 convert.py,改 CONFIG 里的六行,然后 python convert.py。xlsx / xls / csv / shp / gpkg / geojson 都吃。
# -*- coding: utf-8 -*-”””批量坐标转换:xlsx / csv / shp 一把梭。改完 CONFIG 直接 python convert.py”””from pathlib import Pathimport pandas as pdfrom pyproj import CRS, TransformerCONFIG = { ”in_dir”: r”D:\成果\原始”, ”out_dir”: r”D:\成果\转换后”, ”src”: ”EPSG:4490”,# 源坐标系 ”dst”: ”EPSG:4547”,# 目标坐标系 ”lon_col”: ”经度”,# 表格里的经度列名 ”lat_col”: ”纬度”,# 表格里的纬度列名 ”recursive”: True,# 是否连子文件夹一起翻 ”digits”: 3,# 输出保留小数位}TABLE_EXT = {”.xlsx”, ”.xls”, ”.csv”}VECTOR_EXT = {”.shp”, ”.gpkg”, ”.geojson”, ”.json”}def log(tag, msg): print(f”[{tag}] {msg}”)def convert_table(f, out_dir, tf, cfg): df = pd.read_csv(f, encoding=”utf-8-sig”) if f.suffix.lower() == ”.csv” else pd.read_excel(f) lon, lat = cfg[”lon_col”], cfg[”lat_col”] if lon not in df.columns or lat not in df.columns: log(”跳过”, f”{f.name} 找不到 {lon} / {lat} 列,实际列:{list(df.columns)}”) return bad = df[[lon, lat]].isna().any(axis=1).sum() if bad: log(”警告”, f”{f.name} 有 {bad} 行坐标为空,将输出为空值”) x, y = tf.transform(df[lon].values, df[lat].values) df[”X”] = pd.Series(x).round(cfg[”digits”]) df[”Y”] = pd.Series(y).round(cfg[”digits”]) out = out_dir / f”{f.stem}_{cfg['dst'].split(':')[-1]}{f.suffix}” if f.suffix.lower() == ”.csv”: df.to_csv(out, index=False, encoding=”utf-8-sig”) else: df.to_excel(out, index=False) log(”完成”, f”{f.name} -> {out.name} {len(df)} 行”)def convert_vector(f, out_dir, cfg): import geopandas as gpd gdf = gpd.read_file(f) if gdf.crs is None: log(”补 CRS”, f”{f.name} 缺 .prj,按 {cfg['src']} 处理”) gdf = gdf.set_crs(cfg[”src”]) out = out_dir / f”{f.stem}_{cfg['dst'].split(':')[-1]}{f.suffix}” gdf.to_crs(cfg[”dst”]).to_file(out, encoding=”utf-8”) log(”完成”, f”{f.name} {gdf.crs.to_string()} -> {cfg['dst']} {len(gdf)} 个要素”)def main(cfg=CONFIG): in_dir, out_dir = Path(cfg[”in_dir”]), Path(cfg[”out_dir”]) out_dir.mkdir(parents=True, exist_ok=True) log(”源”, CRS.from_user_input(cfg[”src”]).name) log(”目标”, CRS.from_user_input(cfg[”dst”]).name) tf = Transformer.from_crs(cfg[”src”], cfg[”dst”], always_xy=True) files = (in_dir.rglob(”*”) if cfg[”recursive”] else in_dir.glob(”*”)) n = 0 for f in sorted(files): ext = f.suffix.lower() if not f.is_file() or ext not in TABLE_EXT | VECTOR_EXT: continue try: if ext in TABLE_EXT: convert_table(f, out_dir, tf, cfg) else: convert_vector(f, out_dir, cfg) n += 1 except Exception as e: log(”失败”, f”{f.name}: {e}”) log(”汇总”, f”共处理 {n} 个文件,输出在 {out_dir}”)if __name__ == ”__main__”: main()
实际跑起来长这样(真实输出):
[源] China Geodetic Coordinate System 2000[目标] CGCS2000 / 3-degree Gauss-Kruger CM 114E[完成] A区隐患点.xlsx -> A区隐患点_4547.xlsx 3 行[完成] B区隐患点.xlsx -> B区隐患点_4547.xlsx 3 行[完成] C区.csv -> C区_4547.csv 3 行[跳过] 列名不对.xlsx 找不到 经度 / 纬度 列,实际列:['点号', '名称', 'X坐标', '纬度'][完成] 地块A.shp EPSG:4490 -> EPSG:4547 2 个要素[完成] 地块B.shp EPSG:4490 -> EPSG:4547 2 个要素[汇总] 共处理 6 个文件,输出在 D:\成果\转换后
注意那条「跳过」——列名不匹配的文件不会静默出错,而是明确告诉你实际列名是什么。批量作业里,能报错比能跑完更重要。
常见坑
1.always_xy不写,坐标直接飞。pyproj 遵守 EPSG 官方轴序,4490 和 4326 的官方轴序都是「纬度在前」。不加 always_xy=True,你传 (经度, 纬度) 它按 (纬度, 经度) 读。这个坑每个人都要踩一次,早点知道能省一晚上。
2. 带号前缀:38500000 和 500000。EPSG:4526(zone 38)的 Y 值是 38547999.761,EPSG:4547(CM 114E)是 547999.761。同一个点,同一个投影,就差前面那个 38。甲方要哪种一定问清楚,不确定就看他给的样例数据 Y 值有几位数。
3. WGS84 转 CGCS2000,pyproj 会「原样返回」。实测 EPSG:4326 → EPSG:4490 传入 (114.5, 30.5),输出还是 (114.5, 30.5)。因为 EPSG 库里没给这两者定义变换参数,pyproj 按空变换处理。这不代表两者等价——它们参考框架和历元不同,实际存在厘米级差异。日常项目当同一套用问题不大,但高精度控制、变形监测这类场合,成果里必须写清楚到底是哪个,别糊弄过去。
4. 地方独立坐标系,没有参数就是转不了。"城建坐标""某某市独立坐标系"这类,本质是在国家坐标系上做了平移旋转缩放。没有官方给的四参数或七参数,谁也变不出来。网上搜来的参数绝对不能用——参数是分区域的,用错了偏移几米还没有任何报错。正规做法是向当地测绘主管部门申请。拿到参数后用 proj 的 pipeline 写:
Transformer.from_pipeline( ”+proj=pipeline +step +proj=helmert +x=... +y=... +z=... ” ”+rx=... +ry=... +rz=... +s=... +convention=coordinate_frame”)
5. shp 缺.prj,或者属性表中文乱码。.shp 其实是一组文件,.prj 存坐标系、.dbf 存属性、.cpg 存编码。只拷了 .shp 过来,geopandas 读出来 crs 是 None——脚本里已经做了兜底(按 src 处理并打印提示),但你得确认这个兜底假设是对的。中文乱码就在读写时显式指定 encoding="utf-8" 或 "gbk",跟源数据保持一致。
6. 高程不会跟着平面一起转。平面坐标转换动的是 X、Y。1985 国家高程基准是正常高,GNSS 直接测的是大地高,两者差一个高程异常,得靠似大地水准面模型改正——比如内蒙古今年 8 月刚通过验收的那套全域 ±5 厘米精度的 GNSS 似大地水准面模型,走的就是这条路。别指望 Transformer 帮你处理 Z。
7. Excel 里的坐标被存成了文本。表格里坐标列前面带个绿三角,或者显示成 1.145E+02,读进 pandas 就是字符串,transform 会报类型错误。转换前加一行:
df[lon] = pd.to_numeric(df[lon], errors=”coerce”)df[lat] = pd.to_numeric(df[lat], errors=”coerce”)
errors="coerce" 会把转不了的变成 NaN,配合脚本里那个空值统计,一眼就知道有多少行有问题。