cgcs2000坐标系的epsg计算
cgcs2000坐标系的epsg计算
ytkz之前写了一篇关于wgs84影像转为cgcs2000的影像的文章：《wgs84影像直接转换为cgcs2000影像》
在里面我写道：诸如遥感数据处理的东西，其实是没有技术含量的，本质上都是调用gdal的东西，你把gdal的各个接口弄明白了就掌握这个数据处理的技能，唯手熟尔。真正有难度的是算法、加速处理速度。我其实很久不做数据数据处理的东西了，因为我有足够的代码积累，我建立了一个代码库，需要哪个模块，我直接搜索。
cscs2000坐标系的EPSG代码计算
cscs2000坐标系其实和wgs84一样，分为两种：地理坐标系和投影坐标系。
用大白话来区分地理和投影，是看它们的单位。
单位若是度分秒，就是地理坐标系；单位是米，就是投影坐标系。
除此之外，我们还要获取坐标系的EPSG号码。
这个参数，我们需要填写到gdal的Warp进行重投影。
wgs84地理坐标系的EPSG号码4326，cgcs2000地理坐标系的EPSG号码4490。
如果你是想转换为cgcs2000投影坐标系，还得继续计算。
投影坐标系分为3度带，和6度带。
以3度带为例，不同的带号的计算规则如下：
例如，经度为105°E，计算过程是 105 / 3 = 35，得到带号 35。如果有小数点，则进行四舍五入。
我们继续计算 EPSG = 4513 + (带号 - 25)
代码如下：
def get_cgcs2000_3degree_epsg(lon):
"""
根据经度计算 CGCS2000 3度带投影的 EPSG 代码。
计算公式: Zone = int(round(lon / 3.0))
EPSG = 4513 + (Zone - 25)
"""
zone = int(round(lon / 3.0))
epsg = 4513 + (zone - 25)
return epsg, zone6 度带适用于对精度要求相对较低，或者需要覆盖更大范围的项目。计算6度带的epsg号码的代码如下：
def get_cgcs2000_6degree_epsg(lon):
if lon < 75 or lon > 135:
# 超出中国主要范围，需要谨慎处理或使用地理坐标系
pass
zone_6 = int((lon + 1.5) / 6) + 1 # 带号
epsg = 4474 + zone_6
zone_6_simple = int((lon + 6) / 6) + 13
epsg_simple = 4490 + zone_6_simple
return epsg_simple, zone_6_simpleGDAL 的核心：重投影操作 (gdal.Warp)、
一旦我们确定了目标 EPSG 码（无论是地理坐标系 4490，还是投影坐标系如 4528），重投影的操作在 GDAL 中变得异常简单。
核心函数是 gdal.Warp()，它负责将一个数据集（数据源）重投影到另一个坐标系。
from osgeo import gdal
def warp_wgs84_to_cgcs2000(input_file, output_file, target_epsg):
"""
使用 gdal.Warp 实现 WGS84 到 CGCS2000 的重投影。
:param input_file: WGS84 格式的输入影像路径 (默认 EPSG:4326)
:param output_file: CGCS2000 格式的输出影像路径
:param target_epsg: 目标 CGCS2000 的 EPSG 码 (如 4490, 4528 等)
"""
print(f"开始重投影：从 WGS84 (EPSG:4326) 到 CGCS2000 (EPSG:{target_epsg})")
try:
# 调用 gdal.Warp 函数
gdal.Warp(
destNameOrDestDS=output_file, # 输出文件路径
srcDSOrSrcDSTab=input_file, # 输入文件路径
srcSRS='EPSG:4326', # 明确指定源坐标系（WGS84）
dstSRS=f'EPSG:{target_epsg}', # 目标坐标系
resampleAlg=gdal.GRA_Bilinear, # 重采样方法（双线性插值，适用于影像）
creationOptions=['COMPRESS=LZW', 'BIGTIFF=IF_NEEDED'] # 输出文件参数
)
print(f"重投影成功！新影像已保存至: {output_file}")
except Exception as e:
print(f"GDAL 重投影失败: {e}")
# 1. 假设影像中心经度为 105.0
center_lon = 105.0
# 2. 计算 3 度带的 EPSG 码
target_epsg, zone = get_cgcs2000_3degree_epsg(center_lon)
# 3. 执行重投影
input_path = "path/to/your/wgs84_image.tif"
output_path = f"path/to/output/cgcs2000_{target_epsg}.tif"
# warp_wgs84_to_cgcs2000(input_path, output_path, target_epsg)总结：技术含量在哪？
正如文章开头所言：诸如遥感数据处理的东西，其实是没有技术含量的，本质上都是调用 gdal 的东西，你把 gdal 的各个接口弄明白了就掌握这个数据处理的技能，唯手熟尔。
真正的技术含量体现在以下几个方面：
1.算法设计与优化
2.性能加速与资源管理
3.代码积累与模块化复用
4.实际应用中的问题诊断与创新
希望这篇文章能帮你快速上手 WGS84 到 CGCS2000 的转换，有问题欢迎讨论！
