**
使用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坐標投影轉換內容請搜索腳本之家以前的文章或繼續瀏覽下面的相關文章希望大家以后多多支持腳本之家！
您可能感興趣的文章:- 解決python gdal投影坐標系轉換的問題
- python安裝gdal的兩種方法
- 利用python GDAL庫讀寫geotiff格式的遙感影像方法
- 在python中利用GDAL對tif文件進行讀寫的方法
- Python的地形三維可視化Matplotlib和gdal使用實例
