使用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["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",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)
forkinrange(image_num):
ifimage_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["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",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
defdef_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'inimg_array.dtype.name:
datatype=gdal.GDT_Byte
elif'int16'inimg_array.dtype.name:
datatype=gdal.GDT_UInt16
else:
datatype=gdal.GDT_Float32
iflen(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)#写入投影
#写入影像数据
ifimg_bands==1:
dataset.GetRasterBand(1).WriteArray(img_array)
else:
foriinrange(img_bands):
dataset.GetRasterBand(i+1).WriteArray(img_array[i])
print(read_path,'geoCoordSysget!')
到此这篇关于使用Python和GDAL给图片加坐标系的实现思路(坐标投影转换)的文章就介绍到这了,更多相关PythonGDAL坐标投影转换内容请搜索毛票票以前的文章或继续浏览下面的相关文章希望大家以后多多支持毛票票！
声明：本文内容来源于网络，版权归原作者所有，内容由互联网用户自发贡献自行上传，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任。如果您发现有涉嫌版权的内容，欢迎发送邮件至：czq8825#qq.com（发邮件时，请将#更换为@）进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。