118 lines
4.0 KiB
Python
118 lines
4.0 KiB
Python
|
# -*- encoding: utf-8 -*-
|
|||
|
|
# code from https://blog.csdn.net/theonegis/article/details/54427906
|
|||
|
|
from osgeo import gdal
|
|||
|
|
from osgeo import osr
|
|||
|
|
import numpy as np
|
|||
|
|
|
|||
|
|
def getSRSPair(dataset):
|
|||
|
|
"""
|
|||
|
|
获得给定数据的投影参考系和地理参考系
|
|||
|
|
:param dataset: GDAL地理数据
|
|||
|
|
:return: 投影参考系和地理参考系
|
|||
|
|
"""
|
|||
|
|
prosrs = osr.SpatialReference()
|
|||
|
|
prosrs.ImportFromWkt(dataset.GetProjection())
|
|||
|
|
geosrs = prosrs.CloneGeogCS()
|
|||
|
|
return prosrs, geosrs
|
|||
|
|
|
|||
|
|
|
|||
|
|
def geo2lonlat(dataset, x, y):
|
|||
|
|
"""
|
|||
|
|
将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
|
|||
|
|
:param dataset: GDAL地理数据
|
|||
|
|
:param x: 投影坐标x
|
|||
|
|
:param y: 投影坐标y
|
|||
|
|
:return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
|
|||
|
|
"""
|
|||
|
|
prosrs, geosrs = getSRSPair(dataset)
|
|||
|
|
ct = osr.CoordinateTransformation(prosrs, geosrs)
|
|||
|
|
coords = ct.TransformPoint(x, y)
|
|||
|
|
return coords[:2]
|
|||
|
|
|
|||
|
|
|
|||
|
|
def lonlat2geo(dataset, lon, lat):
|
|||
|
|
"""
|
|||
|
|
将经纬度坐标转为投影坐标（具体的投影坐标系由给定数据确定）
|
|||
|
|
:param dataset: GDAL地理数据
|
|||
|
|
:param lon: 地理坐标lon经度
|
|||
|
|
:param lat: 地理坐标lat纬度
|
|||
|
|
:return: 经纬度坐标(lon, lat)对应的投影坐标
|
|||
|
|
"""
|
|||
|
|
prosrs, geosrs = getSRSPair(dataset)
|
|||
|
|
ct = osr.CoordinateTransformation(geosrs, prosrs)
|
|||
|
|
coords = ct.TransformPoint(lat, lon)
|
|||
|
|
return coords[:2]
|
|||
|
|
|
|||
|
|
|
|||
|
|
def imagexy2geo(dataset, row, col):
|
|||
|
|
"""
|
|||
|
|
根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
|
|||
|
|
:param dataset: GDAL地理数据
|
|||
|
|
:param row: 像素的行号
|
|||
|
|
:param col: 像素的列号
|
|||
|
|
:return: 行列号(row, col)对应的投影坐标或地理坐标(x, y)
|
|||
|
|
"""
|
|||
|
|
trans = dataset.GetGeoTransform()
|
|||
|
|
px = trans[0] + col * trans[1] + row * trans[2]
|
|||
|
|
py = trans[3] + col * trans[4] + row * trans[5]
|
|||
|
|
return px, py
|
|||
|
|
|
|||
|
|
|
|||
|
|
def geo2imagexy(dataset, x, y):
|
|||
|
|
"""
|
|||
|
|
根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
|
|||
|
|
:param dataset: GDAL地理数据
|
|||
|
|
:param x: 投影或地理坐标x
|
|||
|
|
:param y: 投影或地理坐标y
|
|||
|
|
:return: 影坐标或地理坐标(x, y)对应的影像图上行列号(col, row)
|
|||
|
|
"""
|
|||
|
|
trans = dataset.GetGeoTransform()
|
|||
|
|
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
|
|||
|
|
b = np.array([x - trans[0], y - trans[3]])
|
|||
|
|
return np.linalg.solve(a, b) # 使用numpy的linalg.solve进行二元一次方程的求解
|
|||
|
|
|
|||
|
|
|
|||
|
|
def test1():
|
|||
|
|
gdal.AllRegister()
|
|||
|
|
tif = 'D:/DATA/testdata/GLCFCS30_E110N25.tif'
|
|||
|
|
# dataset = gdal.Open(r"D:\\DATA\\雷达测试\\GaoFen3_20200528_HH_DB.tif")
|
|||
|
|
dataset = gdal.Open(tif)
|
|||
|
|
|
|||
|
|
print('数据投影：')
|
|||
|
|
print(dataset.GetProjection())
|
|||
|
|
print('数据的大小（行，列）：')
|
|||
|
|
print('(%s %s)' % (dataset.RasterYSize, dataset.RasterXSize))
|
|||
|
|
|
|||
|
|
x = 793214.118
|
|||
|
|
y = 2485865.527
|
|||
|
|
lon = 113.84897082317516
|
|||
|
|
lat = 22.453998686022448
|
|||
|
|
row = 24576
|
|||
|
|
col = 22540
|
|||
|
|
|
|||
|
|
print('图上坐标 -> 投影坐标：')
|
|||
|
|
coords = imagexy2geo(dataset, row, col)
|
|||
|
|
print('(%s, %s)->(%s, %s)' % (row, col, coords[0], coords[1]))
|
|||
|
|
print('投影坐标 -> 图上坐标：')
|
|||
|
|
coords = geo2imagexy(dataset, x, y)
|
|||
|
|
col = coords[0]
|
|||
|
|
row = coords[1]
|
|||
|
|
print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))
|
|||
|
|
|
|||
|
|
print('投影坐标 -> 经纬度：')
|
|||
|
|
coords = geo2lonlat(dataset, x, y)
|
|||
|
|
print('(%s, %s)->(%s, %s)' % (x, y, coords[0], coords[1]))
|
|||
|
|
print('经纬度 -> 投影坐标：')
|
|||
|
|
coords = lonlat2geo(dataset, lon, lat)
|
|||
|
|
print('(%s, %s)->(%s, %s)' % (lon, lat, coords[0], coords[1]))
|
|||
|
|
|
|||
|
|
coords1 = geo2lonlat(dataset, 657974.118, 2633321.527)
|
|||
|
|
print(coords1)
|
|||
|
|
coords2 = geo2lonlat(dataset, 793214.118, 2485865.527)
|
|||
|
|
print(coords2)
|
|||
|
|
pass
|
|||
|
|
|
|||
|
|
# if __name__ == '__main__':
|
|||
|
|
#
|
|||
|
|
# print('done')
|