python将tif从wgs84转gcj02

dragonzoebai

已于 2023-06-14 16:25:48 修改

阅读量1.2k

点赞数

分类专栏： GIS 文章标签： python gdal tif wgs84 gcj02

于 2023-06-14 16:20:02 首次发布

本文链接：https://blog.csdn.net/dragonzoebai/article/details/131210433

版权

GIS 专栏收录该内容

25 篇文章 0 订阅

订阅专栏

# mark: 主要是dem从wgs84转成gcj02，步骤为wgs84的4326转成3857，之后投影进行偏移，再将偏移后的3857转成4326

from osgeo import gdal, osr

# 经纬度转投影
def tif4326To3857(input_file, output_file):
    options = gdal.WarpOptions(format='GTiff', srcSRS='EPSG:4326', dstSRS='EPSG:3857', width=1452)
    gdal.Warp(output_file, input_file, options=options)
    print('转换完成')

# 投影转经纬度
def tif3857To4326(input_file, output_file):
    options = gdal.WarpOptions(format='GTiff', srcSRS='EPSG:3857', dstSRS='EPSG:4326', width=1452)
    gdal.Warp(output_file, input_file, options=options)
    print('转换完成')

# wgs84投影偏移至gcj02
def tifWgs84ToGcj02(input_file, output_file):
    dataset = gdal.Open(input_file)
    im_width = dataset.RasterXSize #栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的列数
    im_array = dataset.GetRasterBand(1).ReadAsArray()
    im_projection = dataset.GetProjection()
    im_GeoTransform = dataset.GetGeoTransform()
    im_GeoTransform_new = (im_GeoTransform[0]+740.4591671, im_GeoTransform[1], im_GeoTransform[2], im_GeoTransform[3]+24.60901496, im_GeoTransform[4], im_GeoTransform[5])

    # 创建新的tif文件
    driver = gdal.GetDriverByName('GTiff')
    out_tif = driver.Create(output_file, im_width, im_height, 1, gdal.GDT_Float32)
    out_tif.SetGeoTransform(im_GeoTransform_new)
    out_tif.SetProjection(im_projection)  # 给新建图层赋予投影信息
    # # 获取地理坐标系统信息，用于选取需要的地理坐标系统
    # srs = osr.SpatialReference()
    # # 导入tif文件的坐标系
    # srs.ImportFromEPSG(3857)  # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
    # out_tif.SetProjection(srs.ExportToWkt())  # 给新建图层赋予投影信息

    # 数据写出
    out_tif.GetRasterBand(1).WriteArray(im_array)  # 将数据写入内存，此时没有写入硬盘
    out_tif.FlushCache()  # 将数据写入硬盘
    out_tif = None  # 注意必须关闭tif文件
    print('偏移完成')


if __name__ == '__main__':
    file1 = 'C:\\Users\\ZJZN\\Desktop\\tiff\\4326.tif'
    file2 = 'C:\\Users\\ZJZN\\Desktop\\tiff\\3857.tiff'
    file3 = 'C:\\Users\\ZJZN\\Desktop\\tiff\\3857-gcj02.tiff'
    file4 = 'C:\\Users\\ZJZN\\Desktop\\tiff\\4326-gcj02.tiff'
    tif4326To3857(file1, file2)
    tifWgs84ToGcj02(file2, file3)
    tif3857To4326(file3, file4)

tif4326To3857(file1, file2) 是将wgs84的经纬度转成投影
tifWgs84ToGcj02(file2, file3) 是将投影的进行偏移
tif3857To4326(file3, file4) 再将偏移后的投影转成经纬度

注意：im_GeoTransform_new 中的偏移量是各地根据实际计算得到的。