使用Python和GDAL给图片加其原有坐标系，进行坐标投影

原理：通过arcgis软件对一张图片进行地理配准，配准后通过GDAL的方法获取仿射矩阵参数，其他图片根据仿射矩阵参数进行投影。（PS：手动配准的图片，仿射矩阵参数精度有限，可能造成图片拉伸扭曲！）
**
-背景和目的
1.从地图网站获得一张PNG格式的图片，已知坐标系（WGS84）和左上角坐标。arcgis地理配准再定义投影即可给它加上原图的坐标系。
2.假设有上千张图片，可用Python和GDAL给图片加坐标系。

-实现思路

1.使用GDAL需要知道待投影图片的地理坐标信息、仿射矩阵参数。

仿射矩阵参数是干什么的？见：https://zhuanlan.zhihu.com/p/72184440
主要含义：
1）不同坐标系的转换，涉及到仿射变换，又称仿射映射，是指在几何中，一个向量空间进行一次线性变换 并接上一个平移，变换为另一个向量空间。
2）仿射矩阵信息有六个参数，描述的是栅格行列号和地理坐标之间的关系：
‘’’
0：图像左上角的X坐标；
1：图像东西方向分辨率；
2：旋转角度，如果图像北方朝上，该值为0；
3：图像左上角的Y坐标；
4：旋转角度，如果图像北方朝上，该值为0；
5：图像南北方向分辨率；
‘’’
由此可见，同一批图片只有左上角XY坐标不同，其他相同

2.在arcgis使用一张图片和三个角点的坐标进行地理配准，再定义投影完成坐标转换；

再使用下面的代码获取仿射矩阵和投影参数：

 dataset = gdal.Open('a.png')
 print （dataset.GetGeoTransform()）#仿射矩阵
 print （dataset.GetProjection()）#地图投影信息
# 打印结果为：
# (116.33333, 8.321688443e-05, 0.0, 39.99999, 0.0, -6.223016769e-05)
# 'GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG", "6326"]], PRIMEM["Greenwich", 0], UNIT["degree", 0.0174532925199433], AUTHORITY["EPSG", "4326"]]'

3.批量获取图片的仿射矩阵

# coors是用来存储图片对应左上角坐标的字典。格式为{‘a.png‘’：[116.33333,39.6]，}
    image_list = os.listdir('D:\\dd')
    image_num = len(image_list)
    for k in range(image_num):
        if image_list[k].endswith('.png'):
            img_name = img_none_path + '/' + image_list[k]
            img_pos_transf = (float(coors[image_list[k]][0]), 8.321688443e-05,
                              0.0, float(coors[image_list[k]][1]), 0.0, -6.223016769e-05)#根据第二步获得像元分辨率和投影
            print(img_pos_transf)
            img_pos_proj = 'GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG", "6326"]], PRIMEM["Greenwich", 0], UNIT["degree", 0.0174532925199433], AUTHORITY["EPSG", "4326"]]'
            def_geoCoordSys(img_name, img_pos_transf, img_pos_proj)#坐标转换的函数

4.给图片加坐标系的主要函数如下

来自文章 ：https://blog.csdn.net/nominior/article/details/102737294

def def_geoCoordSys(read_path, img_transf, img_proj):
        array_dataset = gdal.Open(read_path)
        img_array = array_dataset.ReadAsArray(
            0, 0, array_dataset.RasterXSize, array_dataset.RasterYSize)
        if 'int8' in img_array.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in img_array.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32

        if len(img_array.shape) == 3:
            img_bands, im_height, im_width = img_array.shape
        else:
            img_bands, (im_height, im_width) = 1, img_array.shape

        filename = read_path[:-4] + '_proj' + '.tif'
        driver = gdal.GetDriverByName("GTiff")  # 创建文件驱动
        dataset = driver.Create(
            filename, im_width, im_height, img_bands, datatype)
        dataset.SetGeoTransform(img_transf)  # 写入仿射变换参数
        dataset.SetProjection(img_proj)  # 写入投影

        # 写入影像数据
        if img_bands == 1:
            dataset.GetRasterBand(1).WriteArray(img_array)
        else:
            for i in range(img_bands):
                dataset.GetRasterBand(i + 1).WriteArray(img_array[i])
        print(read_path, 'geoCoordSys get!')

5.第一次发文章，如有侵权，欢迎沟通。