使用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!')到此这篇关于使用Python和GDAL给图片加坐标系的实现思路(坐标投影转换)的文章就介绍到这了,更多相关Python GDAL坐标投影转换内容请搜索编程网以前的文章或继续浏览下面的相关文章希望大家以后多多支持编程网！
免责声明：
① 本站未注明“稿件来源”的信息均来自网络整理。其文字、图片和音视频稿件的所属权归原作者所有。本站收集整理出于非商业性的教育和科研之目的，并不意味着本站赞同其观点或证实其内容的真实性。仅作为临时的测试数据，供内部测试之用。本站并未授权任何人以任何方式主动获取本站任何信息。
② 本站未注明“稿件来源”的临时测试数据将在测试完成后最终做删除处理。有问题或投稿请发送至: 邮箱/279061341@qq.com QQ/279061341
