个人记录：TIF文件内部坐标系wgs84转gcj02

第一步安装Anaconda 这里就不赘述了

第二步创建环境

创建python环境，指定版本号
conda create --name test python=3.12.3 test指的是环境名，python指的是当前python的系统版本
激活python环境
activate test
安装gdal
conda install -c conda-forge gdal

第三步复制代码

import math
from osgeo import gdal, osr

# GCJ-02 坐标系转换算法（由 WGS-84 转为 GCJ-02）
def wgs84_to_gcj02(lng, lat):
    a = 6378245.0
    ee = 0.00669342162296594323
    pi = math.pi

    def transform_lat(x, y):
        ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * pi) + 20.0 * math.sin(2.0 * x * pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(y * pi) + 40.0 * math.sin(y / 3.0 * pi)) * 2.0 / 3.0
        ret += (160.0 * math.sin(y / 12.0 * pi) + 320 * math.sin(y * pi / 30.0)) * 2.0 / 3.0
        return ret

    def transform_lng(x, y):
        ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * math.sqrt(abs(x))
        ret += (20.0 * math.sin(6.0 * x * pi) + 20.0 * math.sin(2.0 * x * pi)) * 2.0 / 3.0
        ret += (20.0 * math.sin(x * pi) + 40.0 * math.sin(x / 3.0 * pi)) * 2.0 / 3.0
        ret += (150.0 * math.sin(x / 12.0 * pi) + 300.0 * math.sin(x / 30.0 * pi)) * 2.0 / 3.0
        return ret

    def out_of_china(lng, lat):
        return not (73.66 < lng < 135.05 and 3.86 < lat < 53.55)

    if out_of_china(lng, lat):
        return lng, lat

    dlat = transform_lat(lng - 105.0, lat - 35.0)
    dlng = transform_lng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return mglng, mglat

# 使用 GDAL 读取和转换 GeoTIFF 坐标系
def convert_tif_to_gcj02(input_tif, output_tif):
    # 读取 .tif 文件
    dataset = gdal.Open(input_tif)
    if dataset is None:
        print("无法打开文件")
        return

    # 获取原始投影和地理变换
    proj_wkt = dataset.GetProjection()
    geotransform = dataset.GetGeoTransform()

    # 原始坐标系（假设为 WGS84）
    source = osr.SpatialReference()
    source.ImportFromWkt(proj_wkt)

    # 创建 GCJ-02 的自定义投影 (以 WGS-84 为基础)
    target = osr.SpatialReference()
    target.ImportFromEPSG(4326)  # EPSG 4326 为 WGS84

    # 获取图像的宽度、高度
    width = dataset.RasterXSize
    height = dataset.RasterYSize

    # 获取原始图像的地理范围
    minx = geotransform[0]
    maxy = geotransform[3]
    maxx = minx + width * geotransform[1]
    miny = maxy + height * geotransform[5]

    # 转换左上角和右下角的坐标
    min_gcj02 = wgs84_to_gcj02(minx, miny)
    max_gcj02 = wgs84_to_gcj02(maxx, maxy)

    # 设置新的地理变换信息
    new_geotransform = (min_gcj02[0], geotransform[1], geotransform[2], max_gcj02[1], geotransform[4], geotransform[5])

    # 创建新的 TIFF 文件
    driver = gdal.GetDriverByName('GTiff')
    new_dataset = driver.CreateCopy(output_tif, dataset)

    # 设置新的投影和地理变换
    new_dataset.SetGeoTransform(new_geotransform)
    new_dataset.SetProjection(target.ExportToWkt())

    # 关闭文件
    new_dataset.FlushCache()
    new_dataset = None
    dataset = None

    print("转换完成，保存为: " + output_tif)

# 示例：转换文件
input_tif = r'D:\BaiduNetdiskDownload\123.tif'
output_tif = r'D:\BaiduNetdiskDownload\123_gcj03.tif'
convert_tif_to_gcj02(input_tif, output_tif)