使用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_Byteelif 'int16' in img_array.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(img_array.shape) == 3:img_bands, im_height, im_width = img_array.shapeelse:img_bands, (im_height, im_width) = 1, img_array.shapefilename = 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.第一次发文章，如有侵权，欢迎沟通。

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！