python gdal 影像投影坐标系转对应的地理坐标系

	image1 = gdal.Open("qian2_cgcs2000_zone39.tif")  # CGCS2000_3_Degree_GK_Zone_39
    # image2 = gdal.Open("/home/algo/b/zhangmin/code/hou1.tif")  # GCS_WGS_1984

    geotransform1 = image1.GetGeoTransform()

    projection1 = image1.GetProjection()

    srs1 = osr.SpatialReference()
    srs1.ImportFromWkt(projection1)

    # CloneGeoCS()函数可以获取投影坐标系对应的地理坐标系
    geocs1 = srs1.CloneGeogCS()
    # 坐标转换
    cor_tran = osr.CoordinateTransformation(srs1, geocs1)
    # 转换原点左上角坐标点
    coords = cor_tran.TransformPoint(geotransform1[0], geotransform1[3])
    # print(coords) # (116.13803729817306, 35.65263915023701, 0.0)

    # 计算右下角坐标
    right = geotransform1[0] + geotransform1[1] * image1.RasterXSize
    bottom = geotransform1[3] + geotransform1[5] * image1.RasterYSize

    # 转换右下角坐标点
    coords1 = cor_tran.TransformPoint(right, bottom)

    # 保持行列数不变的情况下，重新计算空间分辨率
    resolutionX = (coords1[0] - coords[0]) / image1.RasterXSize
    resolutionY = (coords1[1] - coords[1]) / image1.RasterYSize
    print(resolutionX, resolutionY)

    output_ds = gdal.GetDriverByName("GTiff").Create("qian_proj2geo2.tif",
                                                     image1.RasterXSize,  # 列数
                                                     image1.RasterYSize,  # 行数
                                                     image1.RasterCount,
                                                     image1.GetRasterBand(1).DataType)

    output_ds.SetGeoTransform((coords[0], resolutionX, 0, coords[1], 0, resolutionY))

    # 将空间参考对象（SpatialReference）转换为Well - KnownText（WKT）格式的字符串
    dstSRS = geocs1.ExportToWkt()
    output_ds.SetProjection(dstSRS)

    # 读取和写入重叠区域的数据
    for i in range(1, 2):
        image1_band = image1.GetRasterBand(i)
        data1 = image1_band.ReadAsArray()
        output_band = output_ds.GetRasterBand(i)
        output_band.WriteArray(data1)