**
使用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使用一张图片和三个角点的坐标进行地理配准，再定义投影完成坐标转换。
使用下面的代码获取仿射矩阵和投影参数：
|
1
2
3
4
5
6
|
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.批量获取图片的仿射矩阵
|
1
2
3
4
5
6
7
8
9
10
11
|
# 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
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
|
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!')
|
到此这篇关于使用Python和GDAL给图片加坐标系的实现思路(坐标投影转换)的文章就介绍到这了,更多相关Python GDAL坐标投影转换内容请搜索服务器之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持服务器之家！
原文链接：https://blog.csdn.net/qq_41440099/article/details/114793272