使用Python和GDAL给图片加坐标系的实现思路(坐标投影转换)

**

使用Python和GDAL给图片加坐标系

**
假设你已经知道arcgis地理配准（如下图内容），懂一点python。

**

-目的和背景

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：图像南北方向分辨率；
‘''

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!')

来源：https://blog.csdn.net/qq_41440099/article/details/114793272