wgs84影像直接转换为cgcs2000影像
wgs84影像直接转换为cgcs2000影像
ytkz最近写遥感的东西比较少，今天写一写如何将wgs84影像直接转换为cgcs2000影像，我们直接上代码就好了。
说白了这东西没什么技术含量，都是调用gdal写好的接口，以后我们再讨论一下wgs84的tif和cgcs2000的tif本质上有什么区别，这块我后面再好好研究一下。
至于wgs84影像直接转换为cgcs2000影像的代码，我直接放后面。此外，我还想再唠上一唠。
诸如遥感数据处理的东西，其实是没有技术含量的，本质上都是调用gdal的东西，你把gdal的各个接口弄明白了就掌握这个数据处理的技能，唯手熟尔。真正有难度的是算法、加速处理速度。我其实很久不做数据数据处理的东西了，因为我有足够的代码积累，我建立了一个代码库，需要哪个模块，我直接搜索。
当然，这些代码库，随着这两年的公众号的文章更新，早已被我耗尽。不然我去年不可能做到高强度的更新。
我自认为部分的技术文章还是写的有一点水平的，哈哈我这种黄婆卖瓜自卖自夸行为真不要脸。
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2025/10/26 21:30
# @File : python tif_wgs84_to_cgcs2000.py
'''
TIF影像坐标系转换工具
支持从WGS84转换到CGCS2000等坐标系
'''
from osgeo import gdal, osr
import os
def get_image_info(input_tif):
"""
获取影像的基本信息和坐标系信息
参数:
input_tif: 输入TIF文件路径
返回:
dict: 包含影像信息的字典
"""
dataset = gdal.Open(input_tif)
if dataset is None:
raise Exception(f"无法打开输入影像文件：{input_tif}")
# 获取投影信息
projection = dataset.GetProjection()
geotransform = dataset.GetGeoTransform()
# 获取影像尺寸
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount
# 获取数据类型
datatype = dataset.GetRasterBand(1).DataType
# 解析坐标系
srs = osr.SpatialReference()
srs.ImportFromWkt(projection)
info = {
'width': cols,
'height': rows,
'bands': bands,
'datatype': datatype,
'projection': projection,
'geotransform': geotransform,
'epsg': srs.GetAttrValue('AUTHORITY', 1) if srs.GetAttrValue('AUTHORITY', 1) else 'Unknown',
'proj_name': srs.GetAttrValue('PROJCS') if srs.GetAttrValue('PROJCS') else srs.GetAttrValue('GEOGCS')
}
dataset = None
return info
def reproject_tif_simple(input_tif, output_tif, target_epsg, resample_method='bilinear'):
"""
使用GDAL Warp进行坐标转换（推荐方法）
参数:
input_tif: 输入TIF文件路径
output_tif: 输出TIF文件路径
target_epsg: 目标坐标系的EPSG代码
resample_method: 重采样方法 ('nearest', 'bilinear', 'cubic', 'cubicspline', 'lanczos')
"""
# 检查输入文件
if not os.path.exists(input_tif):
raise Exception(f"输入文件不存在：{input_tif}")
# 获取输入影像信息
print("正在获取输入影像信息...")
info = get_image_info(input_tif)
print(f"输入影像信息：")
print(f" 尺寸: {info['width']} x {info['height']}")
print(f" 波段数: {info['bands']}")
print(f" 当前坐标系: {info['proj_name']} (EPSG:{info['epsg']})")
# 设置重采样方法映射
resample_dict = {
'nearest': gdal.GRA_NearestNeighbour,
'bilinear': gdal.GRA_Bilinear,
'cubic': gdal.GRA_Cubic,
'cubicspline': gdal.GRA_CubicSpline,
'lanczos': gdal.GRA_Lanczos
}
resample_alg = resample_dict.get(resample_method.lower(), gdal.GRA_Bilinear)
# 创建目标坐标系
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(target_epsg)
print(f"目标坐标系: {target_srs.GetAttrValue('PROJCS') or target_srs.GetAttrValue('GEOGCS')} (EPSG:{target_epsg})")
print("正在执行坐标转换...")
# 使用GDAL Warp进行重投影
warp_options = gdal.WarpOptions(
dstSRS=f'EPSG:{target_epsg}',
resampleAlg=resample_alg,
format='GTiff',
creationOptions=['COMPRESS=LZW', 'TILED=YES'] # 添加压缩和分块选项
)
# 执行重投影
result = gdal.Warp(output_tif, input_tif, options=warp_options)
if result is None:
raise Exception("坐标转换失败！")
result = None
print(f"坐标转换完成！输出文件：{output_tif}")
# 验证输出文件
output_info = get_image_info(output_tif)
print(f"输出影像信息：")
print(f" 尺寸: {output_info['width']} x {output_info['height']}")
print(f" 坐标系: {output_info['proj_name']} (EPSG:{output_info['epsg']})")
def reproject_tif_advanced(input_tif, output_tif, target_epsg, pixel_size=None, extent=None):
"""
高级坐标转换方法，支持自定义像素大小和范围
参数:
input_tif: 输入TIF文件路径
output_tif: 输出TIF文件路径
target_epsg: 目标坐标系的EPSG代码
pixel_size: 输出像素大小 (xres, yres)，如果为None则自动计算
extent: 输出范围 (minx, miny, maxx, maxy)，如果为None则自动计算
"""
# 检查输入文件
if not os.path.exists(input_tif):
raise Exception(f"输入文件不存在：{input_tif}")
print("正在执行高级坐标转换...")
# 构建warp选项
warp_options_dict = {
'dstSRS': f'EPSG:{target_epsg}',
'resampleAlg': gdal.GRA_Bilinear,
'format': 'GTiff',
'creationOptions': ['COMPRESS=LZW', 'TILED=YES']
}
# 设置像素大小
if pixel_size:
warp_options_dict['xRes'] = pixel_size[0]
warp_options_dict['yRes'] = pixel_size[1]
# 设置输出范围
if extent:
warp_options_dict['outputBounds'] = extent
warp_options = gdal.WarpOptions(**warp_options_dict)
# 执行重投影
result = gdal.Warp(output_tif, input_tif, options=warp_options)
if result is None:
raise Exception("高级坐标转换失败！")
result = None
print(f"高级坐标转换完成！输出文件：{output_tif}")
def get_suitable_cgcs2000_epsg(lon, lat):
"""
根据经纬度获取合适的CGCS2000投影坐标系EPSG代码
参数:
lon: 经度
lat: 纬度
返回:
int: 合适的EPSG代码
"""
# CGCS2000 3度分带投影坐标系
# 中央经线从75°开始，每3度一个分带
zone = int((lon + 1.5) / 3) + 25
if zone < 25 or zone > 45:
# 如果超出范围，使用6度分带
zone_6 = int((lon + 6) / 6) + 13
epsg = 4490 + zone_6 # CGCS2000 6度分带
print(f"使用CGCS2000 6度分带，带号：{zone_6}，EPSG：{epsg}")
else:
epsg = 4513 + zone - 25 # CGCS2000 3度分带
print(f"使用CGCS2000 3度分带，带号：{zone}，EPSG：{epsg}")
return epsg
def auto_reproject_to_cgcs2000(input_tif, output_tif=None):
"""
自动选择合适的CGCS2000坐标系进行转换
参数:
input_tif: 输入TIF文件路径
output_tif: 输出TIF文件路径，如果为None则自动生成
"""
if output_tif is None:
base_name = os.path.splitext(input_tif)[0]
output_tif = f"{base_name}_cgcs2000.tif"
# 获取影像中心点坐标
dataset = gdal.Open(input_tif)
if dataset is None:
raise Exception(f"无法打开输入影像文件：{input_tif}")
geotransform = dataset.GetGeoTransform()
width = dataset.RasterXSize
height = dataset.RasterYSize
# 计算影像中心点
center_x = geotransform[0] + width * geotransform[1] / 2
center_y = geotransform[3] + height * geotransform[5] / 2
dataset = None
print(f"影像中心点坐标：({center_x:.6f}, {center_y:.6f})")
# 获取合适的CGCS2000 EPSG代码
target_epsg = get_suitable_cgcs2000_epsg(center_x, center_y)
# 执行转换
reproject_tif_simple(input_tif, output_tif, target_epsg)
return output_tif
if __name__ == "__main__":
# 使用示例
input_tif = r"D:\temp\rgb_uint8.tif" # 输入TIF文件路径
# 方法1：自动选择合适的CGCS2000坐标系（推荐）
print("=== 方法1：自动转换到CGCS2000 ===")
try:
output_file = auto_reproject_to_cgcs2000(input_tif)
print(f"自动转换完成！输出文件：{output_file}")
except Exception as e:
print(f"自动转换失败：{e}")
# 方法2：手动指定目标坐标系
print("\n=== 方法2：手动指定坐标系 ===")
output_tif = r"D:\temp\rgb_uint8_cgcs2000_manual.tif"
target_epsg = 4526 # CGCS2000 3度分带 120°E中央经线
try:
reproject_tif_simple(input_tif, output_tif, target_epsg)
print(f"手动转换完成！输出文件：{output_tif}")
except Exception as e:
print(f"手动转换失败：{e}")
# 方法3：高级转换（自定义参数）
print("\n=== 方法3：高级转换 ===")
output_tif_advanced = r"D:\temp\rgb_uint8_cgcs2000_advanced.tif"
try:
# 自定义像素大小为2米
reproject_tif_advanced(input_tif, output_tif_advanced, target_epsg, pixel_size=(2.0, 2.0))
print(f"高级转换完成！输出文件：{output_tif_advanced}")
except Exception as e:
print(f"高级转换失败：{e}")
print("\n=== 转换完成 ===")
print("注意事项：")
print("1. 确保输入文件路径正确")
print("2. CGCS2000坐标系会根据影像位置自动选择合适的分带")
print("3. 输出文件会自动压缩以节省空间")
print("4. 如果转换失败，请检查GDAL是否正确安装")
