python将tif从wgs84转gcj02

news/2024/10/17 23:31:58/

# mark: 主要是dem从wgs84转成gcj02，步骤为wgs84的4326转成3857，之后投影进行偏移，再将偏移后的3857转成4326from 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 中的偏移量是各地根据实际计算得到的。