python GDAL 经纬度转像素坐标（包括投影坐标）

大江哥

1. # -*- encoding: utf-8 -*-
3. from osgeo import gdal
4. from osgeo import osr
5. import numpy as np
7. def getSRSPair(dataset):
8.     '''
9.     获得给定数据的投影参考系和地理参考系
10.     :param dataset: GDAL地理数据
11.     :return: 投影参考系和地理参考系
12.     '''
13.     prosrs = osr.SpatialReference()
14.     prosrs.ImportFromWkt(dataset.GetProjection())
15.     geosrs = prosrs.CloneGeogCS()
16.     return prosrs, geosrs
18. def geo2lonlat(dataset, x, y):
19.     '''
20.     将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
21.     :param dataset: GDAL地理数据
22.     :param x: 投影坐标x
23.     :param y: 投影坐标y
24.     :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
25.     '''
26.     prosrs, geosrs = getSRSPair(dataset)
27.     ct = osr.CoordinateTransformation(prosrs, geosrs)
28.     coords = ct.TransformPoint(x, y)
29.     return coords[:2]
32. def lonlat2geo(dataset, lon, lat):
33.     '''
34.     将经纬度坐标转为投影坐标（具体的投影坐标系由给定数据确定）
35.     :param dataset: GDAL地理数据
36.     :param lon: 地理坐标lon经度
37.     :param lat: 地理坐标lat纬度
38.     :return: 经纬度坐标(lon, lat)对应的投影坐标
39.     '''
40.     prosrs, geosrs = getSRSPair(dataset)
41.     ct = osr.CoordinateTransformation(geosrs, prosrs)
42.     coords = ct.TransformPoint(lon, lat)
43.     return coords[:2]
45. def imagexy2geo(dataset, row, col):
46.     '''
47.     根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
48.     :param dataset: GDAL地理数据
49.     :param row: 像素的行号
50.     :param col: 像素的列号
51.     :return: 行列号(row, col)对应的投影坐标或地理坐标(x, y)
52.     '''
53.     trans = dataset.GetGeoTransform()
54.     px = trans[0] + col * trans[1] + row * trans[2]
55.     py = trans[3] + col * trans[4] + row * trans[5]
56.     return px, py
59. def geo2imagexy(dataset, x, y):
60.     '''
61.     根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
62.     :param dataset: GDAL地理数据
63.     :param x: 投影或地理坐标x
64.     :param y: 投影或地理坐标y
65.     :return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
66.     '''
67.     trans = dataset.GetGeoTransform()
68.     a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
69.     b = np.array([x - trans[0], y - trans[3]])
70.     return np.linalg.solve(a, b)  # 使用numpy的linalg.solve进行二元一次方程的求解
73. if __name__ == '__main__':
74.     gdal.AllRegister()
75.     dataset = gdal.Open(r"F:\2016\Data\Great Khingan\DEM\Projection\strm_6102_UTM.tif")
76.     print('数据投影：')
77.     print(dataset.GetProjection())
78.     print('数据的大小（行，列）：')
79.     print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))
81.     x = 464201
82.     y = 5818760
83.     lon = 122.47242
84.     lat = 52.51778
85.     row = 2399
86.     col = 3751
88.     print('投影坐标 -> 经纬度：')
89.     coords = geo2lonlat(dataset, x, y)
90.     print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))
91.     print('经纬度 -> 投影坐标：')
92.     coords = lonlat2geo(dataset, lon, lat)
93.     print('(%s, %s)->(%s, %s)' % (lon, lat, coords[0], coords[1]))
95.     print('像素坐标 -> 投影坐标：')
96.     coords = imagexy2geo(dataset, row, col)
97.     print('(%s, %s)->(%s, %s)' % (row, col, coords[0], coords[1]))
98.     print('投影坐标 -> 像素坐标：')
99.     coords = geo2imagexy(dataset, x, y)
100.     print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))

复制代码

 

来源:https://www.cnblogs.com/Swlymbcty/p/17943215
免责声明：由于采集信息均来自互联网，如果侵犯了您的权益，请联系我们【E-Mail:cb@itdo.tech】 我们会及时删除侵权内容，谢谢合作！